[mlpack-git] master: WIP: fix a few more bugs (c2df6c3)
gitdub at big.cc.gt.atl.ga.us
gitdub at big.cc.gt.atl.ga.us
Thu Mar 5 22:12:43 EST 2015
Repository : https://github.com/mlpack/mlpack
On branch : master
Link : https://github.com/mlpack/mlpack/compare/904762495c039e345beba14c1142fd719b3bd50e...f94823c800ad6f7266995c700b1b630d5ffdcf40
>---------------------------------------------------------------
commit c2df6c3c68dc94c9327bc142db698416264fdc9e
Author: Stephen Tu <stephent at berkeley.edu>
Date: Thu Jan 15 17:23:37 2015 -0800
WIP: fix a few more bugs
>---------------------------------------------------------------
c2df6c3c68dc94c9327bc142db698416264fdc9e
src/mlpack/core/optimizers/sdp/primal_dual.cpp | 53 ++++++++++++++++++++++----
src/mlpack/core/optimizers/sdp/primal_dual.hpp | 3 ++
src/mlpack/tests/sdp_primal_dual_test.cpp | 35 ++++++++++++++---
3 files changed, 78 insertions(+), 13 deletions(-)
diff --git a/src/mlpack/core/optimizers/sdp/primal_dual.cpp b/src/mlpack/core/optimizers/sdp/primal_dual.cpp
index ebe100a..36c8ba8 100644
--- a/src/mlpack/core/optimizers/sdp/primal_dual.cpp
+++ b/src/mlpack/core/optimizers/sdp/primal_dual.cpp
@@ -13,7 +13,8 @@ PrimalDualSolver::PrimalDualSolver(const SDP& sdp)
tau(0.5),
normXzTol(1e-7),
primalInfeasTol(1e-7),
- dualInfeasTol(1e-7)
+ dualInfeasTol(1e-7),
+ maxIterations(1000)
{
}
@@ -32,7 +33,8 @@ PrimalDualSolver::PrimalDualSolver(const SDP& sdp,
tau(0.5),
normXzTol(1e-7),
primalInfeasTol(1e-7),
- dualInfeasTol(1e-7)
+ dualInfeasTol(1e-7),
+ maxIterations(1000)
{
if (X0.n_rows != sdp.N() || X0.n_cols != sdp.N())
Log::Fatal << "PrimalDualSolver::PrimalDualSolver(): "
@@ -136,7 +138,7 @@ double PrimalDualSolver::Optimize(arma::mat& X,
arma::vec rp, rd, rc, gk, Einv_Frd_rc, Einv_Frd_ATdy_rc, rhs;
arma::mat Rc, E, F, Einv_F_AsparseT, Einv_F_AdenseT, Gk,
- M, Einv_Frd_rc_Mat, Frd_rc_Mat, Frd_ATdy_Mat,
+ M, Einv_Frd_rc_Mat, Frd_rc_Mat, Frd_ATdy_rc_Mat,
Einv_Frd_ATdy_rc_Mat;
rp.set_size(sdp.NumConstraints());
@@ -146,7 +148,7 @@ double PrimalDualSolver::Optimize(arma::mat& X,
Einv_F_AdenseT.set_size(n2bar, sdp.NumDenseConstraints());
M.set_size(sdp.NumConstraints(), sdp.NumConstraints());
- for (;;)
+ for (size_t iteration = 0; iteration < maxIterations; iteration++)
{
const double mu = sigma * arma::dot(sx, sz) / n;
@@ -158,18 +160,33 @@ double PrimalDualSolver::Optimize(arma::mat& X,
rp(arma::span(sdp.NumSparseConstraints(), sdp.NumConstraints() - 1)) =
sdp.DenseB() - Adense * sx;
+ //std::cout << "rp" << std::endl;
+ //std::cout << rp << std::endl;
+
rd = - sz - Asparse.t() * ysparse - Adense.t() * ydense;
if (sdp.HasSparseObjective())
rd += scsparse;
if (sdp.HasDenseObjective())
rd += scdense;
+ //std::cout << "rd" << std::endl;
+ //std::cout << rd << std::endl;
+
Rc = mu*arma::eye<arma::mat>(n, n) - 0.5*(X*Z + Z*X);
math::Svec(Rc, rc);
+ //std::cout << "rc" << std::endl;
+ //std::cout << rc << std::endl;
+
math::SymKronId(Z, E);
math::SymKronId(X, F);
+ //std::cout << "E" << std::endl;
+ //std::cout << E << std::endl;
+
+ //std::cout << "F" << std::endl;
+ //std::cout << F << std::endl;
+
for (size_t i = 0; i < sdp.NumSparseConstraints(); i++)
{
SolveLyapunov(Gk, Z, X * sdp.SparseA()[i] + sdp.SparseA()[i] * X);
@@ -177,6 +194,9 @@ double PrimalDualSolver::Optimize(arma::mat& X,
Einv_F_AsparseT.col(i) = gk;
}
+ //std::cout << "Einv_F_AsparseT" << std::endl;
+ //std::cout << Einv_F_AsparseT << std::endl;
+
for (size_t i = 0; i < sdp.NumDenseConstraints(); i++)
{
SolveLyapunov(Gk, Z, X * sdp.DenseA()[i] + sdp.DenseA()[i] * X);
@@ -184,7 +204,6 @@ double PrimalDualSolver::Optimize(arma::mat& X,
Einv_F_AdenseT.col(i) = gk;
}
-
if (sdp.NumSparseConstraints())
{
M.submat(arma::span(0,
@@ -218,6 +237,8 @@ double PrimalDualSolver::Optimize(arma::mat& X,
Adense * Einv_F_AdenseT;
}
+ //std::cout << "M" << std::endl;
+ //std::cout << M << std::endl;
math::Smat(F * rd - rc, Frd_rc_Mat);
SolveLyapunov(Einv_Frd_rc_Mat, Z, 2. * Frd_rc_Mat);
@@ -230,18 +251,30 @@ double PrimalDualSolver::Optimize(arma::mat& X,
rhs(arma::span(sdp.NumSparseConstraints(), sdp.NumConstraints() - 1)) += Adense * Einv_Frd_rc;
// TODO(stephentu): use a more efficient method (e.g. LU decomposition)
+ //std::cout << "rhs" << std::endl;
+ //std::cout << rhs << std::endl;
+
arma::solve(dy, M, rhs);
if (sdp.NumSparseConstraints())
dysparse = dy(arma::span(0, sdp.NumSparseConstraints() - 1));
if (sdp.NumDenseConstraints())
dydense = dy(arma::span(sdp.NumSparseConstraints(), sdp.NumConstraints() - 1));
+ //std::cout << "dy" << std::endl;
+ //std::cout << dy << std::endl;
- math::Smat(F * (rd - Asparse.t() * dysparse - Adense.t() * dydense), Frd_ATdy_Mat);
- SolveLyapunov(Einv_Frd_ATdy_rc_Mat, Z, 2. * Frd_ATdy_Mat);
+ math::Smat(F * (rd - Asparse.t() * dysparse - Adense.t() * dydense) - rc, Frd_ATdy_rc_Mat);
+ SolveLyapunov(Einv_Frd_ATdy_rc_Mat, Z, 2. * Frd_ATdy_rc_Mat);
math::Svec(Einv_Frd_ATdy_rc_Mat, Einv_Frd_ATdy_rc);
dsx = -Einv_Frd_ATdy_rc;
+
+ //std::cout << "dsx" << std::endl;
+ //std::cout << dsx << std::endl;
+
dsz = rd - Asparse.t() * dysparse - Adense.t() * dydense;
+ //std::cout << "dsz" << std::endl;
+ //std::cout << dsz << std::endl;
+
math::Smat(dsx, dX);
math::Smat(dsz, dZ);
@@ -283,9 +316,10 @@ double PrimalDualSolver::Optimize(arma::mat& X,
arma::dot(sdp.SparseB(), ysparse) +
arma::dot(sdp.DenseB(), ydense);
- const double duality_gap = dual_obj - primal_obj;
+ const double duality_gap = primal_obj - dual_obj;
Log::Debug
+ << "iter=" << iteration + 1 << ", "
<< "primal=" << primal_obj << ", "
<< "dual=" << dual_obj << ", "
<< "gap=" << duality_gap << ", "
@@ -298,6 +332,9 @@ double PrimalDualSolver::Optimize(arma::mat& X,
primal_infeas <= primalInfeasTol)
return primal_obj;
}
+
+ Log::Warn << "did not converge" << std::endl;
+ return -1;
}
} // namespace optimization
diff --git a/src/mlpack/core/optimizers/sdp/primal_dual.hpp b/src/mlpack/core/optimizers/sdp/primal_dual.hpp
index 9f4d038..d01a177 100644
--- a/src/mlpack/core/optimizers/sdp/primal_dual.hpp
+++ b/src/mlpack/core/optimizers/sdp/primal_dual.hpp
@@ -39,6 +39,7 @@ class PrimalDualSolver {
double& NormXzTol() { return normXzTol; }
double& PrimalInfeasTol() { return primalInfeasTol; }
double& DualInfeasTol() { return dualInfeasTol; }
+ size_t& MaxIterations() { return maxIterations; }
private:
SDP sdp;
@@ -54,6 +55,8 @@ class PrimalDualSolver {
double primalInfeasTol;
double dualInfeasTol;
+ size_t maxIterations;
+
};
} // namespace optimization
diff --git a/src/mlpack/tests/sdp_primal_dual_test.cpp b/src/mlpack/tests/sdp_primal_dual_test.cpp
index 3f00cd4..bd4ff16 100644
--- a/src/mlpack/tests/sdp_primal_dual_test.cpp
+++ b/src/mlpack/tests/sdp_primal_dual_test.cpp
@@ -111,7 +111,7 @@ class UndirectedGraph
};
static inline SDP
-ConstructMaxCutSDP(const UndirectedGraph& g)
+ConstructMaxCutSDPFromGraph(const UndirectedGraph& g)
{
SDP sdp(g.NumVertices(), g.NumVertices(), 0);
g.Laplacian(sdp.SparseC());
@@ -138,6 +138,25 @@ Diag(const arma::vec& diag)
return ret;
}
+static inline SDP
+ConstructMaxCutSDPFromLaplacian(const std::string& laplacianFilename)
+{
+ arma::mat laplacian;
+ data::Load(laplacianFilename, laplacian, false);
+ if (laplacian.n_rows != laplacian.n_cols)
+ Log::Fatal << "laplacian not square" << std::endl;
+ SDP sdp(laplacian.n_rows, laplacian.n_rows, 0);
+ sdp.SparseC() = -arma::sp_mat(laplacian);
+ for (size_t i = 0; i < laplacian.n_rows; i++)
+ {
+ sdp.SparseA()[i].zeros(laplacian.n_rows, laplacian.n_rows);
+ sdp.SparseA()[i](i, i) = 1.;
+ }
+ sdp.SparseB().ones();
+ return sdp;
+}
+
+
BOOST_AUTO_TEST_SUITE(SdpPrimalDualTest);
/**
@@ -145,19 +164,25 @@ BOOST_AUTO_TEST_SUITE(SdpPrimalDualTest);
*/
BOOST_AUTO_TEST_CASE(SmallMaxCutFeasibleSdp)
{
- UndirectedGraph g;
- UndirectedGraph::ErdosRenyiRandomGraph(g, 10, 0.3, true);
- SDP sdp = ConstructMaxCutSDP(g);
+ //UndirectedGraph g;
+ //UndirectedGraph::ErdosRenyiRandomGraph(g, 10, 0.3, true);
+ //SDP sdp = ConstructMaxCutSDPFromGraph(g);
+
+ SDP sdp = ConstructMaxCutSDPFromLaplacian("r10.txt");
arma::mat X0, Z0;
arma::vec ysparse0, ydense0;
ydense0.set_size(0);
- X0.eye(g.NumVertices(), g.NumVertices());
+ X0.eye(sdp.N(), sdp.N());
ysparse0 = -1.1 * arma::vec(arma::sum(arma::abs(sdp.SparseC()), 0).t());
Z0 = -Diag(ysparse0) + sdp.SparseC();
+ //std::cout << "ysparse0" << std::endl;
+ //std::cout << ysparse0 << std::endl;
+
PrimalDualSolver solver(sdp, X0, ysparse0, ydense0, Z0);
+ //solver.MaxIterations() = 1;
arma::mat X, Z;
arma::vec ysparse, ydense;
More information about the mlpack-git
mailing list