Take the 2-minute tour ×
Computational Science Stack Exchange is a question and answer site for scientists using computers to solve scientific problems. It's 100% free, no registration required.

I have a block matrix (either 2x2 blocks or 3x3 blocks) which is the covariance matrix for a joint space of two or three multivariate normal variables. ie

C = [Cxx Cxy; Cxy' Cyy];

I need to compute the cholesky factorisation of this matrix (C), as well as the diagonal blocks (Cxx, Cyy, the covariance matrices of the individual multivariate normals), and I would like to do this as fast as possible. At the moment I am doing three chol decompositions. I was wondering if it would be possible to obtain chol(Cxx) and chol(Cyy) from chol(C) (i.e. from extracting subblocks of the full decomposition) or if there would be any other trick to help do this faster.

(I have looked at QR factorisation instead of explicitly calculating the covariance but for my case it is many times slower)

share|improve this question
add comment

1 Answer

There is some redundant work, but not as much as you were hoping. Recall the $LDL^T$ factorization $$C = \begin{pmatrix} A & B^T \\ B & D \end{pmatrix} =\begin{pmatrix} 1 & 0 \\ B A^{-1} & 1 \end{pmatrix} \begin{pmatrix} A & 0 \\ 0 & S \end{pmatrix} \begin{pmatrix} 1 & A^{-1} B^T \\ 0 & 1 \end{pmatrix}$$ where $S = D - B A^{-1} B^T$ is the Schur complement. Clearly you need $A^{-1}$ and $S^{-1}$ to solve using this factorization, with $D^{-1}$ providing no value.

Substituting $A = L_A L_A^T$ and $S = L_S L_S^T$ into the factorization above, we obtain the Cholesky factorization [1] $$C = \begin{pmatrix} A & B^T \\ B & D \end{pmatrix} =\begin{pmatrix} L_A & 0 \\ B L_A^{-T} & L_S \end{pmatrix} \begin{pmatrix} L_A^T & L_A^{-1} B^T \\ 0 & L_S^T \end{pmatrix} = L_C L_C^T.$$ Thus you can simply extract the factors of $A$ from the factors of $C$, but will have to factor $D$ separately. If the block sizes are equal, $C$ takes eight times more work to factor, so this is little savings.

If you have a sequence of problems and can partition your system so that $A$ does not change between subsequent solves, its factorization could be reused, but $S$ will change and have to be refactored. Depending on the decomposition and conditioning, iterative methods can be useful to reuse a partial factorization.

[1] Note that this "block Cholesky" factorization is also valid with $L_A$ replaced by a $A^{1/2}$ or another non-triangular factorization of $A$.

share|improve this answer
1  
Cough. The Cholesky factor of Cxx is the upper-left triangle of the Cholesky factor of C. –  Jack Poulson Jan 19 '13 at 19:13
    
Yup, thanks for reminding me what the original question was asking for. –  Jed Brown Jan 19 '13 at 20:51
    
Thanks - its much clearer to me with the edit - so I can take chol(Cxx) from a subblock but not chol(Cyy). I don't really understand where the asymmtetry comes from... at least with simple numerical examples I also get Cyy from lower right block: pastebin.com/GM7Xr3Xj –  robince Jan 21 '13 at 11:23
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.