quad.form: Evaluate a quadratic form efficiently

Description Usage Arguments Details Note Author(s) See Also Examples

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).

These functions invoke the following lower-level calls:

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

Usage

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
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)
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

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
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

emulator documentation built on May 19, 2017, 10:26 p.m.

Search within the emulator package
Search all R packages, documentation and source code