What algorithm does CVXPY actually use to solve semidefinite programs with the constraints of the form $\sum\limits_i E_iXE_i^T \succ B$?

Question

_{Crossposted on Mathematics SE}

---

CVXPY is a famous software as a solver for optimization problems. Nowadays, I use it to run a program presented in a paper, the Example 7.1, and the program runs as follows (I haven't found out how to insert my codes here, otherwise I would post them) all matrices are $4 \times 4$: $O_0$ and $O_3$ are varibles, the target function is $\min\text{tr}(O_0+O_3)$ and the constraints are

$$O_0 \succeq 0, \qquad O_3 \succeq 0, \qquad O_0 \preceq I, \qquad O_3 \preceq I$$

$$\left( \begin{matrix}0&0.707&0&0.707\\0&0&0&0\\0&0.707&0&-0.707\\0&0&0&0 \end{matrix}\right)O_0\left( \begin{matrix}0&0&0&0\\0.707&0&0.707&0\\0&0&0&0\\0.707&0&-0.707&0 \end{matrix}\right)+\left(\begin{matrix}0&0&0&0\\0&1&0&0\\0&0&0&0\\0&0&0&1 \end{matrix}\right)O_3\left(\begin{matrix}0&0&0&0\\0&1&0&0\\0&0&0&0\\0&0&0&1 \end{matrix}\right)\succeq\left(\begin{matrix}1&0&0&0\\0&0&0&0\\0&0&0&0\\0&0&0&0 \end{matrix}\right)\\ \left( \begin{matrix}0&0.707&0&0.707\\0&0&0&0\\0&0.707&0&-0.707\\0&0&0&0 \end{matrix}\right)O_0\left( \begin{matrix}0&0&0&0\\0.707&0&0.707&0\\0&0&0&0\\0.707&0&-0.707&0 \end{matrix}\right)-O_0+\left(\begin{matrix}0&0&0&0\\0&1&0&0\\0&0&0&0\\0&0&0&1 \end{matrix}\right)O_3\left(\begin{matrix}0&0&0&0\\0&1&0&0\\0&0&0&0\\0&0&0&1 \end{matrix}\right)\succeq 0$$

The output of the program is: the optimal value is $1.9999988116689815$, and

$$O_0=O_3=\left(\begin{matrix}0&0&0&0\\0&0.5&0&0.5\\0&0&0&0\\0&0.5&0&0.5 \end{matrix}\right)$$

It is not hard to see that the constraints are all in the form of $$\sum\limits_i\sum\limits_j A_{ij} X_i A_{ij}^T \preceq \text{(or} \succeq\text{)} B$$ where $X_i$ are variables. However, for SDP algorithms(such as interior point method and cutting plane method) what we are dealing with are SDP is in the form of the standard or the dual form.However the above CVXPY solver can directly give the output without converting to standard form or dual form, so I wonder what algorithm does CVXPY actually take. I have read the guidebook of CVXPY but haven't found out where it is mentioned.I am interested in this problem because I once studied such semidefinite programs.

Answer 1

Cvxpy cannot solve SDPs by it itself. It feeds the problem into an optimizer such as Mosek. Therefore, you should consult the documentation of the optimizer you are using.

Btw it is trivial to convert your SDP to any standard form. That is precisely the thing that cvxpy does.

Answer 2

I assume your question is "which optimizer did cvxpy choose for my problem", as cvxpy can use a number of optimizers that it selects based on the problem, and indeed the output doesn't specify which one was used.

To find out more, you can run the solver with the verbose argument. That will tell you the algorithm, and it even sometimes gives links to the publications and authors etc...

problem.solve(verbose=True)

Here's an example of the relevant part of the logs for a quadratic optimization but you'll need to run it for your problem:

-------------------------------------------------------------------------------
                                Numerical solver                               
-------------------------------------------------------------------------------
(CVXPY) Aug 03 01:01:33 PM: Invoking solver OSQP  to obtain a solution.
-----------------------------------------------------------------
           OSQP v0.6.3  -  Operator Splitting QP Solver
              (c) Bartolomeo Stellato,  Goran Banjac
        University of Oxford  -  Stanford University 2021
-----------------------------------------------------------------
problem:  variables n = 6, constraints m = 1
          nnz(P) + nnz(A) = 21
settings: linear system solver = qdldl,
          eps_abs = 1.0e-05, eps_rel = 1.0e-05,
          eps_prim_inf = 1.0e-04, eps_dual_inf = 1.0e-04,
          rho = 1.00e-01 (adaptive),
          sigma = 1.00e-06, alpha = 1.60, max_iter = 10000
          check_termination: on (interval 25),
          scaling: on, scaled_termination: off
          warm start: on, polish: on, time_limit: off

iter   objective    pri res    dua res    rho        time
   1  -2.0032e-03   0.00e+00   1.63e-03   1.00e-01   7.09e-05s
  25  -3.1300e-03   0.00e+00   7.70e-09   1.00e-01   9.38e-05s
plsh  -3.1300e-03   0.00e+00   4.69e-19   --------   1.36e-04s

status:               solved
solution polish:      successful
number of iterations: 25
optimal objective:    -0.0031
run time:             1.36e-04s
optimal rho estimate: 1.00e-06
```