# Modified Cholesky factorization and retrieving the usual LT matrix

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:

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

• 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. – J. M. isn't a mathematician Dec 6 '11 at 2:03
• @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$. – Samik R Dec 6 '11 at 3:16
• I need to step out for now. I'll try to look at my implementation later and get back to you. – J. M. isn't a mathematician Dec 6 '11 at 3:22
• 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. – J. M. isn't a mathematician Dec 6 '11 at 9:00
• @J. M.: I just added an example of a 4X4 matrix. Please let me know if it makes sense now. Thanks again. – Samik R Dec 6 '11 at 17:35

## 3 Answers

Here is an example of how the effect of columnwise rotation occurs. See the protocol using my MatMate-program. The "rot"-command performs a column-wise rotation to triangular shape, where the order of the invovled rows (the rotation-citeria) is given as list, and the triangular form is so just permuted. However, the permutation includes also a re-computation of the diagonal (and the other entries, too)

;===========================================
;MatMate-Listing vom:06.12.2011 19:40:03
;============================================
 A =  mk(4,4,6,3,4,8, 3,6,5,1, 4,5,10,7, 8,1,7, 25)
 L = cholesky(A)
A :
6.00        3.00        4.00        8.00
3.00        6.00        5.00        1.00
4.00        5.00       10.00        7.00
8.00        1.00        7.00       25.00

L :
2.45        0.00        0.00        0.00
1.22        2.12        0.00        0.00
1.63        1.41        2.31        0.00
3.27       -1.41        1.59        3.13

 M = rot(L,"drei",1´2´3´4)


Command  does not change anything because the order given is the default order. In

 M = rot(L,"drei",2´3´4´1)
M :
1.22        0.62        1.38       -1.48
2.45        0.00        0.00        0.00
2.04        2.42       -0.00        0.00
0.41        2.55        4.28       -0.00

 M = rot(L,"drei",3´4´1´2)
M :
1.26        1.16        1.75        0.00
1.58       -0.56        0.94        1.52
3.16       -0.00       -0.00       -0.00
2.21        4.48        0.00       -0.00


in command  we find the solution with the number 5 in the first column:

 M = rot(L,"drei",4´1´2´3)
M :
1.60        1.85       -0.00       -0.00
0.20        1.44        1.97        0.00
1.40        0.95        1.70       -2.06
5.00        0.00       -0.00       -0.00


After that I found the "modified cholesky" solution by the rotation in the following order:

 M = rot(L,"drei",4´3´2´1)
M :
1.60        0.62        0.92        1.48
0.20        1.66        1.79        0.00
1.40        2.84        0.00       -0.00
5.00        0.00        0.00        0.00


only that the rows are permuted.

So we can even systematically recover the original vanilla-cholesky matrix from that of the modified-cholesky-procedure: if we have r rows, then we need at most r(r-1)/2 rotations to triangular shape to recover the original vanilla-cholesky matrix.

[update 2] A remark about "rotation of columns". Rotation of columns can be thought as postmultiplication of the cholesky L by some rotation-matrix T. T itself is generated by successively rotating pairs of columns, so T can be seen as product of elementary rotation-matrices $\small t_{k,j}=\begin{array} {rrrr} \cos(\varphi_{k,j}) & \sin(\varphi_{k,j}) \\ -\sin(\varphi_{k,j}) & \cos(\varphi_{k,j}) \end{array}$ where that k,j'th rotation-parameters are inserted in the ID-matrix at the appropriate pair of rows k and j and the same columns. Then $\small T = t_{1,2} \cdot t_{1,3} \cdot t_{1,4} \cdot \ldots t_{2,3} \cdot t_{2,4} \cdot \ldots t_{3,4} \cdot \ldots \cdot t_{n-1,n}$ where n is the number of columns of the matrix which is to be rotated. In particular, T is not a permuationmatrix!
The different rotations mentioned in the example above occur, because for the triangular rotation we need to find a rotationmatrix $\small T_1$ which rotates the matrix L such that the entries in row 1 are collected in column 1, call this version of L $\small L_1$. Next we need to find the rotationmatrix $\small T_2$ which rotates then $\small L_1$ such that the entries in row 2 are collected in column 2 but where we do not touch the column 1. Save this in matrix $\small L_2$ and so on until $\small L_n = T_1 \cdot T_2 \cdot T_3 \cdot \ldots \cdot T_n$ is a triangular matrix. The different version in the example above is then simply, that the order of the $\small T_k$ in that product is changed (according to the list, given in the rot(...)-command in MatMate. (you may try this using MatMate yourself)
[update] The "modified cholesky" seems to proceed by extracting the row/column (thr "factor") which has the highest entry on the diagonal ("variance") first . Then proceeds with the reduced matrix after extraction of the next row/column. This explains the final permutation order between the rows 1 and 2, which have initially the same diagonal value 6, and the factor in row 1 has even a greater "individual variance" than the factor in row 2. See the following protocol, where I used a copy "chk" of the original matrix A and proceeded to extract always that factor with the highest variance in the (residual) covariance-matrix (MatMate-command "RemVar" ("RemoveVariance of one factor"). Negative signs at the zero is spurious numerical machine-epsilon:

 chk = A
chk :
6.00        3.00        4.00        8.00
3.00        6.00        5.00        1.00
4.00        5.00       10.00        7.00
8.00        1.00        7.00       25.00

 chk = remvar(chk,4)
chk :
3.44        2.68        1.76        0.00
2.68        5.96        4.72        0.00
1.76        4.72        8.04        0.00
0.00        0.00        0.00        0.00

 chk = remvar(chk,3)
chk :
3.05        1.65       -0.00        0.00
1.65        3.19       -0.00        0.00
-0.00       -0.00       -0.00        0.00
0.00        0.00        0.00        0.00

 chk = remvar(chk,2)
chk :
2.20        0.00        0.00        0.00
0.00        0.00        0.00        0.00
0.00        0.00       -0.00        0.00
0.00        0.00        0.00        0.00

• Thanks for your detailed comments. However, I am feeling slightly lost. The first part of your post seem to suggest that there is a way to recover the vanilla Cholesky factor given the modified Cholesky factor and the permutation matrix. But I did not get the procedure (or "steps") to do that. Regarding the second part of your post: you are correct, that is how the algorithm proceeds (by pivoting on the maximum diagonal), but I am not sure what your example demonstrates. Sorry, can you please help me out some more. – Samik R Dec 6 '11 at 23:24
• @SamikR: at the second part: I just did not really know the implementation of the modified cholesky-decomposition, so I tried just to re-engineer it. at the first part: you can bring any matrix to lower triangular form just by rotating between pairs of columns: first you focus on the first row and do pairwise rotations until all except the first entries are zero. Then you focus on the next row and do the same, only that you do not touch the first column. So you proceed until your matrix is triangular. If you do the same but with arbitrary change of focus to the rows in another order ... – Gottfried Helms Dec 6 '11 at 23:35
• you get a permutation of the triangular shape: you can do that the "top left" edge occurs in, say, row 3, and the (full) "last line" in, say, row 1 - just by these column-rotations. All possible (modified) cholesky-solutions of the same matrix A are one of that rotated scrumbled triangular rotations, so there are r(r-1)/2 possible permuted (modified) cholesky-decompositions; one of it is selected by your software. My goal with the first part of my post was a) to show that rotational effect and b) to re-engineer the modified version from the vanilla one to uncover the applied permutation. – Gottfried Helms Dec 6 '11 at 23:40
• So, you give me a matrix A, the modified cholesky-matrix and I find the permutation in it by doing the vanilla cholesky on my own and column-rotation until I get the same coefficients -only in row-permuted order. This exhibits then the hidden/unknown "modified-cholesky-permutation-matrix" – Gottfried Helms Dec 6 '11 at 23:47
• Thanks for the further explanations. Let me digest and get back - I think I now have enough information. Summary: I need to apply rotation in order to get back the vanilla Cholesky factor. The rotations that are needed are, I am guessing, to be somehow found from the permutation matrix, ostensibly from the order of pivoting that has been performed in the algorithm. Note that: the algorithm in the paper produces both the L (modified Cholesky factor matrix) and P (the permutation matrix). The P matrix is not unknown, as you mention in the last line. – Samik R Dec 7 '11 at 0:30

Perhaps off topic, but I hope it helps. Matrix decompositions often include permutations (also called pivoting) for stability and these are not mentioned. For instance, when computational people talk about an LU-decomposition of a matrix $A$ they usually mean a decomposition $A = PLU$ (or even $A = PLUQ$) where $P$ is a permutation matrix.

The reason for this is that permutation matrices often do not matter. In your case, you have Cholesky with symmetric pivoting: $A = PLL^TP^T$ (forgetting about $E$). If I understand you correctly, you want to use the decomposition to solve a system like $Ax=b$. The solution is $b = (PLL^TP^T)^{-1} x$, which you can compute by solving linear systems with the permutation matrix $P$ and the triangular matrix $L$ and their transpose. Solving a system with $P$ is simply a matter of finding the inverse permutation (the inverse of a permutation matrix is the matrix corresponding to the inverse of the underlying permutation), which does not cost a lot of effort.

• "Solving a system with $\mathbf P$ is simply a matter of finding the inverse permutation (the inverse of a permutation matrix is the matrix corresponding to the inverse of the underlying permutation), which does not cost a lot of effort." - indeed, one just swaps components of your vectors around for this. If the permutations are stored in an appropriate format, one can easily do either $\mathbf P$ or $\mathbf P^\top\mathbf x$. – J. M. isn't a mathematician Dec 7 '11 at 12:19

If I may suggest an alternative approach. There is an algorithm for decomposing a symmetric matrix (not necessarily positive definite) in LDL^T form, as efficient and stable as Cholesky. L is lower diagonal, and D a diagonal of entries 1 or -1. The number of -1 entries equals the number of negative eigenvalues. If the matrix is positive-definite, D is the unit matrix and L your usual Cholesky result.

Originally (as far as I know) published in a nuclear-physics paper, it is in use in physics and maybe it is worth including in packages. Original publication with algorithm included: http://iopscience.iop.org/article/10.1209/0295-5075/78/12001/meta

For free here: https://arxiv.org/abs/nucl-th/0702031

Screenshot of relevant part: • Interesting!! I haven't been working in this area for a few years now, but thanks for adding the response - surely it will help others looking for this somewhat niche topic. – Samik R Nov 13 '18 at 8:31