11.10 Semidefinite Relaxation of MIQCQO Problems

In this case study we will discuss a fairly common application for Semidefinite Optimization: to define a continuous semidefinite relaxation of a mixed-integer quadratic optimization problem. This section is based on the method by Park and Boyd [PB15].

We will focus on problems of the form:

(11.40)\[\begin{split}\begin{array}{ll} \minimize & x^T P x + 2q^T x \\ \st &x\in \integral^n \end{array}\end{split}\]

where \(q\in \real^n\) and \(P\in \PSD^{n\times n}\) is positive semidefinite. There are many important problems that can be reformulated as (11.40), for example:

• integer least squares: minimize \(\|Ax -b\|^2_2\) subject to \(x\in \integral^n\),

• closest vector problem: minimize \(\|v - z\|_2\) subject to \(z\in \{ Bx~|~x\in \integral^n\}\).

Following [PB15], we can derive a relaxed continuous model. We first relax the integrality constraint

\[\begin{split}\begin{array}{lll} \minimize & x^T P x + 2q^T x &\\ \st &x_i(x_i-1) \geq 0 & i=1,\dots,n. \end{array}\end{split}\]

The last constraint is still non-convex. We introduce a new variable \(X\in \real^{n\times n}\), such that \(X = x\cdot x^T\). This allows us to write an equivalent formulation:

\[\begin{split}\begin{array}{ll} \minimize & \trace(PX) + 2q^T x \\ \st & \diag(X) \geq x, \\ & X = x\cdot x^T. \end{array}\end{split}\]

To get a conic problem we relax the last constraint and apply the Schur complement. The final relaxation follows:

(11.41)\[\begin{split}\begin{array}{ll} \minimize & \trace(PX) + 2q^T x \\ \st & \diag(X) \geq x, \\ & \left[ \begin{array} {cc} X & x \\ x^T & 1 \end{array}\right] \in \PSD^{n+1}. \end{array}\end{split}\]

Fusion Implementation

Implementing model (11.41) in Fusion is very simple. We assume the input \(n\), \(P\) and \(q\). Then we proceed creating the optimization model

    M = Model()

The important step is to define a single PSD variable

\[\begin{split}Z = \left[ \begin{array} {cc} X & x \\ x^T & 1 \end{array}\right] \in \PSD^{n+1}.\end{split}\]

Our code will create \(Z\) and two slices that correspond to \(X\) and \(x\):

    Z = M.variable("Z", Domain.inPSDCone(n+1))
    X = Z[:-1, :-1]
    x = Z[:-1, -1:]

Then we define the constraints:

    M.constraint( X.diag() >= x )
    M.constraint( Z[n, n] == 1 )

The objective function uses several available linear expressions:

    M.objective( ObjectiveSense.Minimize, 
        Expr.dot( P, X ) + 2.0 * Expr.dot(x, q) )

Note that the trace operator is not directly available in Fusion, but it can easily be defined from scratch.

Complete code

Listing 11.22 Fusion implementation of model (11.41). Click here to download.

def miqcqp_sdo_relaxation(n,P,q):
    M = Model()

    M.setLogHandler(sys.stdout) 

    Z = M.variable("Z", Domain.inPSDCone(n+1))
    X = Z[:-1, :-1]
    x = Z[:-1, -1:]

    M.constraint( X.diag() >= x )
    M.constraint( Z[n, n] == 1 )

    M.objective( ObjectiveSense.Minimize, 
        Expr.dot( P, X ) + 2.0 * Expr.dot(x, q) )

    return M

Numerical Examples

We present now some simple numerical experiments for the integer least squares problem:

(11.42)\[\begin{split}\begin{array}{ll} \minimize & \|Ax-b\|_2^2 \\ \st &x\in \integral^n. \end{array}\end{split}\]

It corresponds to the problem (11.40) with \(P=A^TA\) and \(q=-A^Tb\). Following [PB15] we will generate the input data by taking all entries of \(A\) from the normal distribution \(\mathcal{N}(0,1)\) and setting \(b=Ac\) where \(c\) comes from the uniform distribution on \([0,1]\).

An integer rounding xRound of the solution to (11.41) is a feasible integer solution to (11.42). We can compare it to the actual optimal integer solution xOpt, whenever the latter is available. Of course it is very simple to formulate the integer least squares problem in Fusion:

def int_least_squares(n, A, b):
    M = Model()

    M.setLogHandler(sys.stdout) 

    x = M.variable("x", n, Domain.integral(Domain.unbounded()))
    t = M.variable("t", 1, Domain.unbounded())

    M.constraint( Expr.vstack(t, A @ x - b), Domain.inQCone() )
    M.objective( ObjectiveSense.Minimize, t )

    return M

All that remains is to compare the values of the objective function \(\|Ax-b\|_2\) for the two solutions.

Listing 11.23 The comparison of two solutions. Click here to download.

# problem dimensions
n = 20
m = 2*n

# problem data
A = numpy.reshape(numpy.random.normal(0., 1.0, n*m), (m,n))
c = numpy.random.uniform(0., 1.0, n)
P = A.transpose().dot(A)
q = - P.dot(c)
b = A.dot(c)

# solve the problems
M = miqcqp_sdo_relaxation(n, P, q)
Mint = int_least_squares(n, A, b)
M.solve()
Mint.solve()

M.writeTask('M.ptf')
Mint.writeTask('Mint.ptf')

# rounded and optimal solution
xRound = numpy.rint(M.getVariable("Z").slice([0,n], [n,n+1]).level())
xOpt = numpy.rint(Mint.getVariable("x").level())

print(M.getSolverDoubleInfo("optimizerTime"), Mint.getSolverDoubleInfo("optimizerTime"))
print(numpy.linalg.norm(A.dot(xRound)-b), numpy.linalg.norm(A.dot(xOpt)-b)) 

Experimentally the objective value for xRound approximates the optimal solution with a factor of \(1.1\)-\(1.4\). We refer to [PB15] for a more involved iterative rounding procedure, producing integer solutions of even better quality, and for a detailed discussion of test results.