3
$\begingroup$

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]
$\endgroup$
1
  • 1
    $\begingroup$ 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 $\endgroup$ Sep 21 '21 at 23:46
2
$\begingroup$

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]

enter image description here

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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