Description Usage Arguments Details Value Symmetric form Multiple B matrices References See Also Examples
View source: R/WoodburyMatrix.R
Creates an implicitly defined matrix representing the equation
A^{1} + U B^{1} V,
where A, U, B and V are n x n, n x p, p x p and p x n matrices, respectively. A symmetric special case is also possible with
A^{1} + X B^{1} X',
where X is n x p and A and B are additionally symmetric. The available methods are described in WoodburyMatrixclass and in solve. Multiple B / U / V / X matrices are also supported; see below
1  WoodburyMatrix(A, B, U, V, X, O, symmetric)

A 
Matrix A in the definition above. 
B 
Matrix B in the definition above, or list of matrices. 
U 
Matrix U in the definition above, or list of matrices. Defaults to a diagonal matrix/matrices. 
V 
Matrix V in the definition above, or list of matrices. Defaults to a diagonal matrix/matrices. 
X 
Matrix X in the definition above, or list of matrices. Defaults to a diagonal matrix/matrices. 
O 
Optional, precomputed value of O, as defined above. THIS IS NOT CHECKED FOR CORRECTNESS, and this argument is only provided for advanced users who have precomputed the matrix for other purposes. 
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
A  A U O^{1} V A
where O = B + VAU, and its determinant by
det(O) det(A)^{1} det(B)^{1}.
These relationships are often called the Woodbury matrix identity and the matrix determinant lemma, respectively. If A and B are sparse or otherwise easy to deal with, and/or when p < n, manipulating the matrices via these relationships rather than forming W directly can have huge advantageous because it avoids having to create the (typically dense) matrix
A^{1} + U B^{1} V
directly.
A GWoodburyMatrix
object for a nonsymmetric
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:
A^{1} + ∑_{i = 1}^n U_i B_i^{1} V_i,
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 viceversa).
More information on the underlying linear algebra can be found in Harville, D. A. (1997) <doi:10.1007/b98818>.
WoodburyMatrix, solve, instantiate
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  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)))

Formal class 'dgeMatrix' [package "Matrix"] with 4 slots
..@ x : num [1:100] 0.905 1.224 2.072 1.015 1.75 ...
..@ Dim : int [1:2] 100 1
..@ Dimnames:List of 2
.. ..$ : NULL
.. ..$ : NULL
..@ factors : list()
$modulus
[1] 2.192931
attr(,"logarithm")
[1] TRUE
$sign
[1] 1
attr(,"class")
[1] "det"
Formal class 'dgeMatrix' [package "Matrix"] with 4 slots
..@ x : num [1:100] 0.138 1.853 0.281 0.238 0.229 ...
..@ Dim : int [1:2] 100 1
..@ Dimnames:List of 2
.. ..$ : NULL
.. ..$ : NULL
..@ factors : list()
Formal class 'dgeMatrix' [package "Matrix"] with 4 slots
..@ x : num [1:100] 11.35 6.63 4.18 16.9 11.17 ...
..@ Dim : int [1:2] 100 1
..@ Dimnames:List of 2
.. ..$ : NULL
.. ..$ : NULL
..@ factors : list()
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.