7
$\begingroup$

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)

$\endgroup$
  • 1
    $\begingroup$ This is coincidentially solved in the SLAM++ library. See sf.net/p/slam-plus-plus and the associated publications. $\endgroup$ – the swine Nov 3 '16 at 8:45
6
$\begingroup$

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$.

$\endgroup$
  • 1
    $\begingroup$ Cough. The Cholesky factor of Cxx is the upper-left triangle of the Cholesky factor of C. $\endgroup$ – Jack Poulson Jan 19 '13 at 19:13
  • $\begingroup$ Yup, thanks for reminding me what the original question was asking for. $\endgroup$ – Jed Brown Jan 19 '13 at 20:51
  • $\begingroup$ 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 $\endgroup$ – robince Jan 21 '13 at 11:23
  • 1
    $\begingroup$ hi, just bumped on this now. For my application $A$ is huge, $B$ is just a vector which keeps changing, and $D$ is zero. Thus, this is really valuable, as updation seems really straightforward. However, $A$ in my application is positive semi-definite, so how can I extend this to that? $\endgroup$ – dineshdileep May 14 '16 at 2:38
  • $\begingroup$ @dineshdileep You can use the pseudoinverse instead of the inverse for $A$, which will have the same result. $\endgroup$ – Guillermo Angeris Jun 3 '18 at 1:53

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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