Storage complexity of SDP solver SCS

Question

This is a follow up question to this question.

Consider the following SDP in standard form:

\begin{align} &\min_{X\in S^n, X>0} \operatorname{tr}(AX)\\ &\mbox{subject to}\; \operatorname{tr}(B_iX)\leq b_i, i=1,\dots,m\\ \end{align}

Here, $A\in S^n,B_i \in S^n$, $b_i\in \mathbb{R}$ and $m=O(n)$.

Can anyone please explain the expected storage complexity of this problem in the general case?

My understanding is that first order solver like SCS will form a solve a $n^2$ sized linear system using gradient information only, and hence the storage complexity will be $O(n^4)$ which reflects the size of gradient matrix which will be $n^2\times n^2$.

Is this understanding correct? My confusion is regarding the possibility of dual problem being smaller size. Is that taken care of by CVX automatically? For example, if $n=1000$ and $m=1000$, it seems like primal form as mentioned above would form the smaller problem.

---

Edit After Johan's Answer: I am currently using CVX. I tried an example similar to what you gave. Using notation in my post above, I am using $n=100,m=10$. $X\in S^n$. If I submit my problem in format above, I get following output:

---

Calling SCS 1.1.7: 5050 variables, 10 equality constraints

Lin-sys: sparse-indirect, nnz in A = 50500, CG tol ~ 1/iter^(2.00) eps = 1.35e-06, alpha = 1.50, max_iters = 10000, normalize = 1, scale = 1.00 Variables n = 10, constraints m = 5050 Cones: sd vars: 5050, sd blks: 1 Setup time: 3.02e-03s

---

Now if I submit my problem by forming a dual of this, I get the message that CVX will solve the dual of my dual (i.e. the original problem), and I get exactly same number of constraints and variables as above.

So my question is whether CVX is really solving the way you have described ? What is the expected space complexity given the above variables ?

Answer 1

No, your understanding is not correct. By introducing a slack $s\in R^{m}_+$ you can write the constraints in the form $\operatorname{tr}(B_iX) + s_i = b_i$. You now have a standard semidefinite program, in primal form, over the product of a semidefinite cone of size $n$, and an LP cone of size $m$. The linear system will be a linear system of size $m \times m$ as there are $m$ equalities in the primal ($m$ variables in the dual). A very bad model would be to do what you implicitly hint, i.e., to introduce the elements of $X$ as variables, leading to a linear system of size $O(n^2)$. When solving a primal-dual pair in semidefinite programming, you never explicitly solve for the elements of the primal cone, the matrix as a whole is reconstructed from the solution of the line-search for the duals.

The stuff here might be a relevant read (and the linked paper there)

EDIT: Trivial example in YALMIP (disclaimer, developed by me)

You don't want YALMIP to interpret this as a problem defined by the individual variables in $X$ and fit data into a dual description, which is default (leads to 125250 variables in dual /125250 equalities in primal)

X = sdpvar(500);
optimize([X>=0,trace(X)==1],trace(randn(500)*X),sdpsettings('solver','scs'))

Instead, you want YALMIP to interpret it from a primal point of view (1 equality in primal / 1 variable in dual)

optimize([X>=0,trace(X)==1],trace(randn(500)*X),sdpsettings('solver','scs','dualize',1))