Cholesky with complete pivoting?

Question

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

Answer 1

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)