[mlpack] sdp primal dual solver is very slow
Joe Dinius
josephwdinius at gmail.com
Mon Sep 16 12:23:47 EDT 2019
I have attached code (see bottom) utilizing ensmallen's primal dual sdp
solver to solve a moderately sized SDP problem (257x257) with 67
constraints. My problem is this: The code is painfully slow (anywhere
from 48-60 seconds per iteration), even when compiling with -O3 and openmp
enabled. As a comparison, there is a matlab package called SDPNAL+(
http://www.math.nus.edu.sg/~mattohkc/SDPNALplus.html) that was able to
solve the same-sized problem as the one below in 5 iterations in under 10
seconds. SDPNAL+ is implemented in matlab, and I have a very hard time
believing that there is such a performance discrepancy, even discounting
the difference in underlying algorithm, between these two implementations.
My intuition is that I am doing something very, very foolish and I am
hoping someone here can shed some light on what it may be. Here is the
code:
//! c/c++ headers
#include <algorithm>
#include <cmath>
#include <limits>
#include <map>
#include <string>
#include <utility>
#include <vector>
//! dependency headers
#include <ensmallen.hpp>
//! project headers
/**
* @brief SdpConstraintMap
*
* Setup sequence of matrices {Ai} and sequences of scalars {Bi} s.t. <Ai,
Y> == Bi for all i,
* where 0 \le i \lt (number of constraints).
*
* @note Sequential ordering is:
* 1) inequality constraints from genSDPConstraintMap (modeled as
equality constraints with slack variables)
* 2) equality constraints from genSDPConstraintMap
* 3) Y(N, N) == 1
* 4) Tr(Y) == k+1
* 5) slack variable constraints (all 1's in row/column pivot except
where row_id == col_id, at which the value is 0)
*/
class SdpConstraintMap {
private:
size_t num_constraints_;
std::vector<arma::sp_mat> A_;
arma::vec B_;
public:
/**
* @brief Constructor
*
* @param [in] N (square of) number of subsampling points for sdp
problem definition
* @param [in] k (square of) number of point matches to pick within
subsampled subset
* @return constructed SdpConstraintMap object
*/
SdpConstraintMap(size_t const & N, size_t const & k) {
//! get size of W/Y matrices from value-copy in scope
size_t const n = std::sqrt(N);
//! there are 2*n inequality constraints, so 2*n slack variables
are needed
size_t const n_slack = 2*n;
//! total size of Y = N + n_slack + 1
size_t const size_Y = N + n_slack + 1;
//! initialize containers
num_constraints_ = 4*n+3;
A_.resize(num_constraints_);
B_.resize(num_constraints_);
B_.zeros();
//! CONSTRAINTS
size_t i_constr = 0;
//! Ai*Y <= bi
for (size_t i = 0; i < n; ++i) {
arma::sp_mat A(size_Y, size_Y);
for (size_t j = i; j < N ; j+=n) {
A(N, j) = 1.0;
}
//! \sum_i Y(N, :) + si = 1, si >= 0
A(N+i+1, N+i+1) = 1.0;
B_(i_constr) = 1.0;
A_[i_constr++] = A;
}
for (size_t i = 0; i < n; ++i) {
arma::sp_mat A(size_Y, size_Y);
for (size_t j = 0; j < n ; ++j) {
A(N, i*n+j) = 1.0;
}
//! \sum_i Y(N, :) + si = 1, si >= 0
A(N+i+n+1, N+i+n+1) = 1.0;
B_(i_constr) = 1.0;
A_[i_constr++] = A;
}
//! Aeq*Y == beq
arma::mat Aeq(size_Y, size_Y, arma::fill::zeros);
Aeq(N, arma::span(0, N-1)).fill(1.0);
B_(i_constr) = static_cast<double>(k);
A_[i_constr++] = arma::sp_mat(Aeq);
//! Y(N+1, N+1) == 1
arma::sp_mat A0(size_Y, size_Y);
A0(N, N) = 1.0;
B_(i_constr) = 1.0;
A_[i_constr++] = A0;
//! Trace(Y) == (< I, Y >) == k + 1
arma::mat A1(size_Y, size_Y, arma::fill::eye);
A1(arma::span(N+1, size_Y-1), arma::span(N+1, size_Y-1)).fill(0.0);
B_(i_constr) = static_cast<double>(k+1);
A_[i_constr++] = arma::sp_mat(A1);
/**
* Constraint matrices for the slack variables look like the
following:
*
* |0 0 0 ... 1 0 0 0|
* |0 0 0 ... 1 0 0 0|
* |0 0 0 ... 1 0 0 0|
* | . |
* | . |
* | . |
* |1 1 1 ... 0 1 1 1|
* |0 0 0 ... 1 0 0 0|
* |0 0 0 ... 1 0 0 0|
* | . |
* | . |
* |0 0 0 ... 1 0 0 0|
*
* which ensures positive semidefinite values of the slack
variables
*/
for (size_t i = 0; i < n_slack; ++i) {
arma::mat Aslack(size_Y, size_Y, arma::fill::zeros);
Aslack.row(N+1+i).fill(1.0);
Aslack.col(N+1+i).fill(1.0);
Aslack(N+1+i, N+1+i) = 0.0;
B_(i_constr) = 0.0;
A_[i_constr++] = arma::sp_mat(Aslack);
}
}
/**
* @brief Get number of constraints
*
* @note const method
*
* @return num_constraints_ (copy of) number of constraints
*/
size_t const getNumConstraints() const noexcept { return
num_constraints_; }
/**
* @brief Get sequence of lhs constraint matrices
*
* @note const method
*
* @return A_ (copy of) sequence of lhs constraint matrices
*/
std::vector< arma::sp_mat > const getConstraintMapA() const noexcept {
return A_; }
/**
* @brief Get sequence of rhs constraints
*
* @note const method
*
* @return B_ (copy of) sequence of rhs constraints
*/
arma::vec const getConstraintMapB() const noexcept { return B_; }
};
/**
* @brief Solve SDP:
* max <W, Y> subject to constraints <Ai, Y> == bi
*
* @note ensmallen doesn't support type templating for containers, so just
use doubles
*
* @param [in] W weight matrix
* @param [in] constr_map sdp constraint mapping
* @param [in][out] Y optimal solution
* @return true if converged, false otherwise
*/
bool sdp(arma::mat const & W, SdpConstraintMap const & constr_map,
arma::mat & Y) {
size_t const size_Y = constr_map.getConstraintMapA()[0].n_rows;
ens::SDP<arma::sp_mat> sdp_obj(size_Y, constr_map.getNumConstraints(),
0);
//! objective: maximize <C, Y>
sdp_obj.C()(arma::span(0, W.n_rows-1), arma::span(0, W.n_cols-1)) = W;
//! subject to <Ai, Y> == bi for 0 \le i \lt num_constraints
sdp_obj.SparseA() = constr_map.getConstraintMapA();
sdp_obj.SparseB() = constr_map.getConstraintMapB();
//! setup solver
arma::mat _Y(size_Y, size_Y);
ens::PrimalDualSolver<ens::SDP<arma::sp_mat>> solver(sdp_obj);
solver.MaxIterations() = 10; // XXX leave at `1` until I figure out
why the solver is so slow
solver.Tau() = 0.99;
solver.PrimalInfeasTol() = 0.1;
solver.DualInfeasTol() = 0.1;
static_cast<void>(solver.Optimize(_Y));
//! remove slack variables
Y = _Y(arma::span(0, W.n_rows), arma::span(0, W.n_cols));
return true; // XXX figure out criteria for quality of solution before
returning true
}
int main(int argc, char* argv[]) {
auto const constraints = SdpConstraintMap(256, 4);
arma::Mat<T> const W(257, 257, arma::fill::eye);
arma::Mat<T> solution(257, 257);
static_cast<void>(sdp(W, constraints, solution));
return 0;
}
Thanks!
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://knife.lugatgt.org/pipermail/mlpack/attachments/20190916/a38fb24f/attachment.html>
More information about the mlpack
mailing list