help smwoodbury

Syntax

real matrix smwoodbury_lu(real matrix Ainv, real matrix U, real matrix Cinv, real matrix V)

real matrix smwoodbury_lu_solve(real matrix Ainv, real matrix U, real matrix Cinv, real matrix V, real matrix AinvB)

real matrix smwoodbury_qr(real matrix Ainv, real matrix U, real matrix Cinv, real matrix V)

real matrix smwoodbury_qr_solve(real matrix Ainv, real matrix U, real matrix Cinv, real matrix V, real matrix AinvB)

real matrix smwoodbury_sym1_lu(real matrix Ainv, real matrix U)

real matrix smwoodbury_sym1_lu_solve(real matrix Ainv, real matrix U, real matrix AinvB)

real matrix smwoodbury_sym1_qr(real matrix Ainv, real matrix U)

real matrix smwoodbury_sym1_qr_solve(real matrix Ainv, real matrix U, real matrix AinvB)

real matrix smwoodbury_sym1_posdef(real matrix Ainv, real matrix U)

real matrix smwoodbury_sym1_posdef_solve(real matrix Ainv, real matrix U, real matrix AinvB)

real matrix smwoodbury_sym2_lu(real scalar s, real matrix U, real matrix Cinv)

real matrix smwoodbury_sym2_lu_solve(real scalar s, real matrix U, real matrix Cinv, real matrix sB)

real matrix smwoodbury_sym2_qr(real scalar s, real matrix U, real matrix Cinv)

real matrix smwoodbury_sym2_qr_solve(real scalar s, real matrix U, real matrix Cinv, real matrix sB)

real matrix smwoodbury_sym2_posdef(real scalar s, real matrix U, real matrix Cinv)

real matrix smwoodbury_sym2_posdef_solve(real scalar s, real matrix U, real matrix Cinv, real matrix sB)

Description

This package provides a set of Mata functions that implement the Sherman-Morrison-Woodbury formula for a rank-k update to a matrix inverse and that use the formula to efficiently solve matrix equations. The formula is:

inv(A+U*C*V) = inv(A) - inv(A)*U*inv(inv(C)+V*inv(A)*U)*V*inv(A)

where A is NxN, U is NxK, C is KxK, V is KxN.
Use of this formula can save computation time when N is much larger than K and either inv(A) is easy to calculate or inv(A+U*C*V) must be calculated many times for different values of U, C and V.
Such situations arise, for example, in Markov chain Monte Carlo estimation of mixed models.

Notice that the functions require as input Ainv and Cinv, the inverses of A and C, rather than A and C themselves.

smwoodbury_lu() and smwoodbury_qr() implement the full formula.
The former uses lusolve() and thus assumes Cinv+V*Ainv*U has full rank.
The latter uses qrsolve() and will work even when Cinv+V*Ainv*U does not have full rank (or is poorly conditioned), but is typically slower.

smwoodbury_sym1_lu() and smwoodbury_sym1_qr() implement the special case where C=I(K), V=U', and A is symmetric.
The former uses lusolve() and thus assumes I(K)+U'*Ainv*U has full rank. This assumption generally holds but can fail for unusual values of A and U.
The latter uses qrsolve() and will work even when I(K)+U'*Ainv*U does not have full rank (or is poorly conditioned), but is typically slower.
To maximize speed, smwoodbury_sym1_lu() and smwoodbury_sym1_qr() do not verify that Ainv is symmetric.
You may get unexpected results if Ainv is not symmetric.

smwoodbury_sym1_posdef() implements the special case where C=I(K), V=U', A is symmetric, and I(K)+U'*Ainv*U is positive definite.
To maximize speed, the function does not verify that Ainv is symmetric or that I(K)+U'*Ainv*U is positive definite. You may get unexpected results if these conditions fail.

smwoodbury_sym2_lu() and smwoodbury_sym2_qr() implement the special case where V=U' and A=I(N)/s for a nonzero scalar s.
The former uses lusolve() and thus assumes Cinv+s*U'*U has full rank. This assumption generally holds but can fail for unusual values of s and U.
The latter uses qrsolve() and will work even when Cinv+s*U'*U does not have full rank (or is poorly conditioned), but is typically slower.

smwoodbury_sym2_posdef() implements the special case where V=U', A=I(N)/s for a nonzero scalar s, Cinv is symmetric, and Cinv+s*U'*U is positive definite.
To maximize speed, the function does not verify that Cinv is symmetric or that Cinv+s*U'*U is positive definite. You may get unexpected results if these conditions fail.

Functions smwoodbury_*_solve() are equivalent to the functions without solve in the name but are designed to efficiently solve the equation (A+U*C*V)X=B for X. These functions require as input either AinvB=inv(A)*B or sB=(s*I(N))*B.

Author

Sam Schulhofer-Wohl

Federal Reserve Bank of Minneapolis

90 Hennepin Ave.

Minneapolis MN 55480-0291

wohls@minneapolisfed.org

The views expressed herein are those of the author and not necessarily those of the Federal Reserve Bank of Minneapolis or the Federal Reserve System.