Multiprecision and Python 3.11 support for SDPA (#2097)

* Extends our SDPA interface to use either ``sdpa-python`` or ``sdpa-multiprecision``.

* Add an ill-posed LP (test_lp_7) to solver_test_helpers.py to test SDPA's multiprecision support.

---------

Co-authored-by: Riley Murray <rileyjmurray@users.noreply.github.com>

continuous_integration/install_dependencies.sh

@@ -32,11 +32,22 @@ fi
 
 if [[ "$PYTHON_VERSION" == "3.11" ]]; then
   python -m pip install gurobipy clarabel osqp
+  if [[ "$RUNNER_OS" == "Windows" ]]; then
+    # SDPA with OpenBLAS backend does not pass LP5 on Windows
+    python -m pip install sdpa-multiprecision
+  else
+    python -m pip install sdpa-python
+  fi
 # Python 3.8 on Windows will uninstall NumPy 1.16 and install NumPy 1.24 without the exception.
 elif [[ "$RUNNER_OS" == "Windows" ]] && [[ "$PYTHON_VERSION" == "3.8" ]]; then
   python -m pip install gurobipy clarabel osqp
 else
-  python -m pip install "ortools>=9.3,<9.5" coptpy cplex sdpa-python diffcp gurobipy xpress clarabel sdpa-python
+  python -m pip install "ortools>=9.3,<9.5" coptpy cplex diffcp gurobipy xpress clarabel
+  if [[ "$RUNNER_OS" == "Windows" ]]; then
+    python -m pip install sdpa-multiprecision
+  else
+    python -m pip install sdpa-python
+  fi
 fi
 
 # cylp has wheels for all versions 3.7 - 3.10, except for 3.7 on Windows

cvxpy/reductions/solvers/conic_solvers/sdpa_conif.py

@@ -84,52 +84,6 @@ def apply(self, problem):
         tuple
             (dict of arguments needed for the solver, inverse data)
         """
-
-        # CVXPY represents cone programs as
-        #     (P) min_x { c^T x : A x + b \in K } + d,
-        # or, using Dualize
-        #     (D) max_y { -b^T y : c - A^T y = 0, y \in K^* } + d.
-
-        # SDPAP takes a conic program in CLP form:
-        #     (P) min_x { c^T x : A x - b \in J, x \in K }
-        # or
-        #     (D) max_y { b^T y : y \in J^*, c - a^T y \in K^*}
-
-        # SDPA takes a conic program in SeDuMI form:
-        #     (P) min_x { c^T x : A x - b = 0, x \in K }
-        # or
-        #     (D) max_y { b^T y : c - A^T y \in K^*}
-
-        # SDPA always takes an input in SeDuMi form
-        # SDPAP therefore converts the problem from CLP form to SeDuMi form which
-        # involves getting rid of constraints involving cone J (and J^*) of the CLP form
-
-        # We have two ways to solve using SDPA
-
-        # 1. CVXPY (P) -> CLP (P), by
-        #       - flipping sign of b
-        #       - setting J of CLP (P) to K of CVXPY (P)
-        #       - setting K of CLP (P) to a free cone
-        #
-        #       (a) then CLP (P)-> SeDuMi (D) if user specifies `convMethod` option to `clp_toLMI`
-        #       (b) then CLP (P)-> SeDuMi (P) if user specifies `convMethod` option to `clp_toEQ`
-        #
-        # 2. CVXPY (P) -> CVXPY (D), by
-        #       - using Dualize
-        #       - then CVXPY (D) -> SeDuMi (P) by
-        #           - setting c of SeDuMi (P) to -b of CVXPY (D)
-        #           - setting b of SeDuMi (P) to c of CVXPY (D)
-        #           - setting x of SeDuMi (P) to y of CVXPY (D)
-        #           - setting K of SeDuMi (P) to K^* of CVXPY (D)
-        #           - transposing A
-
-        # 2 does not give a benefit over 1(a)
-        #   1 (a) and 2 both flip the primal and dual between CVXPY and SeDuMi
-        # 1 (b) does not flip primal and dual throughout, but requires `clp_toEQ` to be implemented
-        #   TODO: Implement `clp_toEQ` in `sdpa-python`
-        # Thereafter, 1 will allows user to choose between solving as primal or dual
-        # by specifying `convMethod` option in solver options
-
         data = {}
         inv_data = {self.VAR_ID: problem.x.id}
 
@@ -179,6 +133,28 @@ def invert(self, solution, inverse_data):
             return failure_solution(status)
 
     def solve_via_data(self, data, warm_start: bool, verbose: bool, solver_opts, solver_cache=None):
+        """
+        CVXPY represents cone programs as
+            (P) min_x { c^T x : A x + b \in K } + d
+
+        SDPA Python takes a conic program in CLP format:
+            (P) min_x { c^T x : A x - b \in J, x \in K }
+
+        CVXPY (P) -> CLP (P), by
+            - flipping sign of b
+            - setting J of CLP (P) to K of CVXPY (P)
+            - setting K of CLP (P) to a free cone
+
+        CLP format is a generalization of the SeDuMi format. Both formats are explained at
+        https://sdpa-python.github.io/docs/formats/
+
+        Internally, SDPA Python will reduce CLP form to SeDuMi dual form using `clp_toLMI`.
+        In SeDuMi format, the dual is in LMI form. In SDPA format, the primal is in LMI form.
+        The backend (i.e. `libsdpa.a` or `libsdpa_gmp.a`) uses the SDPA format.
+
+        For details on the reverse relationship between SDPA and SeDuMi formats, please see
+        https://sdpa-python.github.io/docs/formats/sdpa_sedumi.html
+        """
         import sdpap
         from scipy import matrix
 
@@ -200,40 +176,13 @@ def solve_via_data(self, data, warm_start: bool, verbose: bool, solver_opts, sol
         x, y, sdpapinfo, timeinfo, sdpainfo = sdpap.solve(
             A, -matrix(b), matrix(c), K, J, solver_opts)
 
-        # This should be set according to the value of `#define REVERSE_PRIMAL_DUAL`
-        # in `sdpa` (not `sdpa-python`) source.
-        # By default it's enabled and hence, the primal problem in SDPA takes
-        # the LMI form (opposite that of SeDuMi).
-        # If, while building `sdpa` from source you disable it, you need to change this to `False`.
-        reverse_primal_dual = True
-        # By disabling `REVERSE_PRIMAL_DUAL`, the primal-dual pair
-        # in SDPA becomes similar to SeDuMi, i.e. the form we want (see long note in `apply` method)
-        # However, with `REVERSE_PRIMAL_DUAL` disabled, we have some
-        # accuracy related issues on unit tests using the default parameters.
-
-        REVERSE_PRIMAL_DUAL = {
-            "noINFO": "noINFO",
-            "pFEAS": "pFEAS",
-            "dFEAS": "dFEAS",
-            "pdFEAS": "pdFEAS",
-            "pdINF": "pdINF",
-            "pFEAS_dINF": "pINF_dFEAS",  # invert flip within sdpa
-            "pINF_dFEAS": "pFEAS_dINF",  # invert flip within sdpa
-            "pdOPT": "pdOPT",
-            "pUNBD": "dUNBD",  # invert flip within sdpa
-            "dUNBD": "pUNBD"  # invert flip within sdpa
-        }
-
-        if (reverse_primal_dual):
-            sdpainfo['phasevalue'] = REVERSE_PRIMAL_DUAL[sdpainfo['phasevalue']]
-
         solution = {}
-        solution[s.STATUS] = self.STATUS_MAP[sdpainfo['phasevalue']]
+        solution[s.STATUS] = self.STATUS_MAP[sdpapinfo['phasevalue']]
 
         if solution[s.STATUS] in s.SOLUTION_PRESENT:
             x = x.toarray()
             y = y.toarray()
-            solution[s.VALUE] = sdpainfo['primalObj']
+            solution[s.VALUE] = sdpapinfo['primalObj']
             solution[s.PRIMAL] = x
             solution[s.EQ_DUAL] = y[:dims['f']]
             solution[s.INEQ_DUAL] = y[dims['f']:]

cvxpy/tests/solver_test_helpers.py

@@ -256,6 +256,31 @@ def lp_6() -> SolverTestHelper:
     return sth
 
 
+def lp_7() -> SolverTestHelper:
+    """
+    An ill-posed problem to test multiprecision ability of solvers.
+
+    This test will not pass on CVXOPT (as of v1.3.1) and on SDPA without GMP support.
+    """
+    n = 50
+    a = cp.Variable((n+1))
+    delta = cp.Variable((n))
+    b = cp.Variable((n+1))
+    objective = cp.Minimize(cp.sum(cp.pos(delta)))
+    constraints = [
+        a[1:] - a[:-1] == delta,
+        a >= cp.pos(b),
+    ]
+    con_pairs = [(constraints[0], None),
+                 (constraints[1], None)]
+    var_pairs = [(a, None),
+                 (delta, None),
+                 (b, None)]
+    obj_pair = (objective, 0.)
+    sth = SolverTestHelper(obj_pair, var_pairs, con_pairs)
+    return sth
+
+
 def qp_0() -> SolverTestHelper:
     # univariate feasible problem
     x = cp.Variable(1)
@@ -1012,6 +1037,15 @@ def test_lp_6(solver, places: int = 4, duals: bool = True, **kwargs) -> SolverTe
             sth.check_dual_domains(places)
         return sth
 
+    @staticmethod
+    def test_lp_7(solver, places: int = 4, duals: bool = True, **kwargs) -> SolverTestHelper:
+        sth = lp_7()
+        import sdpap
+        if sdpap.sdpacall.sdpacall.get_backend_info()["gmp"]:
+            sth.solve(solver, **kwargs)
+            sth.verify_objective(places)
+        return sth
+
     @staticmethod
     def test_mi_lp_0(solver, places: int = 4, **kwargs) -> SolverTestHelper:
         sth = mi_lp_0()

cvxpy/tests/test_conic_solvers.py

@@ -847,11 +847,13 @@ def test_sdpa_lp_3(self) -> None:
     def test_sdpa_lp_4(self) -> None:
         StandardTestLPs.test_lp_4(solver='SDPA')
 
-    @unittest.skip('Known limitation of SDPA for degenerate LPs.')
     def test_sdpa_lp_5(self) -> None:
         # this also tests the ability to pass solver options
         StandardTestLPs.test_lp_5(solver='SDPA',
-                                  gammaStar=0.86, epsilonDash=8.0E-6, betaStar=0.18, betaBar=0.15)
+                                  betaBar=0.1, gammaStar=0.8, epsilonDash=8.0E-6)
+
+    def test_sdpa_lp_7(self) -> None:
+        StandardTestLPs.test_lp_7(solver='SDPA')
 
     def test_sdpa_socp_0(self) -> None:
         StandardTestSOCPs.test_socp_0(solver='SDPA')

cvxpy/tests/test_problem.py

@@ -283,7 +283,7 @@ def test_verbose(self) -> None:
                 # Don't test GLPK because there's a race
                 # condition in setting CVXOPT solver options.
                 if solver in [cp.GLPK, cp.GLPK_MI, cp.MOSEK, cp.CBC,
-                              cp.SCIPY, cp.SDPA, cp.COPT]:
+                              cp.SCIPY, cp.COPT]:
                     continue
                 sys.stdout = StringIO()  # capture output
 

doc/source/tutorial/advanced/index.rst

@@ -460,7 +460,7 @@ The table below shows the types of problems the supported solvers can handle.
 +----------------+----+----+------+-----+-----+-----+-----+
 | `CVXOPT`_      | X  | X  | X    | X   |     |     |     |
 +----------------+----+----+------+-----+-----+-----+-----+
-| `SDPA`_        | X  | X  | X    | X   |     |     |     |
+| `SDPA`_ ***    | X  | X  | X    | X   |     |     |     |
 +----------------+----+----+------+-----+-----+-----+-----+
 | `SCS`_         | X  | X  | X    | X   | X   | X   |     |
 +----------------+----+----+------+-----+-----+-----+-----+
@@ -475,6 +475,8 @@ The table below shows the types of problems the supported solvers can handle.
 
 (**) Except mixed-integer SDP.
 
+(***) Multiprecision support is available on SDPA if the appropriate SDPA package is installed. With multiprecision support, SDPA can solve your problem with much smaller `epsilonDash` and/or `epsilonStar` parameters. These parameters must be manually adjusted to achieve the desired degree of precision. Please see the solver website for details. SDPA can also solve some ill-posed problems with multiprecision support.
+
 Here EXP refers to problems with exponential cone constraints. The exponential cone is defined as
 
     :math:`\{(x,y,z) \mid y > 0, y\exp(x/y) \leq z \} \cup \{ (x,y,z) \mid x \leq 0, y = 0, z \geq 0\}`.