emulator (version 1.2-20)

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)

• Functionquad.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.

See Also

optimize

Examples

# NOT RUN {
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
# }

Run the code above in your browser using DataCamp Workspace