| choldrop {mgcv} | R Documentation |
Given a Cholesky factor, R, of a matrix, A, choldrop finds the Cholesky factor of A[-k,-k],
where k is an integer. cholup finds the factor of A+uu' (update) or A-uu' (downdate).
choldrop(R,k) cholup(R,u,up)
R |
Cholesky factor of a matrix, |
k |
row and column of |
u |
vector defining rank one update. |
up |
if |
First consider choldrop. If 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 simon.wood@r-project.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)