Cholesky of a structured matrix

[markusm73]

I'd like to compute the Cholesky factor of the following obviously positive definite matrix of dimensions MxM:

(I + R*R^t)

where R is an upper trapezoidal matrix with dimensions MxN, and M << N. It seems this could be done more stably and efficiently without actually having to calculate R*R^t (unstable) followed by adding the diagonal of ones and then using dpotrf.

Is there anything in LAPACK that could help with this? Thanks!

[Julien Langou]

Hello,

One thing comes to my mind ... You can use updating / downdating algorithms for Cholesky factorization. I am not an expert so I do not know about the stability. This kind
of functionality was in LINPACK (ZCHUD and ZCHDD) but is not in LAPACK. An example of use is the "cholupdate" routine from Matlab. So you start from the Cholesky factorization of the identity
matrix (which is pretty easy to have ...) and then you update it.

Julien.

[sven]

Just to expand a little on what Julien suggested, if we put Rhat as

Rhat = ( I R ) then Rhat*Rhat^T = I + R*R^T,

so what you need is a QR factorizaton of Rhat^T,

Rhat^T = Q*U, U^T = ( Rtilde^T 0 ) and Rhat* Rhat^T = Rtilde^T*Rtilde,

so that Rtilde is the Cholesky factor of (I + R*R^T). This avoids forming R*R^T and is numerically stable. If your matrices are not too large and you don't mind being somewhat inefficient, then you could simply use the LAPACK QR factorization routine DGEQRF. Daniel Kressner has some nice QR updating routines at

http://www.math.ethz.ch/~kressner/qrupdate.php

Hope that helps,

Sven Hammarling.

[markusm73]

Thanks to everybody for the suggestions. I had thought of attaching the identity matrix to R and use QR factorization, but R already was the result of a QR factorization so I thought there must be a better way of solving my problem. In fact, I later found a better solution to the overall problem (inversion of a huge, sparse matrix) which, after exploiting the Woodbury formula for matrix inversion, allows me to also avoid inner products and use only one QR factorization to achieve the same result.