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.