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{array}{r}\begin{array}{ll}\text{minimize}& x^{T}Px+2q^{T}x\\ \text{subject to}& x\in {\mathbb{Z}}^n\end{array}\end{array}$$

where $q\in {\mathbb{R}}^n$ and $P\in {\mathcal{S}}_{+}^{n\times n}$ is positive semidefinite. There are many important problems that can be reformulated as (11.40), for example:

• integer least squares: minimize $\Vert Ax-b{\Vert }_2^2$ subject to $x\in {\mathbb{Z}}^n$,

• closest vector problem: minimize $\Vert v-z{\Vert }_2$ subject to $z\in \{Bx|x\in {\mathbb{Z}}^n\}$.

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

$$\begin{array}{r}\begin{array}{lll}\text{minimize}& x^{T}Px+2q^{T}x& \\ \text{subject to}& x_i(x_i-1)\ge 0& i=1,\dots ,n.\end{array}\end{array}$$

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

$$\begin{array}{r}\begin{array}{ll}\text{minimize}& \mathrm{Tr}(PX)+2q^{T}x\\ \text{subject to}& \mathbf{diag}(X)\ge x,\\ & X=x\cdot x^T.\end{array}\end{array}$$

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

(11.41)$$\begin{array}{r}\begin{array}{ll}\text{minimize}& \mathrm{Tr}(PX)+2q^{T}x\\ \text{subject to}& \mathbf{diag}(X)\ge x,\\ & \left[\begin{array}{cc}X& x\\ x^T& 1\end{array}\right]\in {\mathcal{S}}_{+}^{n+1}.\end{array}\end{array}$$

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

  Model::t M = new Model();

The important step is to define a single PSD variable

$$\begin{array}{r}Z=\left[\begin{array}{cc}X& x\\ x^T& 1\end{array}\right]\in {\mathcal{S}}_{+}^{n+1}.\end{array}$$

Our code will create $Z$ and two slices that correspond to $X$ and $x$:

  Variable::t Z = M->variable("Z", Domain::inPSDCone(n + 1));

  Variable::t X = Z->slice(new_array_ptr<int, 1>({0, 0}), new_array_ptr<int, 1>({n, n}));
  Variable::t x = Z->slice(new_array_ptr<int, 1>({0, n}), new_array_ptr<int, 1>({n, n + 1}));

Then we define the constraints:

  M->constraint( Expr::sub(X->diag(), x), Domain::greaterThan(0.) );
  M->constraint( Z->index(n, n), Domain::equalsTo(1.) );

The objective function uses several available linear expressions:

  M->objective( ObjectiveSense::Minimize, Expr::add(
                  Expr::sum( Expr::mulElm( P, X ) ),
                  Expr::mul( 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.

Model::t miqcqp_sdo_relaxation(int n, Matrix::t P, const std::shared_ptr<ndarray<double, 1>> & q) {
  Model::t M = new Model();

  Variable::t Z = M->variable("Z", Domain::inPSDCone(n + 1));

  Variable::t X = Z->slice(new_array_ptr<int, 1>({0, 0}), new_array_ptr<int, 1>({n, n}));
  Variable::t x = Z->slice(new_array_ptr<int, 1>({0, n}), new_array_ptr<int, 1>({n, n + 1}));

  M->constraint( Expr::sub(X->diag(), x), Domain::greaterThan(0.) );
  M->constraint( Z->index(n, n), Domain::equalsTo(1.) );

  M->objective( ObjectiveSense::Minimize, Expr::add(
                  Expr::sum( Expr::mulElm( P, X ) ),
                  Expr::mul( 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{array}{r}\begin{array}{ll}\text{minimize}& \Vert Ax-b{\Vert }_2^2\\ \text{subject to}& x\in {\mathbb{Z}}^n.\end{array}\end{array}$$

It corresponds to the problem (11.40) with $P=A^{T}A$ and $q=-A^{T}b$. 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]$.

We implement the linear algebra operations using the LinAlg module available in MOSEK.

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:

Model::t int_least_squares(int n, Matrix::t A, const std::shared_ptr<ndarray<double, 1>> & b) {
  Model::t M = new Model();

  Variable::t x = M->variable("x", n, Domain::integral(Domain::unbounded()));
  Variable::t t = M->variable("t", 1, Domain::unbounded());

  M->constraint( Expr::vstack(t, Expr::sub(Expr::mul(A, x), b)), Domain::inQCone() );
  M->objective( ObjectiveSense::Minimize, t );

  return M;
}

All that remains is to compare the values of the objective function $\Vert Ax-b{\Vert }_2$ for the two solutions.

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

  // problem dimensions
  int n = 20;
  int m = 2 * n;

  auto c = new_array_ptr<double, 1>(n);
  auto A = new_array_ptr<double, 1>(n * m);
  auto P = new_array_ptr<double, 1>(n * n);
  auto b = new_array_ptr<double, 1>(m);
  auto q = new_array_ptr<double, 1>(n);


  std::generate(A->begin(), A->end(), std::bind(normal_distr, generator));
  std::generate(c->begin(), c->end(), std::bind(unif_distr, generator));
  std::fill(b->begin(), b->end(), 0.0);
  std::fill(q->begin(), q->end(), 0.0);

  // P = A^T A
  syrk(MSK_UPLO_LO, MSK_TRANSPOSE_YES,
       n, m, 1.0, A, 0., P);
  for (int j = 0; j < n; j++) for (int i = j + 1; i < n; i++) (*P)[i * n + j] = (*P)[j * n + i];

  // q = -P c, b = A c
  gemv(MSK_TRANSPOSE_NO, n, n, -1.0, P, c, 0., q);
  gemv(MSK_TRANSPOSE_NO, m, n, 1.0, A, c, 0., b);

  // Solve the problems
  {
    Model::t M = miqcqp_sdo_relaxation(n, Matrix::dense(n, n, P), q);
    Model::t Mint = int_least_squares(n, Matrix::dense(n, m, A)->transpose(), b);
    M->solve();
    Mint->solve();

    auto xRound = M->getVariable("Z")->
                  slice(new_array_ptr<int, 1>({0, n}), new_array_ptr<int, 1>({n, n + 1}))->level();
    for (int i = 0; i < n; i++) (*xRound)[i] = round((*xRound)[i]);
    auto yRound = new_array_ptr<double, 1>(m);
    auto xOpt   = Mint->getVariable("x")->level();
    auto yOpt   = new_array_ptr<double, 1>(m);
    std::copy(b->begin(), b->end(), yRound->begin());
    std::copy(b->begin(), b->end(), yOpt->begin());
    gemv(MSK_TRANSPOSE_NO, m, n, 1.0, A, xRound, -1.0, yRound);        // Ax_round-b
    gemv(MSK_TRANSPOSE_NO, m, n, 1.0, A, xOpt, -1.0, yOpt);            // Ax_opt-b

    std::cout << M->getSolverDoubleInfo("optimizerTime") << " " << Mint->getSolverDoubleInfo("optimizerTime") << "\n";
    double valRound, valOpt;
    dot(m, yRound, yRound, valRound); dot(m, yOpt, yOpt, valOpt);
    std::cout << sqrt(valRound) << " " << sqrt(valOpt) << "\n";
  }

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.