3
$\begingroup$

I have been looking at the modified Cholesky decomposition suggested by the following paper: Schnabel and Eskow, A Revised Modified Cholesky Factorization Algorithm, SIAM J. Optim. 9, pp. 1135-1148 (14 pages).

The paper talks about an implementation of the Cholesky decomposition, modified by the fact that when the matrix is not positive definite, a diagonal perturbation is added before the decomposition takes place. The algorithm given in the paper (Algorithm 1) suggests that it finds factorization such that $LL^\top = A + E, E \ge 0$. But, when implemented as stated in the paper, what I get is a lower triangular factor matrix $L$, such that: $$P\cdot(L\cdot L^\top)\cdot P^\top = A + E$$ where, $P$ is a permutation matrix.

This $L$ matrix is not the same as when using the normal Cholesky factorization for a PD matrix, where $A = L_1L_1^\top$, say.

Now, I am using it in optimization context, specifically in trust region methods, where this decomposition is followed by inverting $L$ (to compute the trust region step), so it would be helpful to have factors in lower triangular form. Is there a way to get back $L_1$ (the original Cholesky factor) for a positive definite matrix from the $P$ and $L$ matrix obtained from the modified Cholesky factorization? I am a bit surprised that the algorithm misstates what it would produce, so maybe I am missing some step here.

Related threads I have found so far:

  1. Cholesky factorization
  2. Cholesky decomposition and permutation matrix

The second thread says that the relationship is not possible to find for any given set of matrix (not necessarily for a modified Cholesky decomposition as mentioned here). Not sure if the answer still holds in this case as well.

Example: Consider the following $4\times4$ matrix which is PD. $$A = \begin{bmatrix}6 & 3 & 4 & 8 \\ 3 & 6 & 5 & 1 \\ 4 & 5 & 10 & 7 \\ 8 & 1 & 7 & 25 \end{bmatrix}$$ The vanilla Cholesky factor $L_1$, such that $L_1L_1^\top=A$ is: $$L_1 = \begin{bmatrix}2.4495 & 0 & 0 & 0 \\ 1.2247 & 2.1213 & 0 & 0 \\ 1.6330 & 1.4142 & 2.3094 & 0 \\ 3.2660 & -1.4142 & 1.5877 & 3.1325 \end{bmatrix}$$ Now, if I perform the modified Cholesky, I get: $$L=\begin{bmatrix}5 & 0 & 0 & 0 \\ 1.4 & 2.83549 & 0 & 0 \\ 0.2 & 1.66462 & 1.78579 & 0 \\ 1.6 & 0.62070 & 0.92215 & 1.48471 \end{bmatrix}$$ and $$P=\begin{bmatrix}0 & 0 & 0 & 1 \\ 0 & 0 & 1 & 0 \\ 0 & 1 & 0 & 0 \\ 1 & 0 & 0 & 0 \end{bmatrix}$$ Such that, $P\cdot (L\cdot L^\top)\cdot P^\top = A$. Of course, since $A$ is PD, $E=0$.

  • 0
    If the input matrix is positive definite to begin with, the Cholesky triangle produced by the modified version is the same as the one obtained from vanilla Cholesky.2011-12-06
  • 0
    @J.M.: Thanks for replying. My implementation of the algorithm did not give the same lower triangular matrix as the vanilla Cholesky for pd matrices. I have also checked the implementation with a few other and the answers match. I think the reason is that, PD or not, the algorithm involves pivoting and hence the answer does not remain same. The only relationship that holds is $P*(LL^T)*P^T=A+E$.2011-12-06
  • 0
    I need to step out for now. I'll try to look at my implementation later and get back to you.2011-12-06
  • 0
    Okay, just so we have a nice baseline/diagnostic: could you try out your algorithm on a Hilbert matrix (say, $5\times 5$) and see if the algorithm does pivot? Also maybe include the Cholesky triangle you get for good measure.2011-12-06
  • 0
    @J. M.: I just added an example of a 4X4 matrix. Please let me know if it makes sense now. Thanks again.2011-12-06
  • 0
    If I want to "normalize" the cholesky-decomposition I can apply a (column-wise) rotation to another form, which represents a row-permutation. Here it seems, the rows are permuted to get the diagonal values sorted in descending order?2011-12-06
  • 0
    Hmm, I only saw the paper you are referring to now; it seems that this is a bit different from the perturbed Cholesky implementation I am accustomed to (what they refer to as the GMW version). I'll try a few experiments...2011-12-07

3 Answers 3