Title: | Fast Matrix Operations via the Woodbury Matrix Identity |
---|---|
Description: | A hierarchy of classes and methods for manipulating matrices formed implicitly from the sums of the inverses of other matrices, a situation commonly encountered in spatial statistics and related fields. Enables easy use of the Woodbury matrix identity and the matrix determinant lemma to allow computation (e.g., solving linear systems) without having to form the actual matrix. More information on the underlying linear algebra can be found in Harville, D. A. (1997) <doi:10.1007/b98818>. |
Authors: | Michael Bertolacci [aut, cre, cph]
|
Maintainer: | Michael Bertolacci <[email protected]> |
License: | MIT + file LICENSE |
Version: | 0.0.3 |
Built: | 2025-02-08 04:22:38 UTC |
Source: | https://github.com/mbertolacci/woodburymatrix |
Calculates the (log) determinant of a WoodburyMatrix
using the matrix determinant lemma.
## S4 method for signature 'WoodburyMatrix,logical' determinant(x, logarithm)
## S4 method for signature 'WoodburyMatrix,logical' determinant(x, logarithm)
x |
A object that is a subclass of |
logarithm |
Logical indicating whether to return the logarithm of the matrix. |
Same as base::determinant.
This is a generic to represent direct instantiation of an implicitly defined matrix. In general this is a bad idea, because it removes the whole purpose of using an implicit representation.
instantiate(x) ## S4 method for signature 'GWoodburyMatrix' instantiate(x) ## S4 method for signature 'SWoodburyMatrix' instantiate(x)
instantiate(x) ## S4 method for signature 'GWoodburyMatrix' instantiate(x) ## S4 method for signature 'SWoodburyMatrix' instantiate(x)
x |
Implicit matrix to directly instantiate. |
The directly instantiated matrix.
instantiate,GWoodburyMatrix-method
: Method for general matrices.
instantiate,SWoodburyMatrix-method
: Method for symmetric matrices.
WoodburyMatrix, WoodburyMatrix
Generic for computing the squared Mahalanobis distance.
mahalanobis(x, center, cov, inverted = FALSE, ...) ## S4 method for signature 'ANY,ANY,SWoodburyMatrix' mahalanobis(x, center, cov, inverted = FALSE, ...)
mahalanobis(x, center, cov, inverted = FALSE, ...) ## S4 method for signature 'ANY,ANY,SWoodburyMatrix' mahalanobis(x, center, cov, inverted = FALSE, ...)
x |
Numeric vector or matrix |
center |
Numeric vector representing the mean; if omitted, defaults to zero mean |
cov |
Covariance matrix |
inverted |
Whether to treat |
... |
Passed to the |
x = ANY,center = ANY,cov = SWoodburyMatrix
: Use the Woodbury matrix identity to compute the
squared Mahalanobis distance with the implicit matrix as the covariance.
SWoodburyMatrix
objectsDraw samples and compute density functions for the multivariate normal
distribution with an SWoodburyMatrix
object as its covariance matrix.
dwnorm(x, mean, covariance, log = FALSE) rwnorm(n, mean, covariance)
dwnorm(x, mean, covariance, log = FALSE) rwnorm(n, mean, covariance)
x |
A numeric vector or matrix. |
mean |
Optional mean vector; defaults to zero mean. |
covariance |
|
log |
Logical indicating whether to return log of density. |
n |
Number of samples to return. If |
dwnorm
: Compute the density of the
distribution
rwnorm
: Draw samples from the distribution
library(Matrix) # Trivial example with diagonal covariance matrices W <- WoodburyMatrix(Diagonal(10), Diagonal(10)) x <- rwnorm(10, covariance = W) print(dwnorm(x, covariance = W, log = TRUE))
library(Matrix) # Trivial example with diagonal covariance matrices W <- WoodburyMatrix(Diagonal(10), Diagonal(10)) x <- rwnorm(10, covariance = W) print(dwnorm(x, covariance = W, log = TRUE))
WoodburyMatrix
objectsMethods based on solve
to solve a linear system of equations
involving WoodburyMatrix
objects. These methods take
advantage of the Woodbury matrix identity and therefore can be much more
time and memory efficient than forming the matrix directly.
Calling this function while omitting the b
argument returns the
inverse of a
. This is NOT recommended, since it removes any benefit
from using an implicit representation of a
.
## S4 method for signature 'GWoodburyMatrix,missing' solve(a) ## S4 method for signature 'GWoodburyMatrix,numLike' solve(a, b) ## S4 method for signature 'GWoodburyMatrix,matrix' solve(a, b) ## S4 method for signature 'GWoodburyMatrix,ANY' solve(a, b) ## S4 method for signature 'SWoodburyMatrix,missing' solve(a) ## S4 method for signature 'SWoodburyMatrix,numLike' solve(a, b) ## S4 method for signature 'SWoodburyMatrix,matrix' solve(a, b) ## S4 method for signature 'SWoodburyMatrix,ANY' solve(a, b)
## S4 method for signature 'GWoodburyMatrix,missing' solve(a) ## S4 method for signature 'GWoodburyMatrix,numLike' solve(a, b) ## S4 method for signature 'GWoodburyMatrix,matrix' solve(a, b) ## S4 method for signature 'GWoodburyMatrix,ANY' solve(a, b) ## S4 method for signature 'SWoodburyMatrix,missing' solve(a) ## S4 method for signature 'SWoodburyMatrix,numLike' solve(a, b) ## S4 method for signature 'SWoodburyMatrix,matrix' solve(a, b) ## S4 method for signature 'SWoodburyMatrix,ANY' solve(a, b)
a |
|
b |
Matrix, vector, or similar (needs to be compatible with the
submatrices |
The solution to the linear system, or the inverse of the matrix. The
class of the return value will be a vector if b
is a vector, and may
otherwise be either a regular matrix or a subclass of
Matrix
, with the specific subclass determined by
a
and b
.
solve,GWoodburyMatrix,missing-method
: Invert the matrix
solve,GWoodburyMatrix,numLike-method
: Solve the linear system
solve,GWoodburyMatrix,matrix-method
: Solve the linear system
solve,GWoodburyMatrix,ANY-method
: Solve the linear system
solve,SWoodburyMatrix,missing-method
: Invert the symmetric matrix
solve,SWoodburyMatrix,numLike-method
: Solve the linear system
solve,SWoodburyMatrix,matrix-method
: Solve the linear system
solve,SWoodburyMatrix,ANY-method
: Solve the linear system
WoodburyMatrix, WoodburyMatrix
Creates an implicitly defined matrix representing the equation
where and
are n x n, n x p, p x p and p x n matrices,
respectively. A symmetric special case is also possible with
where is n x p and
and
are additionally symmetric.
The available methods are described in WoodburyMatrix-class
and in solve. Multiple B / U / V / X matrices are also supported; see
below
WoodburyMatrix(A, B, U, V, X, O, symmetric)
WoodburyMatrix(A, B, U, V, X, O, symmetric)
A |
Matrix |
B |
Matrix |
U |
Matrix |
V |
Matrix |
X |
Matrix |
O |
Optional, precomputed value of |
symmetric |
Logical value, whether to create a symmetric or general matrix. See Details section for more information. |
The benefit of using an implicit representation is that the inverse of this matrix can be efficiently calculated via
where , and its determinant by
These relationships are often called the Woodbury matrix identity and the
matrix determinant lemma, respectively. If and
are sparse or
otherwise easy to deal with, and/or when
, manipulating the
matrices via these relationships rather than forming
directly can
have huge advantageous because it avoids having to create the (typically
dense) matrix
directly.
A GWoodburyMatrix
object for a non-symmetric
matrix, SWoodburyMatrix
for a symmetric matrix.
Where applicable, it's worth using the symmetric form of the matrix. This takes advantage of the symmetry where possible to speed up operations, takes less memory, and sometimes has numerical benefits. This function will create the symmetric form in the following circumstances:
symmetry = TRUE
; or
the argument X
is provided; or
A
and B
are symmetric (according to
isSymmetric
) and the arguments U
and V
are
NOT provided.
A more general form allows for multiple B matrices:
and analogously for the symmetric form. You can use this form by providing
a list of matrices as the B
(or U
, V
or X
)
arguments. Internally, this is implemented by converting to the standard form
by letting B = bdiag(...the B matrices...)
,
U = cbind(..the U matrices...)
, and so on.
The B
, U
, V
and X
values are recycled to the
length of the longest list, so you can, for instance, provide multiple B
matrices but only one U matrix (and vice-versa).
More information on the underlying linear algebra can be found in Harville, D. A. (1997) <doi:10.1007/b98818>.
WoodburyMatrix, solve, instantiate
library(Matrix) # Example solving a linear system with general matrices A <- Diagonal(100) B <- rsparsematrix(100, 100, 0.5) W <- WoodburyMatrix(A, B) str(solve(W, rnorm(100))) # Calculating the determinant of a symmetric system A <- Diagonal(100) B <- rsparsematrix(100, 100, 0.5, symmetric = TRUE) W <- WoodburyMatrix(A, B, symmetric = TRUE) print(determinant(W)) # Having a lower rank B matrix and an X matrix A <- Diagonal(100) B <- rsparsematrix(10, 10, 1, symmetric = TRUE) X <- rsparsematrix(100, 10, 1) W <- WoodburyMatrix(A, B, X = X) str(solve(W, rnorm(100))) # Multiple B matrices A <- Diagonal(100) B1 <- rsparsematrix(100, 100, 1, symmetric = TRUE) B2 <- rsparsematrix(100, 100, 1, symmetric = TRUE) W <- WoodburyMatrix(A, B = list(B1, B2)) str(solve(W, rnorm(100)))
library(Matrix) # Example solving a linear system with general matrices A <- Diagonal(100) B <- rsparsematrix(100, 100, 0.5) W <- WoodburyMatrix(A, B) str(solve(W, rnorm(100))) # Calculating the determinant of a symmetric system A <- Diagonal(100) B <- rsparsematrix(100, 100, 0.5, symmetric = TRUE) W <- WoodburyMatrix(A, B, symmetric = TRUE) print(determinant(W)) # Having a lower rank B matrix and an X matrix A <- Diagonal(100) B <- rsparsematrix(10, 10, 1, symmetric = TRUE) X <- rsparsematrix(100, 10, 1) W <- WoodburyMatrix(A, B, X = X) str(solve(W, rnorm(100))) # Multiple B matrices A <- Diagonal(100) B1 <- rsparsematrix(100, 100, 1, symmetric = TRUE) B2 <- rsparsematrix(100, 100, 1, symmetric = TRUE) W <- WoodburyMatrix(A, B = list(B1, B2)) str(solve(W, rnorm(100)))
The WoodburyMatrix
is a virtual class, contained by both
GWoodburyMatrix
(for general matrices) and SWoodburyMatrix
(for symmetric matrices). See WoodburyMatrix for construction
of these classes. The methods available for these classes are described
below; see also the solve methods. This class is itself a subclass of
Matrix
, so basic matrix methods like nrow
,
ncol
, dim
and so on also work.
## S4 method for signature 'GWoodburyMatrix' isSymmetric(object) ## S4 method for signature 'SWoodburyMatrix' isSymmetric(object) ## S4 method for signature 'GWoodburyMatrix,ANY' x %*% y ## S4 method for signature 'SWoodburyMatrix,ANY' x %*% y ## S4 method for signature 'GWoodburyMatrix' t(x) ## S4 method for signature 'SWoodburyMatrix' t(x)
## S4 method for signature 'GWoodburyMatrix' isSymmetric(object) ## S4 method for signature 'SWoodburyMatrix' isSymmetric(object) ## S4 method for signature 'GWoodburyMatrix,ANY' x %*% y ## S4 method for signature 'SWoodburyMatrix,ANY' x %*% y ## S4 method for signature 'GWoodburyMatrix' t(x) ## S4 method for signature 'SWoodburyMatrix' t(x)
object |
|
x |
|
y |
Matrix or vector |
GWoodburyMatrix-class
: Sub-class representing a generic matrix.
SWoodburyMatrix-class
: Sub-class representing a symmetric matrix.
Also subclasses symmetricMatrix.
isSymmetric,GWoodburyMatrix-method
: Check for symmetry of matrix; always returns
FALSE
.
isSymmetric,SWoodburyMatrix-method
: Check for symmetry of matrix; always returns
TRUE
.
%*%,GWoodburyMatrix,ANY-method
: Matrix multiplication (generally fast and
%*%,SWoodburyMatrix,ANY-method
: Matrix multiplication (generally fast and
t,GWoodburyMatrix-method
: Return the transpose of the matrix as
another GWoodburyMatrix.
t,SWoodburyMatrix-method
: Does nothing, just returns x
.
A
n x n subclass of Matrix
(GWoodburyMatrix
) or symmetricMatrix
(SWoodburyMatrix
).
B
p x p subclass of Matrix
(GWoodburyMatrix
) or symmetricMatrix
(SWoodburyMatrix
).
U
n x p subclass of Matrix
(only for
V
p x m subclass of Matrix
(only for
X
n x p subclass of Matrix
(only for
O
p x p subclass of Matrix
WoodburyMatrix for object construction, Matrix (the parent of this class).