Given a Cholesky factor,
R, of a matrix,
choldrop finds the Cholesky factor of
k is an integer.
cholup finds the factor of A+uu' (update) or A-uu' (downdate).
Cholesky factor of a matrix,
row and column of
vector defining rank one update.
R is upper triangular then
t(R[,-k])%*%R[,-k] == A[-k,-k], but
R[,-k] has elements on the first sub-diagonal, from its kth column onwards. To get from this to a triangular Cholesky factor of
A[-k,-k] we can apply a sequence of Givens rotations from the left to eliminate the sub-diagonal elements. The routine does this. If
R is a lower triangular factor then Givens rotations from the right are needed to remove the extra elements. If
n is the dimension of
R then the update has O(n^2) computational cost.
cholup (which assumes
R is upper triangular) updates based on the observation that R'R + uu' = [u,R'][u,R']' = [u,R']Q'Q[u,R']', and therefore we can construct Q so that Q[u,R']'=[0,R1']', where R1 is the modified factor. Q is constructed from a sequence of Givens rotations in order to zero the elements of u. Downdating is similar except that hyperbolic rotations have to be used in place of Givens rotations — see Golub and van Loan (2013, section 6.5.4) for details. Downdating only works if A-uu' is positive definite. Again the computational cost is O(n^2).
Note that the updates are vector oriented, and are hence not susceptible to speed up by use of an optimized BLAS. The updates are set up to be relatively Cache friendly, in that in the upper triangular case successive Givens rotations are stored for sequential application column-wise, rather than being applied row-wise as soon as they are computed. Even so, the upper triangular update is slightly slower than the lower triangular update.
Simon N. Wood firstname.lastname@example.org
Golub GH and CF Van Loan (2013) Matrix Computations (4th edition) Johns Hopkins
require(mgcv) set.seed(0) n <- 6 A <- crossprod(matrix(runif(n*n),n,n)) R0 <- chol(A) k <- 3 Rd <- choldrop(R0,k) range(Rd-chol(A[-k,-k])) Rd;chol(A[-k,-k]) ## same but using lower triangular factor A = LL' L <- t(R0) Ld <- choldrop(L,k) range(Ld-t(chol(A[-k,-k]))) Ld;t(chol(A[-k,-k])) ## Rank one update example u <- runif(n) R <- cholup(R0,u,TRUE) Ru <- chol(A+u %*% t(u)) ## direct for comparison R;Ru range(R-Ru) ## Downdate - just going back from R to R0 Rd <- cholup(R,u,FALSE) R0;Rd range(R-Ru)
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.