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:


*

*Cholesky factorization

*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$.
 A: 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.
A: 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: 

