Cholesky of a structured matrix

Open discussion regarding features, bugs, issues, vendors, etc.

Cholesky of a structured matrix

Postby markusm73 » Wed Sep 19, 2012 12:45 pm

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!
Posts: 18
Joined: Tue Jun 14, 2005 8:10 pm
Location: New York

Re: Cholesky of a structured matrix

Postby Julien Langou » Wed Sep 19, 2012 4:38 pm


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 Langou
Posts: 835
Joined: Thu Dec 09, 2004 12:32 pm
Location: Denver, CO, USA

Re: Cholesky of a structured matrix

Postby sven » Thu Sep 20, 2012 5:22 am

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

Hope that helps,

Sven Hammarling.
Posts: 146
Joined: Wed Dec 22, 2004 4:28 am

Re: Cholesky of a structured matrix

Postby markusm73 » Fri Sep 21, 2012 3:37 pm

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.
Posts: 18
Joined: Tue Jun 14, 2005 8:10 pm
Location: New York

Return to User Discussion

Who is online

Users browsing this forum: No registered users and 7 guests