quad.form: Evaluate a quadratic form efficiently
In quadform: Efficient Evaluation of Quadratic Forms
quad.form R Documentation
Evaluate a quadratic form efficiently
Description
\loadmathjax
Given a square matrix \mjseqnM of size \mjeqnn\times nn*n, and a
matrix \mjseqnx of size \mjeqnn\times pn*p (or a vector of length
\mjseqnn), evaluate various quadratic forms.
The archetype is quad.form(M,x) for real or complex square
matrix M and vector or matrix x. This evaluates
Conj(t(x)) %*% M %*% x but using
crossprod(crossprod(M,Conj(x)),x) thus avoiding taking a
needless transpose.
Usage
quad.form(M, x)
quad.form.inv(M, x)
quad.form.chol(chol,x)
quad.tform(M, x)
quad3.form(M,left,right)
quad3.tform(M,left,right)
quad.tform.inv(M,x)
quad.diag(M,x)
quad.tdiag(M,x)
quad3.diag(M,left,right)
quad3.tdiag(M,left,right)
cprod(x,y)
tcprod(x,y)
ht(x)
Arguments
M
Square matrix of size \mjeqnn\times nn*n
x, y
Matrix of size \mjeqnn\times pn*p, or vector of length \mjseqnn
chol
Lower triangular Cholesky decomposition
of the quadratic form, see details
left, right
In function quad3.form(), matrices with
\mjseqnn rows and arbitrary number of columns
Details
ht(x) \mjeqn\mathbfx^*=\overline\mathbfx^\scriptscriptstyle Tomitted
Conj(t(x))ht(x)
cprod(x,y) \mjeqn\mathbfx^*\mathbfyomitted
ht(x) %*% ycp()
tcprod(x,y) \mjeqn\mathbfx\mathbfy^*omitted
x %*% ht(y)tcp()
quad.form(M,x) \mjeqn\mathbfx^*M\mathbfxomitted
ht(x) %*% M %*% xqf()
quad.form.inv(M,x) \mjeqn\mathbfx^*M^-1\mathbfxomitted
ht(x) %*% solve(M) %*% xqfi()
quad.tform(M,x) \mjeqn\mathbfxM\mathbfx^*omitted
x %*% A %*% ht(x)qt()
quad.tform.inv(M,x) \mjeqn\mathbfxM^-1\mathbfx^*omitted
x %*% solve(M) %*% ht(x)qti()
quad3.form(M,l,r) \mjeqn\mathbfl^*M\mathbfromitted
t(l) %*% M %*% rq3()
quad3.form.inv(M,l,r) \mjeqn\mathbfl^*M^-1\mathbfromitted
t(l) %*% solve(M) %*% rq3i()
quad3.tform(M,l,r) \mjeqn\mathbflM\mathbfr^*omitted
l %*% M %*% t(r)q3t()
quad.diag(M,x) \mjeqn\operatornamediag(\mathbfx^*M\mathbfx)omitted
diag(quad.form(M,x))qd()
quad.tdiag(M,x) \mjeqn\operatornamediag(\mathbfxM\mathbfx^*)omitted
diag(quad.tform(M,x))qtd()
quad3.diag(M,l,r) \mjeqn\operatornamediag(\mathbfl^*M\mathbfr)omitted
diag(quad3.form(M,l,r))q3d()
quad3.tdiag(M,l,r) \mjeqn\operatornamediag(\mathbflM\mathbfr^*)omitted
diag(quad3.tform(M,l,r))q3td()
quad.trace(M,x) \mjeqn\operatornametr(\mathbfx^*M\mathbfx)omitted
tr(quad.form(M,x))qt()
quad.ttrace(M,x) \mjeqn\operatornametr(\mathbfxM\mathbfx^*)omitted
tr(quad.tform(M,x))qtt()
In the above, \mjeqn\mathbfx^*ht(x) denotes the
complex conjugate of the transpose, also known as the Hermitian
transpose (this only matters when considering complex numbers).
These various functions generally avoid taking needless expensive
transposes in favour of using nested crossprod() and
tcrossprod() calls. For example, the “meat” of
quad.form() is just crossprod(crossprod(M,Conj(x)),x).
Functions such as quad.form.inv() avoid taking a matrix
inverse. The meat of quad.form.inv(), for example, is
cprod(x, solve(M, x)). Many people have stated things like
“Never invert a matrix unless absolutely necessary”. But I
have never seen a case where quad.form.inv(M,x) is
faster than quad.form(solve(M),x).
One motivation for the package is to return consistent results with
complex arguments. Note, for example, that
base::crossprod(x,y) evaluates t(x) %*% y and not, as
one would almost always want, Conj(t(x)) %*% y. Function
cprod(), unlike crossprod(), is consistent and returns
Conj(t(x)) %*% y [or ht(x) %*% y]; internally it is
essentially crossprod(Conj(x), y).
Function quad.form.chol() interprets argument chol 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)) in calls. If the Cholesky decomposition of M
is available, then using quad.form.chol() and supplying
chol(M) should generally be faster (for large matrices) than
calling quad.form() and using M directly. The time
saving is negligible for matrices smaller than about 50\times
50, even if the overhead of computing the decomposition is
ignored.
Functions quad3.foo() take three arguments: a matrix M
and two other vectors l and r [or left and
right]. For these functions, M is not necessarily
square although of course the matrices have to be compatible.
Functions quad3.form_ab() and quad3.form_bc() are helper
functions not really intended for the end-user. They return
mathematically identical results but differ in the bracketing order of
their operations: quad3.form_ab(M,l,r) returns
\mjeqn\left(\mathbfl^*M\right)\mathbfromitted and
quad3.form_bc(M,l,r) returns
\mjeqn\mathbfl^*\left(M\mathbfr\right)omitted. The mnemonic
for their names is derived from the first multiplication when
calculating \mjseqn(ab)c and \mjseqna(bc). Note that
quad3.form_ab(M,l,r) returns
crossprod(crossprod(M,Conj(l)),r) rather than the
mathematically equivalent cprod(cprod(M,l),r) on efficiency
grounds (only a single conjugate is taken).
Function quad3.form() dispatches to either
quad3.form_ab() or quad3.form_bc() depending on the
dimensions of its argument as per the efficiency discussion at
inst/quadform3test.Rmd. Similar considerations apply to
quad3.tform(), quad3.tform_ab(), and
quad3.tform_bc().
Terse forms [qf() for quad.form(), qti() for
quad.tform.inv(), etc] are provided for the perl
golfers among us.
Value
Generally, return a (dropped) matrix, real or complex as appropriate
Note
These functions are used extensively in the emulator and
calibrator packages, 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
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.chol(M.lower, x)
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 <- matrix(rnorm(22) + 1i*rnorm(22),2,11)
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
quadform documentation built on May 29, 2024, 1:26 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
| quad.form | R Documentation |
Evaluate a quadratic form efficiently
Description
\loadmathjaxGiven a square matrix \mjseqnM of size \mjeqnn\times nn*n, and a matrix \mjseqnx of size \mjeqnn\times pn*p (or a vector of length \mjseqnn), evaluate various quadratic forms.
The archetype is quad.form(M,x) for real or complex square
matrix M and vector or matrix x. This evaluates
Conj(t(x)) %*% M %*% x but using
crossprod(crossprod(M,Conj(x)),x) thus avoiding taking a
needless transpose.
Usage
quad.form(M, x)
quad.form.inv(M, x)
quad.form.chol(chol,x)
quad.tform(M, x)
quad3.form(M,left,right)
quad3.tform(M,left,right)
quad.tform.inv(M,x)
quad.diag(M,x)
quad.tdiag(M,x)
quad3.diag(M,left,right)
quad3.tdiag(M,left,right)
cprod(x,y)
tcprod(x,y)
ht(x)
Arguments
M |
Square matrix of size \mjeqnn\times nn*n |
x, y |
Matrix of size \mjeqnn\times pn*p, or vector of length \mjseqnn |
chol |
Lower triangular Cholesky decomposition of the quadratic form, see details |
left, right |
In function |
Details
ht(x) | \mjeqn\mathbfx^*=\overline\mathbfx^\scriptscriptstyle Tomitted | Conj(t(x)) | ht(x) |
cprod(x,y) | \mjeqn\mathbfx^*\mathbfyomitted | ht(x) %*% y | cp() |
tcprod(x,y) | \mjeqn\mathbfx\mathbfy^*omitted | x %*% ht(y) | tcp() |
quad.form(M,x) | \mjeqn\mathbfx^*M\mathbfxomitted | ht(x) %*% M %*% x | qf() |
quad.form.inv(M,x) | \mjeqn\mathbfx^*M^-1\mathbfxomitted | ht(x) %*% solve(M) %*% x | qfi() |
quad.tform(M,x) | \mjeqn\mathbfxM\mathbfx^*omitted | x %*% A %*% ht(x) | qt() |
quad.tform.inv(M,x) | \mjeqn\mathbfxM^-1\mathbfx^*omitted | x %*% solve(M) %*% ht(x) | qti() |
quad3.form(M,l,r) | \mjeqn\mathbfl^*M\mathbfromitted | t(l) %*% M %*% r | q3() |
quad3.form.inv(M,l,r) | \mjeqn\mathbfl^*M^-1\mathbfromitted | t(l) %*% solve(M) %*% r | q3i() |
quad3.tform(M,l,r) | \mjeqn\mathbflM\mathbfr^*omitted | l %*% M %*% t(r) | q3t() |
quad.diag(M,x) | \mjeqn\operatornamediag(\mathbfx^*M\mathbfx)omitted | diag(quad.form(M,x)) | qd() |
quad.tdiag(M,x) | \mjeqn\operatornamediag(\mathbfxM\mathbfx^*)omitted | diag(quad.tform(M,x)) | qtd() |
quad3.diag(M,l,r) | \mjeqn\operatornamediag(\mathbfl^*M\mathbfr)omitted | diag(quad3.form(M,l,r)) | q3d() |
quad3.tdiag(M,l,r) | \mjeqn\operatornamediag(\mathbflM\mathbfr^*)omitted | diag(quad3.tform(M,l,r)) | q3td() |
quad.trace(M,x) | \mjeqn\operatornametr(\mathbfx^*M\mathbfx)omitted | tr(quad.form(M,x)) | qt() |
quad.ttrace(M,x) | \mjeqn\operatornametr(\mathbfxM\mathbfx^*)omitted | tr(quad.tform(M,x)) | qtt()
|
In the above, \mjeqn\mathbfx^*ht(x) denotes the complex conjugate of the transpose, also known as the Hermitian transpose (this only matters when considering complex numbers).
These various functions generally avoid taking needless expensive
transposes in favour of using nested crossprod() and
tcrossprod() calls. For example, the “meat” of
quad.form() is just crossprod(crossprod(M,Conj(x)),x).
Functions such as quad.form.inv() avoid taking a matrix
inverse. The meat of quad.form.inv(), for example, is
cprod(x, solve(M, x)). Many people have stated things like
“Never invert a matrix unless absolutely necessary”. But I
have never seen a case where quad.form.inv(M,x) is
faster than quad.form(solve(M),x).
One motivation for the package is to return consistent results with
complex arguments. Note, for example, that
base::crossprod(x,y) evaluates t(x) %*% y and not, as
one would almost always want, Conj(t(x)) %*% y. Function
cprod(), unlike crossprod(), is consistent and returns
Conj(t(x)) %*% y [or ht(x) %*% y]; internally it is
essentially crossprod(Conj(x), y).
Function quad.form.chol() interprets argument chol 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)) in calls. If the Cholesky decomposition of M
is available, then using quad.form.chol() and supplying
chol(M) should generally be faster (for large matrices) than
calling quad.form() and using M directly. The time
saving is negligible for matrices smaller than about 50\times
50, even if the overhead of computing the decomposition is
ignored.
Functions quad3.foo() take three arguments: a matrix M
and two other vectors l and r [or left and
right]. For these functions, M is not necessarily
square although of course the matrices have to be compatible.
Functions quad3.form_ab() and quad3.form_bc() are helper
functions not really intended for the end-user. They return
mathematically identical results but differ in the bracketing order of
their operations: quad3.form_ab(M,l,r) returns
\mjeqn\left(\mathbfl^*M\right)\mathbfromitted and
quad3.form_bc(M,l,r) returns
\mjeqn\mathbfl^*\left(M\mathbfr\right)omitted. The mnemonic
for their names is derived from the first multiplication when
calculating \mjseqn(ab)c and \mjseqna(bc). Note that
quad3.form_ab(M,l,r) returns
crossprod(crossprod(M,Conj(l)),r) rather than the
mathematically equivalent cprod(cprod(M,l),r) on efficiency
grounds (only a single conjugate is taken).
Function quad3.form() dispatches to either
quad3.form_ab() or quad3.form_bc() depending on the
dimensions of its argument as per the efficiency discussion at
inst/quadform3test.Rmd. Similar considerations apply to
quad3.tform(), quad3.tform_ab(), and
quad3.tform_bc().
Terse forms [qf() for quad.form(), qti() for
quad.tform.inv(), etc] are provided for the perl
golfers among us.
Value
Generally, return a (dropped) matrix, real or complex as appropriate
Note
These functions are used extensively in the emulator and
calibrator packages, 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
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.chol(M.lower, x)
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 <- matrix(rnorm(22) + 1i*rnorm(22),2,11)
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
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.