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

Question

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]

Answer 1

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]