"SDP failed" with CSDP on local installations, but not on Macaulay2Web

[LegionMammal978]

If I run M2 locally on Ubuntu 24.04 (I've tried version 1.25.11 from ppa:macaulay2/macaulay2, version 1.25.11-82-ga66365a76b (development) from ppa:profzoom/macaulay2, and version 1.25.11-84-g432b335e4d (development) compiled from source), this command always fails:

i1 : needsPackage "SumsOfSquares"; QQ[x,y];

i3 : lowerBound(x^2+x*y+y^2+x+9, RoundTol=>infinity, Solver=>"CSDP")

o3 = (, MomentMatrix   null                    )
        GramMatrix     null                    
        Monomials      3x1 matrix over QQ[x..y]
        Status         SDP failed              

o3 : Sequence

However, if I run the same command on Macaulay2Web (version 1.25.11-1572-g6198e3f03b-dirty (vanilla)), it successfully gives me a lower bound:

i3 : net lowerBound(x^2+x*y+y^2+x+9, RoundTol=>infinity, Solver=>"CSDP")

o3 = (8.66667, MomentMatrix   3x3 matrix over RR_53           )
               GramMatrix     3x3 matrix over RR_53           
               Monomials      3x1 matrix over RR_53[x..y]     
               Parameters     1x1 matrix over RR_53           
               Status         SDP solved, primal-dual feasible

I am unable to get the corresponding commit (pzinn@6198e3f) to build locally to test it. But I've checked via get "!csdp || true" that both the local and web versions are running CSDP 6.2.0.

---

If I locally patch SumsOfSquares.m2 to print out the underlying optimize() call with printingPrecision = 0, I get:

needsPackage "SemidefiniteProgramming";
optimize(sdp(matrix{{4.499999999999998, .5, 0}, {.5, 1, .5}, {0, .5, 1}}, 1:matrix{{.7071067811865475, 0, 0}, {0, 0, 0}, {0, 0, 0}}, matrix{{.7071067811865475}}), Solver=>"CSDP")

which similarly returns SDP failed. However, if I replace 4.499999999999998 with 4.5, I get SDP solved, primal-dual feasible. So this issue may have to do with floating-point weirdness.

But weirdly, this isn't even a hard boundary for feasibility: replacing the first entry with 3.999999999999998 or 5.499999999999997 (corresponding to +8 or +11 in lowerBound()) also makes it succeed.

[LegionMammal978]

Actually, it seems I can get it to fail on the web version with lowerBound(x^2+x*y+y^2+x+2, RoundTol=>infinity, Solver=>"CSDP"). Maybe this instance is straddling the edge of feasibility, so that changing the constant can make the result go one way or the other within CSDP?

[LegionMammal978]

In particular, the input file for CSDP is:

*SDPA file generated by SemidefiniteProgramming.m2
1 =mdim
1 =nblocks
3
-.7071067811865475
0 1 1 1 -4.499999999999998
0 1 1 2 -.5
0 1 2 2 -1
0 1 2 3 -.5
0 1 3 3 -1
1 1 1 1 -.7071067811865475

CSDP succeeds with the default perturbobj=1, but hits Stuck at edge of primal feasibility, giving up. with the perturbobj=0 that M2 sets:

CSDP 6.2.0
Iter:  0 Ap: 0.00e+00 Pobj: -1.9500000e+02 Ad: 0.00e+00 Dobj:  0.0000000e+00 
Iter:  1 Ap: 1.00e+00 Pobj: -6.1470441e+01 Ad: 1.00e+00 Dobj:  2.7914225e+01 
Iter:  2 Ap: 1.00e+00 Pobj: -1.9064111e+01 Ad: 1.00e+00 Dobj:  3.2820554e+00 
Iter:  3 Ap: 1.00e+00 Pobj: -7.8910277e+00 Ad: 1.00e+00 Dobj: -2.3044861e+00 
Iter:  4 Ap: 1.00e+00 Pobj: -5.0977569e+00 Ad: 1.00e+00 Dobj: -3.7011215e+00 
Iter:  5 Ap: 1.00e+00 Pobj: -4.3994392e+00 Ad: 1.00e+00 Dobj: -4.0502804e+00 
Iter:  6 Ap: 9.00e-01 Pobj: -4.1899439e+00 Ad: 1.00e+00 Dobj: -4.1666667e+00 
Iter:  7 Ap: 3.04e-10 Pobj: -4.1899439e+00 Ad: 1.00e+00 Dobj: -4.1627871e+00 
Iter:  8 Ap: 9.00e-01 Pobj: -4.1689944e+00 Ad: 1.00e+00 Dobj: -4.1666667e+00 
Stuck at edge of primal feasibility, giving up. 
Failure: return code is 5 
Primal objective value: -4.1689944e+00 
Dual objective value: -4.1666667e+00 
Relative primal infeasibility: 0.00e+00 
Relative dual infeasibility: 9.54e-17 
Real Relative Gap: 2.49e-04 
XZ Relative Gap: 2.49e-04 
DIMACS error measures: 0.00e+00 0.00e+00 1.01e-16 0.00e+00 2.49e-04 2.49e-04

Meanwhile, Macaulay2Web has a very slightly different input file, which allows CSDP to succeed with perturbobj=0:

*SDPA file generated by SemidefiniteProgramming.m2
1 =mdim
1 =nblocks
3
-.7071067811865476
0 1 1 1 -4.499999999999998
0 1 1 2 -.5
0 1 2 2 -1
0 1 2 3 -.5
0 1 3 3 -1
1 1 1 1 -.7071067811865475

[d-torrance]

One problem with M2 interfacing with software using floating-point arithmetic is that our binary -> decimal conversion function for printing isn't great (we do something like truncate at 16 decimal places, I think). So the numbers we're sending to CSDP may not actually be same numbers that we think we're supposed to be sending. (I'm not sure why it would be different in webapp mode, though...)

I'd really like to implement a good binary -> decimal algorithm one of these days (Dragon4, Grisu, Errol, Ryu, Dragonbox -- there are lots of choices).

Related: #708