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.