quad.form: Evaluate a quadratic form efficiently

In emulator: Bayesian Emulation of Computer Programs

Description

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.

(in the following, t(x) 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 t(x) %*% M %*% x in an efficient manner

• Function quad.form.inv(M,x) returns t(x) %*% solve(M) %*% x using an efficient method that avoids inverting M

• Function quad.tform(M,x) returns x %*% A %*% t(x) using tcrossprod() without taking a transpose

• Function quad.tform.inv(M,x) returns x %*% solve(M) %*% t(x), although a single transpose is needed

• Function quad.3form(M,l,r) returns t(l) %*% M %*% r 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 t(l) %*% solve(M) %*% r

• Function quad.3tform(M,l,r) returns l %*% M %*% t(r) 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 t(x) %*% y, equivalent to crossprod(Conj(x),y)

• Function tcprod(x,y) returns x %*% t(y), equivalent to crossprod(x,Conj(y))

Note again that in the calls above, “transpose” [that is, t(x)] 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*n
x,y	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

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*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

Example output

Loading required package: mvtnorm
          [,1]      [,2]
[1,]  4.799041 -2.026542
[2,] -2.026542 44.979133
          [,1]      [,2]
[1,]  4.799041 -2.026542
[2,] -2.026542 44.979133
          [,1]      [,2]
[1,]  4.799041 -2.026542
[2,] -2.026542 44.979133
 [1] 1.9752066+0i 1.3578512-0i 1.5768595+0i 1.3471074-0i 1.4214876-0i
 [6] 0.2471074+0i 0.5495868+0i 2.0752066+0i 1.3876033-0i 1.4983471-0i
[11] 1.2479339-0i
[1] 7.314050+0i 7.306612-0i
     [,1]         [,2]
[1,]    0 1.387779e-17
[2,]    0 0.000000e+00
              [,1]         [,2]
[1,]  0.000000e+00 3.552714e-15
[2,] -3.552714e-15 0.000000e+00
[1] 0 0
 [1]  0.000000e+00-5.551115e-17i -2.220446e-16+0.000000e+00i
 [3]  0.000000e+00+2.775558e-17i  0.000000e+00+4.163336e-17i
 [5]  0.000000e+00+0.000000e+00i  0.000000e+00+3.469447e-18i
 [7]  0.000000e+00-6.938894e-18i  0.000000e+00+5.551115e-17i
 [9]  0.000000e+00+5.551115e-17i  0.000000e+00+1.387779e-17i
[11]  0.000000e+00-2.428613e-17i
[1] -8.881784e-16-4.440892e-16i  8.881784e-16+3.191891e-16i

emulator documentation built on April 25, 2021, 9:07 a.m.