quad.form: 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

Author

Robin K. S. Hankin

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.

Examples

Run this code

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