Cholesky-like (rank-revealing?) decomposition on singular matrices

For a singular Gram matrix $$B$$ I need to find some matrix $$A$$ such that $$AA'=B$$. I can add a small multiple of identity matrix, but this could break upstream analysis (application is c-approximation to graph cut here).

Searching around I see "Rank-revealing VSV decomposition" (Ch 2 of UTV expansion pack) and modified Cholesky decomposition implemented here. However nothing readily available in Mathematica, any tips?

Example of failure:

g = RandomGraph[BernoulliGraphDistribution[20, .5]];
setupSDP[g_Graph] := Block[
Sequence[{n, activeVars, objective, m, constraint, A, a, b},
n = Length[
VertexList[g]]; activeVars = ReplaceAll[
EdgeList[g], UndirectedEdge[
Pattern[a,
Blank[]],
Pattern[b,
Blank[]]] -> A[a, b]]; objective = -Total[
Map[(1 - #)/2& , activeVars]]; m = Array[A, {n, n}];
m = UpperTriangularize[m] - DiagonalMatrix[
Diagonal[m]]; m = m + Transpose[m] + IdentityMatrix[n];
constraint =
VectorGreaterEqual[{m, 0}, {"SemidefiniteCone", n}]; {
objective, constraint, m}]];
{objective, constraint, matrix} = setupSDP[g];
solution =
matrix /.
SemidefiniteOptimization[objective, constraint,
Variables[matrix]];
CholeskyDecomposition[solution]

• Could use the resource function RationalCholeskyDecomposition: In[272]:= {ssqrt, d} = ResourceFunction["RationalCholeskyDecomposition"][solution]; Max[Abs[ssqrt . DiagonalMatrix[d] . Transpose[ssqrt] - solution]] Out[273]= 2.22045*10^-16 Commented Sep 21, 2021 at 23:46

The ResourceFunction["RationalCholeskyDecomposition"][solution] seems to do the job, as you can see from the rank-2 solution of SDP relaxation to maxcut of sun graph below

(* Given a Gram matrix of reduced rank, produce low rank square root *)

lowRankRepresentation[gram_] := (
RationalCholesky =
ResourceFunction["RationalCholeskyDecomposition"];
{sqrt, d} = RationalCholesky[gram];
maxEig = Max@Eigenvalues[gram];
cutoff = 10^-7*maxEig;
rank =
Count[RationalCholesky[solution] // Last, _?(# > cutoff &)];
Transpose[Transpose[sqrt][[;; rank]]]
);

setupSDP[g_Graph] :=
Block[{n, activeVars, objective, m, constraint, A, a, b},
n = Length[VertexList[g]];
activeVars = EdgeList[g] /. a_ \[UndirectedEdge] b_ -> A[a, b];
objective = -Total[((1 - #1)/2 &) /@ activeVars];
m = Array[A, {n, n}];
m = UpperTriangularize[m] - DiagonalMatrix[Diagonal[m]];
m = m + Transpose[m] + IdentityMatrix[n]; constraint = m
\!$$\*UnderscriptBox[\(\[VectorGreaterEqual]$$,
TemplateBox[{"n"},
"SemidefiniteConeList"]]\) 0; {objective, constraint, m}
];

createPlotHelpers := (
n = Length[VertexList[g]];
{v1, v2} =
Orthogonalize[RandomVariate[NormalDistribution[], {2, n}]];
vfArrow[{xc_, yc_}, name_, {w_, h_}] := Block[{},
vec = xs[[name]];
{dirx, diry} = Normalize[vec]/4.0;
{Black, Arrow[{{xc - dirx, yc - diry}, {xc + dirx, yc + diry}}],
ColorData["Crayola"]["NavyBlue"], Dashed,
Circle[{xc, yc}, Max[w, h], {0, 2 \[Pi]}]}];

vfLabel[{xc_, yc_}, name_, {w_, h_}] := Block[{},
vec = xs[[name]]; {dirx, diry} = Normalize[vec]/2.5; {Gray,
Text[Style[Sign[dir . vec], Large], {xc, yc}],
ColorData["Crayola"]["NavyBlue"], Dashed,
Circle[{xc, yc}, Max[w, h], {0, 2 \[Pi]}]}];
graphLabels[g_, dd_] := (dir = dd;
GraphPlot[g, VertexShapeFunction -> vfLabel,
VertexSize -> {"Scaled", 0.15}]);
; graphArrows[g_] :=
GraphPlot[g, VertexShapeFunction -> vfArrow,
VertexSize -> {"Scaled", 0.15}];

);

g = GraphData[{"Sun", 4}];
{objective, constraint, matrix} = setupSDP[g];
solution =
matrix /.
SemidefiniteOptimization[objective, constraint,
Variables[matrix]];
xs = lowRankRepresentation[solution];

createPlotHelpers
graphArrows[g]