Discrepances between CLP format and SDPA format

[ganglesmunne]

Hi, I am trying to solve an SDP but I come out with a solution which is a bit weird. The variable "d.feas.error" and "primalError" are very different in order of magnitude:

---------- SDPAP Start ----------
SDPA Multiprecision start at [Fri Jun 14 15:42:33 2024]
Schur computation : DENSE 
Converted to SDPA internal data / Starting SDPA main loop
   mu      thetaP  thetaD  objP      objD      alphaP  alphaD  beta 
 0 1.0e+08 1.0e+00 1.0e+00 +0.00e+00 +0.00e+00 9.0e-01 9.0e-01 2.00e-01
 1 1.4e+07 1.0e-01 1.0e-01 +1.26e+04 -1.90e-01 9.0e-01 9.0e-01 2.00e-01
 2 1.8e+06 1.0e-02 1.0e-02 +1.84e+04 -6.75e-03 9.0e-01 9.0e-01 2.00e-01
 3 2.5e+05 1.0e-03 1.0e-03 +2.51e+04 -1.20e-02 9.0e-01 9.0e-01 2.00e-01
 4 3.4e+04 1.0e-04 1.0e-04 +3.41e+04 -1.15e-02 9.0e-01 9.0e-01 2.00e-01
 5 4.8e+03 1.0e-05 1.0e-05 +4.56e+04 -1.13e-02 9.0e-01 9.0e-01 2.00e-01
 6 6.8e+02 1.0e-06 1.1e-06 +5.23e+04 -1.13e-02 9.0e-01 8.9e-01 2.00e-01
 7 1.0e+02 1.0e-07 1.1e-07 +2.67e+04 -1.11e-02 8.5e-01 9.0e-01 2.00e-01
 8 2.1e+01 1.5e-08 1.1e-08 +6.68e+03 -1.18e-02 9.0e-01 1.0e+00 2.00e-01
 9 3.1e+00 1.5e-09 1.6e-56 +1.15e+03 -1.22e-02 9.0e-01 9.0e-01 2.00e-01
10 4.5e-01 1.5e-10 5.6e-54 +1.86e+02 -1.22e-02 9.0e-01 3.7e+01 2.00e-01
11 6.8e-01 1.5e-11 5.4e-52 +2.87e+01 -1.20e-02 9.0e-01 9.0e-01 2.00e-01
12 1.0e-01 1.5e-12 5.3e-53 +1.82e+01 -1.23e-02 9.0e-01 9.6e+00 2.00e-01
13 4.7e-02 1.6e-13 3.1e-51 +4.07e+00 -6.02e-03 9.0e-01 9.0e-01 2.00e-01
14 7.1e-03 1.6e-14 4.3e-52 +1.48e+00 -7.71e-03 8.9e-01 2.0e+00 2.00e-01
15 1.4e-03 1.7e-15 3.6e-52 +3.12e-01 -3.75e-03 8.7e-01 1.7e+00 2.00e-01
16 3.3e-04 2.3e-16 2.8e-52 +7.00e-02 -1.60e-03 7.4e-01 7.6e-01 2.00e-01
17 1.1e-04 6.2e-17 4.9e-46 +2.36e-02 -6.78e-04 5.6e-01 6.4e-01 2.00e-01
18 6.3e-05 2.7e-17 3.5e-45 +1.17e-02 -1.90e-04 6.5e-01 7.4e-01 2.00e-01
19 2.9e-05 9.5e-18 3.1e-45 +5.00e-03 -2.91e-05 8.4e-01 8.4e-01 2.00e-01
20 6.3e-06 1.5e-18 3.8e-45 +1.36e-03 -3.03e-05 8.0e-01 9.4e-01 2.00e-01
21 1.7e-06 3.0e-19 5.5e-45 +3.82e-04 +4.17e-07 8.2e-01 8.7e-01 1.00e-01
22 3.9e-07 3.0e-19 1.9e-45 +7.59e-05 +4.56e-06 6.8e-01 7.2e-01 1.00e-01
23 1.4e-07 3.0e-19 6.4e-46 +2.75e-05 +6.27e-06 7.3e-01 5.2e-01 1.00e-01
24 7.2e-08 3.0e-19 4.6e-46 +1.23e-05 +9.60e-06 1.7e-01 1.6e-01 1.00e-01
25 6.1e-08 3.0e-19 4.0e-46 +1.11e-05 +1.11e-05 1.7e-01 1.6e-01 1.00e-01

phase.value = pdOPT     
   Iteration = 25
          mu = 6.0909024026615191e-08
relative gap = 3.0473987620798781e-68
         gap = 1.2151350293309731e-04
      digits = 6.2562207098598208e+01
objValPrimal = 1.1120809091905244e-05
objValDual   = 1.1120809091905244e-05
**p.feas.error = 9.7209123950282212e-10
d.feas.error = 5.1572966133778820e-38**
relative eps = 8.6361685550944446e-78
total time   = 470.330
  main loop time = 470.329801
      total time = 470.329998
file  check time = 0.000000
file change time = 0.000197
file   read time = 0.000000
Converting optimal solution to CLP format
SDPA Multiprecision end at [Fri Jun 14 15:50:24 2024]
Start: getCLPresult
Making result infomation...
~/.conda/envs/lab/lib/python3.8/site-packages/scipy/sparse/linalg/_eigen/arpack/arpack.py:1272: RuntimeWarning: k >= N - 1 for N * N square matrix. Attempting to use scipy.linalg.eig instead.
  warnings.warn("k >= N - 1 for N * N square matrix. "
========================================
 SDPAP: Result
========================================
    SDPA.phase = pdOPT
     iteration = 25
    convMethod = LMI
     frvMethod = split
  domainMethod = none
   rangeMethod = none
     primalObj = -1.1120809091905243e-05
       dualObj = -1.1120809091905243e-05
    dualityGap = +2.0000000000000000e+00
   **primalError = +5.8151092469806038e+01
     dualError = +3.6379788070917130e-12**
   convertTime = 0.000627
     solveTime = 470.611907
retrievingTime = 0.000002
     totalTime = 470.612663
---------- SDPAP End ----------

According to the documentation,

SPA-format -> https://citeseerx.ist.psu.edu/document?repid=rep1&type=pdf&doi=146eeab7f3c917c092cbb0062029e92f7ac9e4d8
(pag 14)

CLP-format -> https://master.dl.sourceforge.net/project/sdpa/sdpa-p/sdpap_manual.pdf?viasf=1 (pag 19)

shoudn't these two values be comparable?

Due to this discrepancy I am not sure if I can trust the solution. Is it possible that there is a problem in the conversion from SDPA- to CLP-format?

Thanks!

[usamamuneeb]

Hello, it's the primal and dual objectives which should converge close to each other which is happening in your case. Regarding primal and dual feasibility errors, ideally they should be 0, which is also true in your case. They will never be exactly zero but both e-10 or e-38 are practically considered zero. I am looking at the values in your upper set of highlighting which is the C++ solver (SDPA Multiprecision) output.

The lower highlight represents SDPA for Python output, which, as you pointed out solves a CLP. The feasibility errors are not converted, but recomputed in SDPAP (from the SDPA to CLP converted solutions).

[usamamuneeb]

@ganglesmunne, it seems like my initial analysis was not the case here. I just ran a few examples using the regular (sdpa-python) backend and the feasibility errors are similar between C++ solver and SDPA for Python (while sdpa-multiprecision backend gives orders of magnitude smaller feasibility errors in the C++ solver's log for the same problem).

Is it possible for you to run your example without multiprecision support? If yes, I suspect you should also see similar feasibility errors between C++ solver and SDPAP. In summary, when using the multiprecision backend, the C++ solver has to its advantage arbitrary precision support when calculating feasibility errors.

It also seems like the page 15-16 of SDPAP manual reports incorrect dual feasibility error (either due to a bug or due to numerical limitations in the Python code doing feasibility error calculation). It perhaps should not be zero. I will investigate this issue further. Thank you for raising this issue!

[ganglesmunne]

You are welcome! My example without multiprecision support does not go beyond the first iteration. Even though, the program gives me the feasibility errors and they indeed coincide when I am not using the backend.

Does this mean that the feasibility errors from the multiprecision backend are the ones to take into account?

Thank you very much for the feedback!

[usamamuneeb]

Hi @ganglesmunne, yes, you should use the feasibility errors of the C++ solver (multiprecision in this case), because they are computed using the same precision as that of the solver. Here’s my analysis of the problem:

If you use multiprecision backend, the feasibility errors computed by SDPAP will be orders of magnitude larger than those computed by C++ solver (once the problem is solved by multiprecision backend, the solution will be converted to double precision when it goes back to SDPAP). Moreover, if your problem has a very small feasibility set (which I believe is also your case), SDPAP may show you feasibility errors large enough to cause early termination had such large feasibility errors been detected during C++ solver’s main loop (i.e. if the C++ solver did not have multiprecision support).

If your SDPAP and C++ solver both have the same precision (double), then barring the primal dual flip, you should ideally see exactly the same feasibility errors between SDPAP and C++ solver. While CLP format is more general, in order to do a max between violations of different types of constraints, one needs to have them computed in a comparable or equivalent way.

If you use regular (sdpa-python) backend, you’ll notice that the primal error of SDPAP is correctly computed as the dual error of C++ solver. The SDPAP dual error is the problematic part which I have found to be reliant on eigenvalue computation (which has numerical challenges). As a solution, I have decided to disable (re)computation of feasibility errors in the Python code.

[ganglesmunne]

Thank you very much for you explanation, it was really helpful!

[usamamuneeb]

You are very welcome. I’ll see if I can add a secondary feasibility error computation directly inside SDPA Multiprecision code using double precision so the user also knows feasibility errors after precision reduction.

If precision reduction to double significantly increases feasibility error, the user can be warned and potential fixes suggested (relax optimality criterion epsilonStar or relax the problem by removing some constraints).

This phenomenon however is not uncommon when working with mixed precision.