help smwoodbury-------------------------------------------------------------------------------

Syntax

real matrixsmwoodbury_lu(real matrixAinv,real matrixU,real matrixCinv,real matrixV)

real matrixsmwoodbury_lu_solve(real matrixAinv,real matrixU,real matrixCinv,real matrixV,real matrixAinvB)

real matrixsmwoodbury_qr(real matrixAinv,real matrixU,real matrixCinv,real matrixV)

real matrixsmwoodbury_qr_solve(real matrixAinv,real matrixU,real matrixCinv,real matrixV,real matrixAinvB)

real matrixsmwoodbury_sym1_lu(real matrixAinv,real matrixU)

real matrixsmwoodbury_sym1_lu_solve(real matrixAinv,real matrixU,real matrixAinvB)

real matrixsmwoodbury_sym1_qr(real matrixAinv,real matrixU)

real matrixsmwoodbury_sym1_qr_solve(real matrixAinv,real matrixU,real matrixAinvB)

real matrixsmwoodbury_sym1_posdef(real matrixAinv,real matrixU)

real matrixsmwoodbury_sym1_posdef_solve(real matrixAinv,real matrixU,real matrixAinvB)

real matrixsmwoodbury_sym2_lu(real scalars,real matrixU,real matrixCinv)

real matrixsmwoodbury_sym2_lu_solve(real scalars,real matrixU,real matrixCinv,real matrixsB)

real matrixsmwoodbury_sym2_qr(real scalars,real matrixU,real matrixCinv)

real matrixsmwoodbury_sym2_qr_solve(real scalars,real matrixU,real matrixCinv,real matrixsB)

real matrixsmwoodbury_sym2_posdef(real scalars,real matrixU,real matrixCinv)

real matrixsmwoodbury_sym2_posdef_solve(real scalars,real matrixU,real matrixCinv,real matrixsB)

DescriptionThis 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()andsmwoodbury_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()andsmwoodbury_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()andsmwoodbury_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()andsmwoodbury_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 withoutsolvein 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.

AuthorSam 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.