Cvxpy: MOSEK interface: inefficient encoding for SDPs

[juesato]

Hi,
I'm unable to solve mid-size SDPs with CVXPY+Mosek (~4000x4000 SDP variable). I can solve the equivalent SDP with YALMIP+Mosek in Matlab. I'm wondering if there is a different way to encode the SDP that I'm missing, or some difference in how CVXPY and YALMIP model the problem?

On Cvxpy, I get a mosek.Error: rescode.err_space(1051): Out of space. error. In particular, for Cvxpy, I see 32020003 Constraints (from the Mosek verbose debug output), whereas in Matlab, there are 4001, so it seems the SDP constraint is handled very differently in the two.

Here is the script in Cvxpy:

import cvxpy as cp
import numpy as np

def cvx_toy():
  size = 4000
  P = cp.Variable(shape=(size + 1, size + 1))
  x = P[1:, 0]
  X = P[1:, 1:]
  constraints = [
      P == P.T,
      P[0][0] == 1.0,
      P >> 0,
  ]
  diag = cp.atoms.affine.diag.diag
  # b1, b2 = np.random.normal(size=size), np.random.normal(size=size)
  b1 = np.array(range(size)) + 1
  b2 = np.array(range(size)) + 3
  lb, ub = np.minimum(b1, b2), np.maximum(b1, b2)
  constraints += [diag(X) <= (cp.multiply(lb + ub, x) -
                              cp.multiply(lb, ub))]
  obj_cp = cp.sum(x)
  problem = cp.Problem(cp.Minimize(obj_cp), constraints)
  print('Start')
  # problem.solve(solver=cp.SCS)
  problem.solve(solver=cp.MOSEK, verbose=True)
  print(problem.status) 
  print('Opt:', sum(lb))
  print('SDP:', obj_cp.value)
cvx_toy()

And the equivalent script in YALMIP (Matlab):

function [] = cvxtoy()

addpath(genpath('/home/juesato/code/matlab_packages/yalmip/YALMIP-master'))
addpath(genpath('/home/juesato/install/mosek/9.2/toolbox/r2015a'))

sz = 4000
P = sdpvar(sz+1, sz+1)
x = P(2:sz+1, 1)
X = P(2:sz+1, 2:sz+1)
constraints = [P(1, 1) == 1, P >= 0]
b1 = (1:sz)'
b2 = b1 + 2
lb = min(b1, b2)
ub = max(b1, b2)
constraints = [constraints, diag(X) <= (lb+ub) .* x - lb.*ub]
obj_cp = sum(x)
diagnostics = optimize(constraints, obj_cp, sdpsettings('dualize', 1, 'solver', 'mosek'))
val = value(obj_cp)

[aszekMosek]

It is sort of even worse. Even if I use the syntax which avoids adding explicit P==P.T constraints and write just the simplest model

import cvxpy as cp

P = cp.Variable((100,100), PSD=True)
problem = cp.Problem(cp.Minimize(0), [P[0][0]==1])
problem.solve(solver=cp.MOSEK, verbose=True, save_file="prob.ptf")

then cvxpy still creates n*(n+1)/2 auxiliary scalar variables x and n^2 equality constraints setting the x equal to the P entries (half of which are anyway redundant). Unfortunately that's what kills all of MOSEK semidefinite performance in this case.

The only quick workaround I see for you is to use directly Python Fusion https://docs.mosek.com/9.2/pythonfusion/modeling.html, where this problem does not appear and the code constructing the model will be very similar.

[rileyjmurray]

EDIT: actually hold off on the points below, while I test the first claim.

DOUBLE EDIT: the remarks below have some truth to them, but they aren't be the end of the story. See my next comment.

@juesato this inefficiency should be isolated to the MOSEK interface, because of the way I originally wrote it a couple years ago. I know that MOSEK is CVXPY's only competitive large-scale SDP solver, but you should find that the SCS and CVXOPT interfaces compile into a reasonable size. As @aszekMosek said, the short-term solution here is probably to use MOSEK Fusion.

@aszekMosek can you look at possibly removing the SDP inefficiency in our MOSEK interface? The first line where things get dicey is in MOSEK(ConicSolver).apply, and this has a downstream effect on MOSEK(ConicSolver).solve_via_data. I don't think fixing this will affect SDP dual variable recovery.

[rileyjmurray]

@juesato as @aszekMosek said, avoiding P.T == P can reduce the size of the problem. However our MOSEK interface does compile into standard-form in less efficient ways than our SCS and CVXOPT intefaces.

I ran your problem with three different formulations, and compared with the dimensions of matrices generated for these three solvers. I'll paste the script below, provide its output, and then comment on its output.

import cvxpy as cp
import numpy as np


def prob_v1_setup(size=100):
    P = cp.Variable(shape=(size + 1, size + 1))
    constraints = [P == P.T, P[0][0] == 1.0, P >> 0]
    return P, constraints

def prob_v2_setup(size=100):
    P = cp.Variable(shape=(size + 1, size + 1), symmetric=True)
    constraints = [P[0][0] == 1.0, P >> 0]
    return P, constraints

def prob_v3_setup(size=100):
    P = cp.Variable(shape=(size + 1, size + 1), PSD=True)
    constraints = [P[0][0] == 1.0]
    return P, constraints

def common_setup(P, constraints):
    x = P[1:, 0]
    X = P[1:, 1:]
    size = P.shape[0]-1
    diag = cp.atoms.affine.diag.diag
    b1 = np.array(range(size)) + 1
    b2 = np.array(range(size)) + 3
    lb, ub = np.minimum(b1, b2), np.maximum(b1, b2)
    constraints += [diag(X) <= (cp.multiply(lb + ub, x) -
                              cp.multiply(lb, ub))]
    obj_cp = cp.sum(x)

    problem = cp.Problem(cp.Minimize(obj_cp), constraints)
    return problem

def prob_v1(size=100):
    P, constraints = prob_v1_setup(size)
    prob = common_setup(P, constraints)
    return prob


def prob_v2(size=100):
    P, constraints = prob_v2_setup(size)
    prob = common_setup(P, constraints)
    return prob


def prob_v3(size=100):
    P, constraints = prob_v2_setup(size)
    prob = common_setup(P, constraints)
    return prob

def solver_datas(problem):
    scs_data = problem.get_problem_data(solver=cp.SCS)
    cvxopt_data = problem.get_problem_data(solver=cp.CVXOPT)
    mosek_data = problem.get_problem_data(solver=cp.MOSEK)
    return scs_data, cvxopt_data, mosek_data

def print_matrix_shapes(cvxopt_data, scs_data, mosek_data):
    print('CVXOPT')
    print('\t' + str((cvxopt_data[0]['A'].shape, cvxopt_data[0]['G'].shape)))
    print('SCS')
    print('\t' + str(scs_data[0]['A'].shape))
    print('MOSEK')
    print('\t' + str(mosek_data[0]['G'].shape))
    print()


print('Version 1')
scs_data, cvxopt_data, mosek_data = solver_datas(prob_v1())
print_matrix_shapes(cvxopt_data, scs_data, mosek_data)

print('Version 2')
scs_data, cvxopt_data, mosek_data = solver_datas(prob_v2())
print_matrix_shapes(cvxopt_data, scs_data, mosek_data)

print('Version 3')
scs_data, cvxopt_data, mosek_data = solver_datas(prob_v3())
print_matrix_shapes(cvxopt_data, scs_data, mosek_data)

Here's the output. The tuples printed by CVXOPT give the shapes of an equality constraint matrix A @ x == b and a conic constraint matrix G @ x <=_K h. The tuples for SCS and MOSEK are conic constraint matrices which include equality constraints.

Version 1
CVXOPT
    ((10202, 10201), (1, 10201))
SCS
    (15453, 10201)
MOSEK
    (20503, 10201)
Version 2
CVXOPT
    ((1, 5151), (10301, 5151))
SCS
    (5252, 5151)
MOSEK
    (10302, 5151)
Version 3
CVXOPT
    ((1, 5151), (10301, 5151))
SCS
    (5252, 5151)
MOSEK
    (10302, 5151)

In Versions 2 and 3, every solver declares 5151 scalar variables, which the minimum number of scalar variables necessary to represent a PSD matrix of order 101. In pretty much every version we see that the current MOSEK interface generates twice as many conic constraints as necessary. This is because of the implementation inefficiency of MOSEK(ConicSolver) I mentioned in my earlier comment, and it's something that can probably be fixed with a reasonable amount of effort.

A larger problem with the MOSEK interface, is that the 5151 scalar variables should never be declared in the first place. The best approach is to declare a PSD matrix variable (of size 101) directly with the MOSEK Task. We don't do this right now, because the way that CVXPY compiles problems into standard-form automatically vectorizes all variables. It's plausible that the MOSEK interface could bypass this problem by looking at a vectorized feasible set A @ x + b in K, and then identifying blocks in (A,b,K) which specify that a vectorized variable belongs to the PSD cone. If those blocks of constraints were identified, the MOSEK interface could "de-vectorize" the PSD variables. This might take a significant amount of effort to fix.

[SteveDiamond]

This has been an issue since we first added the MOSEK interface. The MOSEK standard form is quite different from the standard forms CVXPY was originally designed for (SCS and CVXOPT). But I'm sure someone can figure out the correct transformation!

[aszekMosek]

@rileyjmurray @SteveDiamond I will try to take a look at it. It will take some time, also because it was a while since I last understood cvxpy internals. Maybe it can be hacked just inside the Mosek interface.

[rileyjmurray]

@aszekMosek we're in the process of merging in changes to the ConicSolver base class, and we're adding lots of web-documentation for solver interfaces (#1020). Within the next couple weeks it should be easier to devise fixes for the inefficiencies here.

[angeris]

I'm just bumping this thread since a similar case has come up again w/ the Mosek interface for some numerical experiments I'm running on medium-sized SDPs.

Big thanks for working on this @rileyjmurray ! :)

[rileyjmurray]

@aszekMosek @angeris @juesato this is just to let you know that I'm working on this. My next PR will be a thorough overhaul of the MOSEK interace. The main issue is how to efficiently convert from CVXPY's representation of a feasible set into a form exected by MOSEK.

• For mixed-integer problems, the conversion will use slack variables (like the current interface does). Note that MOSEK doesn't support mixed-integer SDP so this situation never happens with SDP.
• For continuous problems, dualize CVXPY's standard form min{ c.T @ x : A @ x + b in K} into MOSEK's standard form max{-b.T @ y : A.T @ y == c, y in K^*}. Once in MOSEK's standard form, adding SDP variables should also be possible without slacks.

As of this moment, I've made changes so continuous non-SDP problems work correctly (all tests passing). The next step is mixed-integer problems, and then finally SDP.

I'll update this thread periodically as I make further changes.

[rileyjmurray]

Mixed-integer now works as well. All that's left is SDP.

[angeris]

Great, thank you @rileyjmurray ! Looking forward to it :)

[rileyjmurray]

@aszekMosek @angeris SDP is taken care of. The MOSEK interface is now completely rewritten. I'll make a PR soon.

[rileyjmurray]

Documentation for the "Dualize" reduction is done. I need to write documentation for the "Slacks" reduction, then I can make the PR.