# Efficiently solving SDP relaxation of an integer quadratic program

I have an integer quadratic program of the form,

\begin{align} \underset{x}{\max}&\;\;\|Ax-b\|_2^2\\ \text{subject to}&\;\;x\in{\bf Z}\geq0 \end{align}

I'm currently using the (admittedly crude) convex approximation which relaxes the integer constraint, and then recovers a feasible point by rounding each component of $$x^*$$ to the nearest integer.

I would like to instead use the more sophisticated semidefinite relaxation presented in equation (11) of (Park & Boyd, 2017),

\begin{align} \underset{X,\hspace{1mm}x}{\max}&\;\;\text{tr}(PX)+2q^Tx\\ \text{subject to}&\;\;\text{diag}\hspace{0.5mm}X\geq x\\[3pt] &\;\;[X\;x;\;x^T\;1] \succcurlyeq 0\\[3pt] &\;\;x\geq0 \end{align}

where $$P=A^TA$$ and $$q=-A^Tb$$.

I've tried several different free solvers, and the SCS solver seems to work the best, however if we suppose that $$P$$ is an $$n\times n$$ matrix, then on a single machine I can only scale the problem up to about $$n=100$$. Ideally I would like to scale the problem up to at least $$n=1000$$, and beyond if possible.

Is this possible? For what it's worth, $$P$$ is approximately a band matrix, and for $$n\geq1000$$, the number of nonzero entries of $$P$$ is less than 1%, and it continues to get sparser as $$n$$ gets larger.

I have access to servers with many CPU cores and GPUs if necessary.

• $n=100$ (or even $n=1000$) seems quite small and well within the range of problems that can routinely be solved to high accuracy using primal-dual interior-point methods. May 10 at 15:10
• @BrianBorchers $n=1000$ is over half a million variables, you're saying that that should be considered a small problem size for a semidefinite program? May 10 at 15:36
• Yes, we regularly solve problems with matrix variables that are on the order of n=10000 or more. The number of constraints is also important, but you've only got 2n inequality constraints. If you'll share an example problem with me, I'd be happy to solve it. May 10 at 16:25
• You have to be very careful in how you interpret and model a large SDP of this type. When you say $n=1000$ means half a million variables, it means you interpret it as a parameterization of the standard dual model of an SDP. However, as Brian hints at, this model should be interpreted from the primal side, and in that view it is of size O(n) yalmip.github.io/tutorial/automaticdualization. May 10 at 16:56

As a follow up to my comment, where I say that the model has to be interpreted from the primal model side, here it implements the model for $$n=1000$$ and then tells the modelling layer YALMIP that it should extract a model interpreting this a representation of the primal. The SDP solver Mosek solves this in 1 minute on an old laptop

n = 1000;
P = sprandn(n,n,1/n);
P = P*P'+eye(n);
q = randn(n,1);
Z = sdpvar(n+1);
X = Z(1:n,1:n);
x = Z(1:n,end);
Model =[Z >= 0,
Z(end,end)==1,
diag(X)>=x,
x>=0];
Objective = P(:)'*X(:)+2*q'*x;
optimize(Model,Objective,sdpsettings('dualize',1))

I did not wait for SCS to terminate, but for $$n=100$$ SCS solved it in 2 seconds.

• I don't have Mosek (or Matlab), however I translated your example into Julia and ran it using SCS for $n=1000$, and it took 90 minutes. So I guess the takeaway is to use a better solver, or if you don't have access to one, then just write the algorithm yourself. May 12 at 4:24
• SCS should not take 90 minutes. As I said, you have to interpret it from the primal side, meaning you either have to use some dualization functionality to tell your modelling language to do so (it now sounds like it interprets it from the dual side i.e. parameterization of the matrix $Z$ with variables instead of viewing it as the cone), or manually derive the dual of the model and input that to the modelling layer. I recommend you to read the paper referenced on the Wiki page on automatic dualization linked above in the comment. May 12 at 6:07
• SCS (both direct and indirect) solves the primal form model in roughly 5 minutes on my machine. May 12 at 6:16