# Cholesky with complete pivoting?

I'm looking for an implementation of Cholesky decomposition with complete pivoting. It is the lev3pchol.f function added in Lapack 3.2 (working note)

It has an advantage over regular Cholesky of working on singular matrices and can be used to upper bound matrix rank, wondering if someone has ported it already.

Higher level description is in Higham's note

• There was a talk at the 2003 developer conference about directly calling precompiled Netlib numerical libraries: https://library.wolfram.com/infocenter/Conferences/4891/. As far as I know, Netlib has Lapack:  netlib.org. Maybe it will helpful to you. Aug 2, 2022 at 14:43

If matrix $$A$$ comes in Gram-factored form $$A=X^T X$$, then pivoted Cholesky can be done through QR like in this answer. QRDecomposition has a Pivoted->True option.

(* given X produces R,p such that p\[Transpose].X\[Transpose].X.p==R\
\[Transpose].R *)
PivotedCholeskyDecompositionOfXTX[X_] := Module[{q, R, p},
{q, R, p} = QRDecomposition[X, Pivoting -> True];
{Map[If[Negative[#], -1, 1] &, Diagonal[R]]*R, p}
];

(* from https://mathematica.stackexchange.com/a/271404/217 *)
CholeskyDecompositionOfXTX[X_] :=
With[{R = Last[QRDecomposition[X]]},
Map[If[Negative[#], -1, 1] &, Diagonal[R]]*R];

d = 5;
numSamples = 5;
evals = Table[1/i, {i, 1, d}];
sigma = DiagonalMatrix[evals];
dist = MultinormalDistribution[sigma];
X = RandomVariate[dist, numSamples];

R = CholeskyDecompositionOfXTX[X];
Print["Unpivoted diagonal is decreasing: ",
Diagonal[R] == ReverseSort@Diagonal@R]  (* False *)
{R, p} = PivotedCholeskyDecompositionOfXTX[X];
Print["Pivoted reconstruction works: ",
Norm[p\[Transpose] . X\[Transpose] . X . p - R\[Transpose] . R] <
10^-10]; (* True *)
Print["Pivoted diagonal is decreasing: ",
Diagonal[R] == ReverseSort@Diagonal@R] (* True *)


If it's not factored, then one approach is this (needs Mathematica implementation, one issue is that QRDecomposition gives reduced decomposition, but that approach needs the full version)

• Sep 2, 2022 at 15:03