[mlpack-git] master: fix up PrimalDualSolver::Optimize() return value (736f62c)
gitdub at big.cc.gt.atl.ga.us
gitdub at big.cc.gt.atl.ga.us
Thu Mar 5 22:14:11 EST 2015
Repository : https://github.com/mlpack/mlpack
On branch : master
Link : https://github.com/mlpack/mlpack/compare/904762495c039e345beba14c1142fd719b3bd50e...f94823c800ad6f7266995c700b1b630d5ffdcf40
>---------------------------------------------------------------
commit 736f62c32d76cdb31ff1b5c84f624fde2d2516cb
Author: Stephen Tu <stephent at berkeley.edu>
Date: Wed Jan 21 18:09:41 2015 -0800
fix up PrimalDualSolver::Optimize() return value
>---------------------------------------------------------------
736f62c32d76cdb31ff1b5c84f624fde2d2516cb
src/mlpack/core/optimizers/sdp/primal_dual.hpp | 75 +++++++++++++++++++---
.../core/optimizers/sdp/primal_dual_impl.hpp | 56 ++++++++--------
src/mlpack/tests/sdp_primal_dual_test.cpp | 70 ++++++++++++++++----
3 files changed, 151 insertions(+), 50 deletions(-)
diff --git a/src/mlpack/core/optimizers/sdp/primal_dual.hpp b/src/mlpack/core/optimizers/sdp/primal_dual.hpp
index 8b110d1..7af3793 100644
--- a/src/mlpack/core/optimizers/sdp/primal_dual.hpp
+++ b/src/mlpack/core/optimizers/sdp/primal_dual.hpp
@@ -12,59 +12,116 @@
namespace mlpack {
namespace optimization {
+/**
+ * Interface to a primal dual interior point solver.
+ *
+ * @tparam SDPType
+ */
template <typename SDPType>
class PrimalDualSolver {
public:
+ /**
+ * Construct a new solver instance from a given SDP instance.
+ * Uses a random, positive initialization point.
+ *
+ * @param sdp
+ */
PrimalDualSolver(const SDPType& sdp);
+ /**
+ * Construct a new solver instance from a given SDP instance. Uses a random,
+ * positive initialization point. Both initialX and initialZ need to be
+ * positive definite matrices.
+ *
+ * @param sdp
+ * @param initialX
+ * @param initialYsparse
+ * @param initialYdense
+ * @param initialZ
+ */
PrimalDualSolver(const SDPType& sdp,
const arma::mat& initialX,
const arma::vec& initialYsparse,
const arma::vec& initialYdense,
const arma::mat& initialZ);
- std::pair<bool, double>
- Optimize(arma::mat& X,
- arma::vec& ysparse,
- arma::vec& ydense,
- arma::mat &Z);
-
- std::pair<bool, double>
- Optimize(arma::mat& X)
+ /**
+ * Invoke the optimization procedure, returning the converged values
+ * for the primal and dual variables.
+ *
+ * @param X
+ * @param ysparse
+ * @param ydense
+ * @param Z
+ */
+ double Optimize(arma::mat& X,
+ arma::vec& ysparse,
+ arma::vec& ydense,
+ arma::mat &Z);
+
+ /**
+ * Invoke the optimization procedure, and only return the primal variable.
+ *
+ * @param X
+ */
+ double Optimize(arma::mat& X)
{
arma::vec ysparse, ydense;
arma::mat Z;
return Optimize(X, ysparse, ydense, Z);
}
+ //! Return the underlying SDP instance.
const SDPType& Sdp() const { return sdp; }
+ //! Modify tau. Typical values are 0.99.
double& Tau() { return tau; }
+ //! Modify the XZ tolerance.
double& NormXzTol() { return normXzTol; }
+ //! Modify the primal infeasibility tolerance.
double& PrimalInfeasTol() { return primalInfeasTol; }
+ //! Modify the dual infeasibility tolerance.
double& DualInfeasTol() { return dualInfeasTol; }
+ //! Modify the maximum number of iterations to run before converging.
size_t& MaxIterations() { return maxIterations; }
private:
+
+ //! The SDP problem instance to optimize
SDPType sdp;
+ //! Starting point for X. Needs to be positive definite.
arma::mat initialX;
+
+ //! Starting lagrange multiplier for the sparse constraints.
arma::vec initialYsparse;
+
+ //! Starting lagrange multiplier for the sparse constraints.
arma::vec initialYdense;
+
+ //! Starting point for Z, the complementary slack variable. Needs to be
+ //positive definite.
arma::mat initialZ;
+ //! The step size modulating factor. Needs to be a scalar in (0, 1).
double tau;
+
+ //! The tolerance on the norm of XZ required before terminating.
double normXzTol;
+
+ //! The tolerance required on the primal constraints required before terminating.
double primalInfeasTol;
+
+ //! The tolerance required on the dual constraint required before terminating.
double dualInfeasTol;
+ //! Maximum number of iterations to run. Set to 0 for no limit.
size_t maxIterations;
-
};
} // namespace optimization
diff --git a/src/mlpack/core/optimizers/sdp/primal_dual_impl.hpp b/src/mlpack/core/optimizers/sdp/primal_dual_impl.hpp
index 2d927b0..98eca00 100644
--- a/src/mlpack/core/optimizers/sdp/primal_dual_impl.hpp
+++ b/src/mlpack/core/optimizers/sdp/primal_dual_impl.hpp
@@ -91,6 +91,9 @@ PrimalDualSolver<SDPType>::PrimalDualSolver(const SDPType& sdp,
<< std::endl;
}
+/**
+ * alphahat = sup{ alphahat : A + dA is psd }
+ */
static inline double
AlphaHat(const arma::mat& A, const arma::mat& dA)
{
@@ -103,6 +106,9 @@ AlphaHat(const arma::mat& A, const arma::mat& dA)
return 1. / alphahatinv;
}
+/**
+ * alpha = min(1, tau * alphahat(A, dA))
+ */
static inline double
Alpha(const arma::mat& A, const arma::mat& dA, double tau)
{
@@ -162,7 +168,7 @@ SolveKKTSystem(const arma::sp_mat& Asparse,
// TODO(stephentu): use a more efficient method (e.g. LU decomposition)
if (!arma::solve(dy, M, rhs))
- Log::Fatal << "PrimalDualSolver::Optimize(): Could not solve KKT system" << std::endl;
+ Log::Fatal << "PrimalDualSolver::SolveKKTSystem(): Could not solve KKT system" << std::endl;
if (Asparse.n_rows)
dysparse = dy(arma::span(0, Asparse.n_rows - 1));
@@ -186,7 +192,7 @@ template <typename eT> struct vectype<arma::SpMat<eT>> { typedef arma::SpCol<eT>
} // namespace private_
template <typename SDPType>
-std::pair<bool, double>
+double
PrimalDualSolver<SDPType>::Optimize(arma::mat& X,
arma::vec& ysparse,
arma::vec& ydense,
@@ -240,7 +246,7 @@ PrimalDualSolver<SDPType>::Optimize(arma::mat& X,
Einv_F_AdenseT.set_size(n2bar, sdp.NumDenseConstraints());
M.set_size(sdp.NumConstraints(), sdp.NumConstraints());
- double primal_obj = 0., alpha, beta;
+ double primalObj = 0., alpha, beta;
for (size_t iteration = 0; iteration != maxIterations; iteration++)
{
if (sdp.NumSparseConstraints())
@@ -334,21 +340,17 @@ PrimalDualSolver<SDPType>::Optimize(arma::mat& X,
Z += beta * dZ;
math::Svec(Z, sz);
- const double norm_XZ = arma::norm(X * Z, "fro");
-
- const double sparse_primal_infeas = arma::norm(sdp.SparseB() - Asparse * sx, 2);
- const double dense_primal_infeas = arma::norm(sdp.DenseB() - Adense * sx, 2);
- const double primal_infeas = sqrt(
- sparse_primal_infeas * sparse_primal_infeas +
- dense_primal_infeas * dense_primal_infeas);
+ const double normXZ = arma::norm(X * Z, "fro");
- primal_obj = arma::dot(sdp.C(), X);
+ const double sparsePrimalInfeas = arma::norm(sdp.SparseB() - Asparse * sx, 2);
+ const double densePrimalInfeas = arma::norm(sdp.DenseB() - Adense * sx, 2);
+ const double primalInfeas = sqrt(sparsePrimalInfeas * sparsePrimalInfeas +
+ densePrimalInfeas * densePrimalInfeas);
- const double dual_obj =
- arma::dot(sdp.SparseB(), ysparse) +
- arma::dot(sdp.DenseB(), ydense);
+ primalObj = arma::dot(sdp.C(), X);
- const double duality_gap = primal_obj - dual_obj;
+ const double dualObj = arma::dot(sdp.SparseB(), ysparse) + arma::dot(sdp.DenseB(), ydense);
+ const double dualityGap = primalObj - dualObj;
// TODO(stephentu): this dual check is quite expensive,
// maybe make it optional?
@@ -357,27 +359,25 @@ PrimalDualSolver<SDPType>::Optimize(arma::mat& X,
DualCheck += ysparse(i) * sdp.SparseA()[i];
for (size_t i = 0; i < sdp.NumDenseConstraints(); i++)
DualCheck += ydense(i) * sdp.DenseA()[i];
- const double dual_infeas = arma::norm(DualCheck, "fro");
+ const double dualInfeas = arma::norm(DualCheck, "fro");
Log::Debug
<< "iter=" << iteration + 1 << ", "
- << "primal=" << primal_obj << ", "
- << "dual=" << dual_obj << ", "
- << "gap=" << duality_gap << ", "
- << "||XZ||=" << norm_XZ << ", "
- << "primal_infeas=" << primal_infeas << ", "
- << "dual_infeas=" << dual_infeas << ", "
+ << "primal=" << primalObj << ", "
+ << "dual=" << dualObj << ", "
+ << "gap=" << dualityGap << ", "
+ << "||XZ||=" << normXZ << ", "
+ << "primalInfeas=" << primalInfeas << ", "
+ << "dualInfeas=" << dualInfeas << ", "
<< "mu=" << mu
<< std::endl;
- if (norm_XZ <= normXzTol &&
- primal_infeas <= primalInfeasTol &&
- dual_infeas <= dualInfeasTol)
- return std::make_pair(true, primal_obj);
+ if (normXZ <= normXzTol && primalInfeas <= primalInfeasTol && dualInfeas <= dualInfeasTol)
+ return primalObj;
}
- Log::Warn << "Did not converge!" << std::endl;
- return std::make_pair(false, primal_obj);
+ Log::Warn << "Did not converge after " << maxIterations << " iterations!" << std::endl;
+ return primalObj;
}
} // namespace optimization
diff --git a/src/mlpack/tests/sdp_primal_dual_test.cpp b/src/mlpack/tests/sdp_primal_dual_test.cpp
index 1af340c..8b0d77c 100644
--- a/src/mlpack/tests/sdp_primal_dual_test.cpp
+++ b/src/mlpack/tests/sdp_primal_dual_test.cpp
@@ -180,6 +180,50 @@ ConstructMaxCutSDPFromLaplacian(const std::string& laplacianFilename)
return sdp;
}
+static void CheckPositiveSemiDefinite(const arma::mat& X)
+{
+ const auto evals = arma::eig_sym(X);
+ BOOST_REQUIRE_GE(evals(0), 1e-20);
+}
+
+template <typename SDPType>
+static void CheckKKT(const SDPType& sdp,
+ const arma::mat& X,
+ const arma::vec& ysparse,
+ const arma::vec& ydense,
+ const arma::mat& Z)
+{
+ // require that the KKT optimality conditions for sdp are satisfied
+ // by the primal-dual pair (X, y, Z)
+
+ CheckPositiveSemiDefinite(X);
+ CheckPositiveSemiDefinite(Z);
+
+ const double normXz = arma::norm(X * Z, "fro");
+ BOOST_REQUIRE_SMALL(normXz, 1e-5);
+
+ for (size_t i = 0; i < sdp.NumSparseConstraints(); i++)
+ {
+ BOOST_REQUIRE_SMALL(
+ fabs(arma::dot(sdp.SparseA()[i], X) - sdp.SparseB()[i]),
+ 1e-5);
+ }
+
+ for (size_t i = 0; i < sdp.NumDenseConstraints(); i++)
+ {
+ BOOST_REQUIRE_SMALL(
+ fabs(arma::dot(sdp.DenseA()[i], X) - sdp.DenseB()[i]),
+ 1e-5);
+ }
+
+ arma::mat dualCheck = Z - sdp.C();
+ for (size_t i = 0; i < sdp.NumSparseConstraints(); i++)
+ dualCheck += ysparse(i) * sdp.SparseA()[i];
+ for (size_t i = 0; i < sdp.NumDenseConstraints(); i++)
+ dualCheck += ydense(i) * sdp.DenseA()[i];
+ const double dualInfeas = arma::norm(dualCheck, "fro");
+ BOOST_REQUIRE_SMALL(dualInfeas, 1e-5);
+}
BOOST_AUTO_TEST_SUITE(SdpPrimalDualTest);
@@ -198,8 +242,8 @@ static void SolveMaxCutFeasibleSDP(const SDP<arma::sp_mat>& sdp)
arma::mat X, Z;
arma::vec ysparse, ydense;
- const auto p = solver.Optimize(X, ysparse, ydense, Z);
- BOOST_REQUIRE(p.first);
+ solver.Optimize(X, ysparse, ydense, Z);
+ CheckKKT(sdp, X, ysparse, ydense, Z);
}
static void SolveMaxCutPositiveSDP(const SDP<arma::sp_mat>& sdp)
@@ -220,8 +264,8 @@ static void SolveMaxCutPositiveSDP(const SDP<arma::sp_mat>& sdp)
arma::mat X, Z;
arma::vec ysparse, ydense;
- const auto p = solver.Optimize(X, ysparse, ydense, Z);
- BOOST_REQUIRE(p.first);
+ solver.Optimize(X, ysparse, ydense, Z);
+ CheckKKT(sdp, X, ysparse, ydense, Z);
}
BOOST_AUTO_TEST_CASE(SmallMaxCutSdp)
@@ -247,8 +291,8 @@ BOOST_AUTO_TEST_CASE(SmallLovaszThetaSdp)
arma::mat X, Z;
arma::vec ysparse, ydense;
- const auto p = solver.Optimize(X, ysparse, ydense, Z);
- BOOST_REQUIRE(p.first);
+ solver.Optimize(X, ysparse, ydense, Z);
+ CheckKKT(sdp, X, ysparse, ydense, Z);
}
static inline arma::sp_mat
@@ -376,8 +420,8 @@ BOOST_AUTO_TEST_CASE(LogChebychevApproxSdp)
PrimalDualSolver<SDP<arma::sp_mat>> solver0(sdp0);
arma::mat X0, Z0;
arma::vec ysparse0, ydense0;
- const auto stat0 = solver0.Optimize(X0, ysparse0, ydense0, Z0);
- BOOST_REQUIRE(stat0.first);
+ solver0.Optimize(X0, ysparse0, ydense0, Z0);
+ CheckKKT(sdp0, X0, ysparse0, ydense0, Z0);
const size_t p1 = 10;
const size_t k1 = 5;
@@ -387,8 +431,8 @@ BOOST_AUTO_TEST_CASE(LogChebychevApproxSdp)
PrimalDualSolver<SDP<arma::sp_mat>> solver1(sdp1);
arma::mat X1, Z1;
arma::vec ysparse1, ydense1;
- const auto stat1 = solver1.Optimize(X1, ysparse1, ydense1, Z1);
- BOOST_REQUIRE(stat1.first);
+ solver1.Optimize(X1, ysparse1, ydense1, Z1);
+ CheckKKT(sdp1, X1, ysparse1, ydense1, Z1);
}
/**
@@ -495,9 +539,9 @@ BOOST_AUTO_TEST_CASE(CorrelationCoeffToySdp)
PrimalDualSolver<SDP<arma::sp_mat>> solver(sdp);
arma::mat X, Z;
arma::vec ysparse, ydense;
- const auto p = solver.Optimize(X, ysparse, ydense, Z);
- BOOST_REQUIRE(p.first);
- BOOST_REQUIRE_CLOSE(p.second, 2 * (-0.978), 1e-3);
+ const double obj = solver.Optimize(X, ysparse, ydense, Z);
+ CheckKKT(sdp, X, ysparse, ydense, Z);
+ BOOST_REQUIRE_CLOSE(obj, 2 * (-0.978), 1e-3);
}
///**
More information about the mlpack-git
mailing list