Is there a way to simplify block Cholesky decomposition if you already have decomposed the submatrices along the leading diagonal?

Question

Let's say we have a block matrix $ M =\left( \begin{array}{ccc} A & B\\ B^{*} & C \end{array} \right)$ where $M$ is positive definite. ($A$ and $C$ are also positive definite.)

There is a formula for carrying out block Cholesky decomposition. See Wikipedia: Block LU decomposition. Summarising we have the following result.

The matrix $M = LU$ can be decomposed in an algebraic manner into

$$L = \begin{pmatrix} A^{\frac{1}{2}} & 0 \\ B^{*} A^{-\frac{*}{2}} & Q^{\frac{1}{2}} \end{pmatrix}$$

where $\begin{matrix} Q = C - B^{*} A^{-1} B \end{matrix}$ ($*$ indicates transpose in this case)

Now lets say we have already carried out the Cholesky decomposition for A, and C. So we have already calculated $A^{1/2}$, and $C^{1/2}$ (It is therefore straightforward to calculate the inverses $A^{-1/2}$, and $C^{-1/2}$ using forward substitution).

Rewriting $Q$ in terms of these quantities we now have

$$Q = Q^{1/2}Q^{*/2} = C^{1/2} C^{*/2} - (B^{*} A^{-*/2})(A^{-1/2} B)$ = $(C^{1/2} + B^{*}A^{-*/2})(C^{1/2} - B^{*}A^{-*/2})^{*}.$$

My question is this: Given this set up is it possible to algebraically calculate $Q^{1/2}$ without having to apply Cholesky decomposition to $Q$. Or, in other words, can I use $C^{1/2}$ to help me in the calculation of $Q^{1/2}$.

Answer 1

If A,C are fixed, and B is variable but nice (low-rank), then you want what is called "Cholesky update". If A,B,C are fixed, then probably you should not be picky about how the blocking is done, and you want to use a standard "block Cholesky". I have not found a clear answer for A,C fixed, B variable and not nice (I can ask around, so let me know if that really is your case).

---

Cholesky update

Rank one updates, chol(A) to chol(A+xx*), are easy and safe. Rank one "downdates", chol(A) to chol(A-xx*), are easy but require a little care: stable algorithms are given in Stewart's Matrix Algorithms Vol 1, Algorithm 4.3.8, p. 347. Chapter 12.5 of Golub–Van Loan has some similar stuff, and Cholesky down-dating in 12.5.4. This function has been widely implemented, and the cholupdate command in matlab dates back to 1979 code from LINPACK.

[0,B*;B,0] is a sum of rank one matrices, and so by updating and downdating those rank one guys, you could probably get what you want, and it might even be faster than chol(Q). However, it can be a lot better to update more ranks at a time.

Apparently this is a common request in machine learning, and M. Seeger wrote a technical report on this problem of low rank updates to a Cholesky factorization, and mentions several common pitfalls, especially as regards to actually doing it with existing software.

A more scholarly (and older) treatment is in section 3 of this article version of Ch. 12.5 of GvL:

Gill, P. E.; Golub, G. H.; Murray, W.; Saunders, M. A. "Methods for modifying matrix factorizations." Math. Comp. 28 (1974), 505–535. MR343558 DOI:10.2307/2005923

Davis and Hager in MR1824053 note that algorithm C1 can be used for a reasonably efficient, multiple rank, single pass, update of a dense matrix (and go on to describe sparse techniques).

Note that these mostly do not take advantage of the block structure of [0,B*;B,0], so you might find something better that is more specialized.

---

Block Cholesky

Blocking the Cholesky decomposition is often done for an arbitrary (symmetric positive definite) matrix. I didn't immediately find a textbook treatment, but the description of the algorithm used in PLAPACK is simple and standard.

In their algorithm they do not use the factorization of C, just of A. That allows them to reduce the problem of chol([A,B*;B,C]) to just chol(A) and chol(Q). The point of the algorithm is that you do not choose A and C to have the same size. You choose A to fit nicely in cache, and do your work at a higher BLAS level. In other words, A is a real block, and C is just leftover garbage you'll need to sweep up next.

In particular, C is discarded and replaced by Q during the algorithm, but chol(Q) is also computed by decomposing Q itself into a block matrix. This means that the algorithm is discarding any information you had about C, so if C is fixed while B varies, this would be quite wasteful.

Answer 2

Thanks to all above who replied to my question. I think I've come to an answer although it is not exactly as I'd hoped.

Removing the machine learning context, my question boiled down to whether knowing $C^{1/2}$ would help in the calculation of $Q^{-1/2}$. I'll go into more detail below but to cut the chase, the answer is yes, but only with respect to stability and not computation (can't prove this to be the case currently, but fairly certain).

For why the answer is yes w.r.t. to stability we look at the definition of $Q$ from the original question (I've updated the original as well):

$Q = C - B^{*} A^{-1} B = (C^{1/2} + B^{*}A^{-*/2})(C^{1/2} - B^{*}A^{-*/2})^{*}$

By knowing $C^{1/2}$ before hand, we can calculate $Q$ without having to invert $A$ directly. Direct inversion is not numerically stable. **It has been pointed out, see Jack's comment below, that both formulations of Q do not require A to be inverted directly. However the new form might be more stable if Q is of approximately low rank.

Sadly, although I have done a fair amount of research on the subject, it does not appear that $C^{1/2}$ helps wrt computation in the exact calculation of $Q^{-1/2}$. The best approach appears to be to calculate $Q$ using $C^{1/2}$ as above and then use cholesky to decompose $Q$ to $Q^{1/2}$ and then forward substitution to calculate $Q^{-1/2}$.

Further Research

One area I did not look into in much detail was whether it was possible to use $C^{1/2}$ to approximate $Q^{-1/2}$. Something along the lines of an iterative method using $C^{1/2}$ as a starting point. I do not know of any such iterative approximation process, but I'll keep searching. I may even start a new question with that as the focus.

Again, thanks to all who helped me getting to this answer. I'll update you all if I have any major breakthroughs.