Posts tagged ‘Matrix’

December 6, 2012

Rank 1 update to Cholesky factorization in R

Rank 1 update can be achieved in Matlab with the built-in function cholupdate(). More details about the function can be found here:

http://www.mathworks.com.au/help/matlab/ref/cholupdate.html

In R, the recommended Matrix package has added a function updown() to handle rank 1 update. However, the documentation of that function is not easy to understand, and I find the example given is misleading. That prompted me to replicate the Matlab example with the updown() function. The resulting code snippet is shown below.

library(matrixcalc)
A <- symmetric.pascal.matrix(4)
R <- chol(A)

x <- c(0, 0, 0, 1)

A + x%*%t(x)

library(Matrix)
R <- Cholesky(Matrix(A, sparse=T))
x <- Matrix(x)
R1.sparse <- updown('+', x, R)
R1 <- as(R1.sparse, 'Matrix')


x <- Matrix(c(0, 0, 0, 1/sqrt(2)))
R1.sparse <- updown('-', x, R)
R1 <- as(R1.sparse, 'Matrix')


Note that matrices have to be in the sparse representation to be used in the updown() function.

Advertisements