11
$\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
  • 2
    $\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, 2016 at 8:45

1 Answer 1

13
$\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]: $$\begin{aligned} C & =\begin{pmatrix} A & B^{T}\\ B & D \end{pmatrix}\\ & =\begin{pmatrix} \mathbb{1} & 0\\ B\left( L_{A} L_{A}^{T}\right)^{-1} & \mathbb{1} \end{pmatrix}\begin{pmatrix} L_{A} L_{A}^{T} & 0\\ 0 & L_{S} L_{S}^{T} \end{pmatrix}\begin{pmatrix} \mathbb{1} & \left( L_{A} L_{A}^{T}\right)^{-1} B^{T}\\ 0 & \mathbb{1} \end{pmatrix}\\ & =\left(\begin{pmatrix} \mathbb{1} & 0\\ B\left( L_{A} L_{A}^{T}\right)^{-1} & \mathbb{1} \end{pmatrix}\begin{pmatrix} L_{A} & 0\\ 0 & L_{S} \end{pmatrix}\right)\left(\begin{pmatrix} L_{A}^{T} & 0\\ 0 & L_{S}^{T} \end{pmatrix}\begin{pmatrix} \mathbb{1} & \left( L_{A} L_{A}^{T}\right)^{-1} B^{T}\\ 0 & \mathbb{1} \end{pmatrix}\right)\\ & =\begin{pmatrix} L_{A} & 0\\ BL_{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} . \end{aligned} $$ 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$
5
  • 1
    $\begingroup$ Cough. The Cholesky factor of Cxx is the upper-left triangle of the Cholesky factor of C. $\endgroup$ Jan 19, 2013 at 19:13
  • $\begingroup$ Yup, thanks for reminding me what the original question was asking for. $\endgroup$
    – Jed Brown
    Jan 19, 2013 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, 2013 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$ May 14, 2016 at 2:38
  • $\begingroup$ @dineshdileep You can use the pseudoinverse instead of the inverse for $A$, which will have the same result. $\endgroup$ Jun 3, 2018 at 1:53

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

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