Is there a way to simplify block Cholesky decomposition If you already have decomposed the sub matrices along the leading diagonal?

Lets 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 pos def)

There is a formula for carrying out block Cholesky decomposition. See http://en.wikipedia.org/wiki/Block_LU_decomposition. Summarising we have the following result.

The matrix $$\begin{matrix}M = LU\end{matrix}$$ can be decomposed in an algebraic manner into $[itex]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 the 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 algebraicly 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}$.

Thanks in advance for any replies.

Matt.

-
Don't you mean "$B$ symmetric definite"? You have a lot of typos (including the title) and I encourage you to correct them. – Wadim Zudilin Jul 1 '10 at 10:49
Matt, reaching the end of your question I would say it's sounds like either homework or a very technical problem. It's not a real math problem! Please check FAQ. – Wadim Zudilin Jul 1 '10 at 10:52
Hi Wadim. Thanks for the reply. The problem is indeed technical in its origin , but I'd hoped (perhaps naively) that the problem would also be of interest to other mathematicians. The problem is related to the training a machine learning algorithm. As part this training a positive-definite matrix (covariance matrix) is decomposed using Cholesky decomposition. I'm trying to work out if I can parallelise the training. In this case A and C would correspond to two different training sets. If the block decomposition can be simplified as described it would make this possible. – Matthew Gretton Jul 1 '10 at 11:57
If you think it really doesn't qualify and would not be of interest to mathematicians on the site then I think it's fair to close it. Thanks for the reply all the same. Matt. – Matthew Gretton Jul 1 '10 at 11:58
It seems a fairly natural question to me when you're doing numerical analysis. In my view, most questions here are more specialized than this one. – Jitse Niesen Jul 1 '10 at 14:32

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.

-
Hey. Thanks for the reply. Sadly it's looking like I can't use C^{-1/2} to help be out in the decomposition of Q. Perhaps some kind of approximation to Q^{1/2} using C^{-1/2} as a starting point. I'll have a think... Should have time over the weekend to look over the links you supplied in more detail. I'll reply in more detail then. Thanks again for the detailed reply. – Matthew Gretton Jul 2 '10 at 13:01
Hi. I've got a bit of time to respond to your answer properly. With respect to your first paragraph, as a result of the way the training works given two sets of training samples you want to combine, A,B, and C are fixed. So as you point out we are updating the decomposition of A, and again as you say, the best way to do this is to do the update in one go (calc Q^{-1/2} using cholesky decomposition). So the one question remaining is whether there is any way I can use C^{-1/2} in any way to speed up the calculation of Q^{-1/2}. – Matthew Gretton Jul 5 '10 at 22:15
I'm still hopeful. I have a gut feeling that the computation of C^{-1/2} must in some way give some relevant information and should therefore be able to speed up the calculation. Do you know any Matrix Algebra wizs who would be able to verify one way or the other? Thanks again. – Matthew Gretton Jul 5 '10 at 22:19
Summer is a slow time; so far none around. However, I am doubtful it would help much in your case to compute C^{1/2} (do you need it for anything else?). Imagine A=1, C=0. Then Q = -B*B is an arbitrary (negative) definite symmetric matrix, completely unrelated to A and C. How can knowing A^{1/2}=1 and C^{1/2}=0 possibly help to find the (imaginary) Cholesky decomposition of the arbitrary matrix Q? – Jack Schmidt Jul 7 '10 at 19:11
Thanks! That exactly answered my earlier question. I would phrase your problem as A,C fixed, B changes, just because A,C came (possibly hours) before B was around, and so pre-computation is completely sensible for A,C and impossible for B. I personally still think C cannot help (well, your stability argument has halfway changed my mind), but I'll definitely ask around once people get back. Our numerical linear algebraists are all out of the country right now, but will be back soon enough. – Jack Schmidt Jul 8 '10 at 16:54

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

-
Are you sure it helps with stability (in general)? You are still using A^-1/2 in your expression for Q. To compute BA^-1B, surely you would use the Cholesky decomposition of A to compute A^-1B (back solve), and so that already accounts for all the instability of A^-1: it occurs in both methods. Of course if Q is approximately low rank, then your factorization might help avoid some instability from the cancellation of C and BA^-1B, but I think that is separate. That said, your formula for Q cannot hurt, and has approximately the same operation count (n^2 more adds, but no big deal). – Jack Schmidt Jul 7 '10 at 19:16
Yes. That's certainly true. Direct inversion of A^-1 is avoided in both cases... Q could potentially be low rank so as you say can't hurt to use the reformulation of Q. – Matthew Gretton Jul 8 '10 at 0:31