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:
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
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:
To get a conic problem we relax the last constraint and apply the Schur complement. The final relaxation follows:
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
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
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:
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]\).
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 \(\|Ax-b\|_2\) for the two solutions.
// 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.
