quad.form {emulator}

Evaluate a quadratic form efficiently

Description

Given a square matrix M of size n\times n, and a matrix x of size n\times p (or a vector of length n), evaluate various quadratic forms.

(in the following, x^T denotes the complex conjugate of the transpose, also known as the Hermitian transpose. This only matters when considering complex numbers).

• Function quad.form(M,x) evaluates x^TMx in an efficient manner

• Function quad.form.inv(M,x) returns x^TM^{-1}x using an efficient method that avoids inverting M

• Function quad.tform(M,x) returns xMx^T using tcrossprod() without taking a transpose

• Function quad.tform.inv(M,x) returns xM^{-1}x^T, although a single transpose is needed

• Function quad.3form(M,l,r) returns l^TMr using nested calls to crossprod(). It's no faster than calling crossprod() directly, but makes code neater and less error-prone (IMHO)

• Function quad.3form.inv(M,l,r) returns l^TM^{-1}r

• Function quad.3tform(M,l,r) returns lMr^T using nested calls to tcrossprod(). Again, this is to make for neater code

• Function quad.diag(M,x) returns the diagonal of the (potentially very large) square matrix quad.form(M,x) without calculating the off diagonal elements

• Function quad.tdiag(M,x) similarly returns the diagonal of quad.tform(M,x)

• Function quad.3diag(M,l,r) returns the diagonal of quad.3form(M,l,r)

• Function quad.3tdiag(M,l,r) returns the diagonal of quad.3tform(M,l,r)

These functions invoke the following lower-level calls:

• Function ht(x) returns the Hermitian transpose, that is, the complex conjugate of the transpose, sometimes written x^*

• Function cprod(x,y) returns x^T y, equivalent to crossprod(Conj(x),y)

• Function tcprod(x,y) returns x y^T, equivalent to crossprod(x,Conj(y))

Note again that in the calls above, “transpose” [that is, x^T] means “Conjugate transpose”, or the Hermitian transpose.

Usage

quad.form(M, x, chol=FALSE)
quad.form.inv(M, x)
quad.tform(M, x)
quad.3form(M,left,right)
quad.3tform(M,left,right)
quad.tform.inv(M,x)
quad.diag(M,x)
quad.tdiag(M,x)
quad.3diag(M,left,right)
quad.3tdiag(M,left,right)
cprod(x,y)
tcprod(x,y)
ht(x)

Arguments

M

Square matrix of size n\times n

x, y

Matrix of size n\times 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

Details

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\times 50, even if the overhead of computing M.upper is ignored.

Note

These functions are used extensively in the emulator and calibrator packages' R code, primarily in the interests of elegant code, but also speed. For the problems I usually consider, the speedup (of quad.form(M,x) over t(x) %*% M %*% x, say) is marginal at best.

Author(s)

Robin K. S. Hankin

See Also

optimize

Examples

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)




## Make two Hermitian positive-definite matrices:
L <- matrix(c(1,0.1i,-0.1i,1),2,2)
LL <- diag(11)
LL[2,1] <- -(LL[1,2] <- 0.1i)


z <- t(latin.hypercube(11,2,complex=TRUE))


quad.diag(L,z)     # elements real because L is HPD
quad.tdiag(LL,z)   # ditto




## 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
quad.diag(M,x) - diag(quad.form(M,x))       # should be zero
diag(quad.form(L,z))   - quad.diag(L,z)     # should be zero
diag(quad.tform(LL,z)) - quad.tdiag(LL,z)   # should be zero