/* rank-k update of a matrix inverse using the sherman-morrison-woodbury formula */ /* sam schulhofer-wohl, federal reserve bank of minneapolis, sschulh1.work@gmail.com, 12/19/2011 */ /* formula: inv(A+U*C*V) = inv(A) - inv(A)*U*inv(inv(C)+V*inv(A)*U)*V*inv(A) assumptions: Ainv is NxN U is NxK C is KxK and full rank V is KxN */ version 12.0 mata: mata set matastrict on mata set matafavor speed /* general case using lusolve */ real matrix smwoodbury_lu(real matrix Ainv, real matrix U, real matrix Cinv, real matrix V) { real matrix work work=V*Ainv return(Ainv-Ainv*U*lusolve(Cinv+work*U,work)) } real matrix smwoodbury_lu_solve(real matrix Ainv, real matrix U, real matrix Cinv, real matrix V, real matrix AinvB) { real matrix work work=Ainv*U return(AinvB-work*lusolve(Cinv+V*work,V*AinvB)) } /* general case using qrsolve */ real matrix smwoodbury_qr(real matrix Ainv, real matrix U, real matrix Cinv, real matrix V) { real matrix work work=V*Ainv return(Ainv-Ainv*U*qrsolve(Cinv+work*U,work)) } real matrix smwoodbury_qr_solve(real matrix Ainv, real matrix U, real matrix Cinv, real matrix V, real matrix AinvB) { real matrix work work=Ainv*U return(AinvB-work*qrsolve(Cinv+V*work,V*AinvB)) } /* special case C=I(K), V=U', Ainv symmetric; using lusolve */ real matrix smwoodbury_sym1_lu(real matrix Ainv, real matrix U) { real matrix work work=Ainv*U return(Ainv-work*lusolve(I(cols(U))+U'*work,work')) } real matrix smwoodbury_sym1_lu_solve(real matrix Ainv, real matrix U, real matrix AinvB) { real matrix work work=Ainv*U return(AinvB-work*lusolve(I(cols(U))+U'*work,U'*AinvB)) } /* special case C=I(K), V=U', Ainv symmetric; using qrsolve */ real matrix smwoodbury_sym1_qr(real matrix Ainv, real matrix U) { real matrix work work=Ainv*U return(Ainv-work*qrsolve(I(cols(U))+U'*work,work')) } real matrix smwoodbury_sym1_qr_solve(real matrix Ainv, real matrix U, real matrix AinvB) { real matrix work work=Ainv*U return(AinvB-work*qrsolve(I(cols(U))+U'*work,U'*AinvB)) } /* special case C=I(K), V=U', Ainv symmetric and positive definite; using cholsolve */ real matrix smwoodbury_sym1_posdef(real matrix Ainv, real matrix U) { real matrix work work=Ainv*U return(Ainv-work*cholsolve(I(cols(U))+U'*work,work')) } real matrix smwoodbury_sym1_posdef_solve(real matrix Ainv, real matrix U, real matrix AinvB) { real matrix work work=Ainv*U return(AinvB-work*cholsolve(I(cols(U))+U'*work,U'*AinvB)) } /* special case V=U', A=I(N)/s; using lusolve */ real matrix smwoodbury_sym2_lu(real scalar s, real matrix U, real matrix Cinv) { real matrix work work=s:*U return(s*I(rows(U))-work*lusolve(Cinv+U'*work,work')) } real matrix smwoodbury_sym2_lu_solve(real scalar s, real matrix U, real matrix Cinv, real matrix sB) { real matrix work work=s:*U return(sB-work*lusolve(Cinv+U'*work,U'*sB)) } /* special case V=U', A=I(N)/s; using qrsolve */ real matrix smwoodbury_sym2_qr(real scalar s, real matrix U, real matrix Cinv) { real matrix work work=s:*U return(s*I(rows(U))-work*qrsolve(Cinv+U'*work,work')) } real matrix smwoodbury_sym2_qr_solve(real scalar s, real matrix U, real matrix Cinv, real matrix sB) { real matrix work work=s:*U return(sB-work*qrsolve(Cinv+U'*work,U'*sB)) } /* special case C symmetric and positive definite, V=U', A=I(N)/s; using cholsolve */ real matrix smwoodbury_sym2_posdef(real scalar s, real matrix U, real matrix Cinv) { real matrix work work=s:*U return(s*I(rows(U))-work*cholsolve(Cinv+U'*work,work')) } real matrix smwoodbury_sym2_posdef_solve(real scalar s, real matrix U, real matrix Cinv, real matrix sB) { real matrix work work=s:*U return(sB-work*cholsolve(Cinv+U'*work,U'*sB)) } mata mlib create lsmwoodbury, replace mata mlib add lsmwoodbury smwoodbury_*() mata mlib index end