1
$\begingroup$

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

$\endgroup$
1
  • 2
    $\begingroup$ 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. $\endgroup$ Aug 2, 2022 at 14:43

1 Answer 1

1
$\begingroup$

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)

$\endgroup$
1

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.