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
choldrop(R,k)
cholup(R,u,up)
Cholesky factor of a matrix, A
.
row and column of A
to drop.
vector defining rank one update.
if TRUE
compute update, otherwise downdate.
Simon N. Wood simon.wood@r-project.org
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
cholup
(which assumes R
is upper triangular) updates based on the observation that
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.
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(R0-Rd)
Run the code above in your browser using DataLab