quad.form {emulator} | R Documentation |
Given a square matrix M of size n*n, and a matrix x of size n*p (or a vector of length n), evaluate various quadratic forms.
quad.form()
evaluates t(x)
%*% M %*% x in an efficient manner
quad.form.inv()
returns t(x)
%*% solve(M) %*% x using an efficient method that avoids
inverting M
quad.3form()
returns t(l) %*'% M
%*% r using nested calls to crossprod()
. It's no faster
than calling crossprod()
directly, but makes code neater
(IMHO)
quad.tform()
returns x %*% A %*%
t(x) using tcrossprod()
in such a way as to not take a
transpose
quad.tform.inv()
returns x
%*% solve(M) %*% t(x), although a single transpose is needed
quad.form(M, x, chol=FALSE) quad.form.inv(M, x) quad.tform(M, x) quad.3form(M,left,right) quad.tform.inv(M,x)
M |
Square matrix of size n*n |
x |
Matrix of size n*p, or vector of length n |
chol |
Boolean, with TRUE meaning to interpret
argument M as the lower triangular Cholesky decomposition
of the quadratic form. Remember that
M.lower %*% M.upper == M , and chol() returns the
upper triangular matrix, so one needs to use the transpose
t(chol(M)) |
left,right |
In function quad.3form() , matrices with
n rows and arbitrary number of columns |
The “meat” of quad.form()
for chol=FALSE
is just
crossprod(crossprod(M, x), x)
, and that of
quad.form.inv()
is crossprod(x, solve(M, x))
.
If the Cholesky decomposition of M
is available, then calling
with chol=TRUE
and supplying M.upper
should generally be
faster (for large matrices) than calling with chol=FALSE
and
using M
directly. The time saving is negligible for matrices
smaller than about 50*50, even if the overhead of
computing M.upper
is ignored.
These functions are used extensively in the emulator and calibrator packages' R code in the interests of speed. They are not really intended for the end user.
Robin K. S. Hankin
jj <- matrix(rnorm(80),20,4) M <- crossprod(jj,jj) M.lower <- t(chol(M)) x <- matrix(rnorm(8),4,2) jj.1 <- t(x) %*% M %*% x jj.2 <- quad.form(M,x) jj.3 <- quad.form(M.lower,x,chol=TRUE) print(jj.1) print(jj.2) print(jj.3) ## Now consider accuracy: quad.form(solve(M),x) - quad.form.inv(M,x) # should be zero quad.form(M,x) - quad.tform(M,t(x)) # Should be zero