[mlpack-git] master: WIP: hammer out a few bugs, still doesn't work (a1f7aae)
gitdub at big.cc.gt.atl.ga.us
gitdub at big.cc.gt.atl.ga.us
Thu Mar 5 22:12:41 EST 2015
Repository : https://github.com/mlpack/mlpack
On branch : master
Link : https://github.com/mlpack/mlpack/compare/904762495c039e345beba14c1142fd719b3bd50e...f94823c800ad6f7266995c700b1b630d5ffdcf40
>---------------------------------------------------------------
commit a1f7aaed2eff62c80aeed8135510cf28b40f6999
Author: Stephen Tu <stephent at berkeley.edu>
Date: Thu Jan 15 16:37:58 2015 -0800
WIP: hammer out a few bugs, still doesn't work
>---------------------------------------------------------------
a1f7aaed2eff62c80aeed8135510cf28b40f6999
src/mlpack/core/math/lin_alg.cpp | 9 +--
src/mlpack/core/math/lin_alg.hpp | 11 +++-
src/mlpack/core/optimizers/sdp/primal_dual.cpp | 10 ++--
src/mlpack/tests/lin_alg_test.cpp | 76 ++++++++++++++++++++++++++
src/mlpack/tests/sdp_primal_dual_test.cpp | 2 +-
5 files changed, 97 insertions(+), 11 deletions(-)
diff --git a/src/mlpack/core/math/lin_alg.cpp b/src/mlpack/core/math/lin_alg.cpp
index 3e78fc8..a4ae426 100644
--- a/src/mlpack/core/math/lin_alg.cpp
+++ b/src/mlpack/core/math/lin_alg.cpp
@@ -270,7 +270,7 @@ void mlpack::math::SymKronId(const arma::mat& A, arma::mat& op)
size_t idx = 0;
for (size_t i = 0; i < n; i++)
{
- for (size_t j = 0; j < i; j++, idx++)
+ for (size_t j = i; j < n; j++)
{
for (size_t k = 0; k < n; k++)
{
@@ -278,10 +278,11 @@ void mlpack::math::SymKronId(const arma::mat& A, arma::mat& op)
((k == j) ? 1. : one_over_sq2) * A(i, k);
op(idx, SvecIndex(i, k, n)) +=
((k == i) ? 1. : one_over_sq2) * A(k, j);
- op.row(idx) *= 0.5;
- if (i != j)
- op.row(idx) *= sq2;
}
+ op.row(idx) *= 0.5;
+ if (i != j)
+ op.row(idx) *= sq2;
+ idx++;
}
}
}
diff --git a/src/mlpack/core/math/lin_alg.hpp b/src/mlpack/core/math/lin_alg.hpp
index 1f76f0c..1930a53 100644
--- a/src/mlpack/core/math/lin_alg.hpp
+++ b/src/mlpack/core/math/lin_alg.hpp
@@ -77,6 +77,10 @@ void RemoveRows(const arma::mat& input,
arma::mat& output);
/**
+ * Upper triangular representation of a symmetric matrix, scaled such that,
+ * dot(Svec(A), Svec(B)) == dot(A, B) for symmetric A, B. Specifically,
+ *
+ * Svec(K) = [ K_11, sqrt(2) K_12, ..., sqrt(2) K_1n, K_22, ..., sqrt(2) K_2n, ..., K_nn ]^T
*
* @param input A symmetric matrix
* @param output
@@ -84,7 +88,7 @@ void RemoveRows(const arma::mat& input,
void Svec(const arma::mat& input, arma::vec& output);
/**
- * The inverse of Svec
+ * The inverse of Svec. That is, Smat(Svec(A)) == A.
*
* @param input
* @param output A symmetric matrix
@@ -102,6 +106,11 @@ void Smat(const arma::vec& input, arma::mat& output);
size_t SvecIndex(size_t i, size_t j, size_t n);
/**
+ * If A is a symmetric matrix, then SymKronId returns an operator Op such that
+ *
+ * Op * svec(X) == svec(0.5 * (AX + XA))
+ *
+ * for every symmetric matrix X
*
* @param A
* @param op
diff --git a/src/mlpack/core/optimizers/sdp/primal_dual.cpp b/src/mlpack/core/optimizers/sdp/primal_dual.cpp
index 2dc0589..ebe100a 100644
--- a/src/mlpack/core/optimizers/sdp/primal_dual.cpp
+++ b/src/mlpack/core/optimizers/sdp/primal_dual.cpp
@@ -103,14 +103,14 @@ double PrimalDualSolver::Optimize(arma::mat& X,
for (size_t i = 0; i < sdp.NumSparseConstraints(); i++)
{
math::Svec(DenseFromSparse(sdp.SparseA()[i]), Ai);
- Asparse.row(i) = Ai;
+ Asparse.row(i) = Ai.t();
}
arma::mat Adense(sdp.NumDenseConstraints(), n2bar);
for (size_t i = 0; i < sdp.NumDenseConstraints(); i++)
{
math::Svec(sdp.DenseA()[i], Ai);
- Adense.row(i) = Ai;
+ Adense.row(i) = Ai.t();
}
arma::vec scsparse;
@@ -143,7 +143,7 @@ double PrimalDualSolver::Optimize(arma::mat& X,
rhs.set_size(sdp.NumConstraints());
Einv_F_AsparseT.set_size(n2bar, sdp.NumSparseConstraints());
- Einv_F_AsparseT.set_size(n2bar, sdp.NumDenseConstraints());
+ Einv_F_AdenseT.set_size(n2bar, sdp.NumDenseConstraints());
M.set_size(sdp.NumConstraints(), sdp.NumConstraints());
for (;;)
@@ -184,6 +184,7 @@ double PrimalDualSolver::Optimize(arma::mat& X,
Einv_F_AdenseT.col(i) = gk;
}
+
if (sdp.NumSparseConstraints())
{
M.submat(arma::span(0,
@@ -284,8 +285,7 @@ double PrimalDualSolver::Optimize(arma::mat& X,
const double duality_gap = dual_obj - primal_obj;
-
- Log::Info
+ Log::Debug
<< "primal=" << primal_obj << ", "
<< "dual=" << dual_obj << ", "
<< "gap=" << duality_gap << ", "
diff --git a/src/mlpack/tests/lin_alg_test.cpp b/src/mlpack/tests/lin_alg_test.cpp
index 3d3d9b6..e368b65 100644
--- a/src/mlpack/tests/lin_alg_test.cpp
+++ b/src/mlpack/tests/lin_alg_test.cpp
@@ -185,4 +185,80 @@ BOOST_AUTO_TEST_CASE(TestRemoveRows)
}
}
+BOOST_AUTO_TEST_CASE(TestSvecSmat)
+{
+ arma::mat X(3, 3);
+ X(0, 0) = 0; X(0, 1) = 1, X(0, 2) = 2;
+ X(1, 0) = 1; X(1, 1) = 3, X(1, 2) = 4;
+ X(2, 0) = 2; X(2, 1) = 4, X(2, 2) = 5;
+
+ arma::vec sx;
+ Svec(X, sx);
+ const double sq2 = sqrt(2.);
+ BOOST_REQUIRE_CLOSE(sx(0), 0, 1e-7);
+ BOOST_REQUIRE_CLOSE(sx(1), sq2 * 1., 1e-7);
+ BOOST_REQUIRE_CLOSE(sx(2), sq2 * 2., 1e-7);
+ BOOST_REQUIRE_CLOSE(sx(3), 3., 1e-7);
+ BOOST_REQUIRE_CLOSE(sx(4), sq2 * 4., 1e-7);
+ BOOST_REQUIRE_CLOSE(sx(5), 5., 1e-7);
+
+ arma::mat Xtest;
+ Smat(sx, Xtest);
+ BOOST_REQUIRE_EQUAL(Xtest.n_rows, 3);
+ BOOST_REQUIRE_EQUAL(Xtest.n_cols, 3);
+ for (size_t i = 0; i < 3; i++)
+ for (size_t j = 0; j < 3; j++)
+ BOOST_REQUIRE_CLOSE(X(i, j), Xtest(i, j), 1e-7);
+}
+
+BOOST_AUTO_TEST_CASE(TestSymKronIdSimple)
+{
+ arma::mat A(3, 3);
+ A(0, 0) = 1; A(0, 1) = 2, A(0, 2) = 3;
+ A(1, 0) = 2; A(1, 1) = 4, A(1, 2) = 5;
+ A(2, 0) = 3; A(2, 1) = 5, A(2, 2) = 6;
+ arma::mat Op;
+ SymKronId(A, Op);
+
+ const arma::mat X = A + arma::ones<arma::mat>(3, 3);
+ arma::vec sx;
+ Svec(X, sx);
+
+ const arma::vec lhs = Op * sx;
+ const arma::mat Rhs = 0.5 * (A * X + X * A);
+ arma::vec rhs;
+ Svec(Rhs, rhs);
+
+ BOOST_REQUIRE_EQUAL(lhs.n_elem, rhs.n_elem);
+ for (size_t j = 0; j < lhs.n_elem; j++)
+ BOOST_REQUIRE_CLOSE(lhs(j), rhs(j), 1e-5);
+}
+
+BOOST_AUTO_TEST_CASE(TestSymKronId)
+{
+ const size_t n = 10;
+ arma::mat A = arma::randu<arma::mat>(n, n);
+ A += A.t();
+
+ arma::mat Op;
+ SymKronId(A, Op);
+
+ for (size_t i = 0; i < 5; i++)
+ {
+ arma::mat X = arma::randu<arma::mat>(n, n);
+ X += X.t();
+ arma::vec sx;
+ Svec(X, sx);
+
+ const arma::vec lhs = Op * sx;
+ const arma::mat Rhs = 0.5 * (A * X + X * A);
+ arma::vec rhs;
+ Svec(Rhs, rhs);
+
+ BOOST_REQUIRE_EQUAL(lhs.n_elem, rhs.n_elem);
+ for (size_t j = 0; j < lhs.n_elem; j++)
+ BOOST_REQUIRE_CLOSE(lhs(j), rhs(j), 1e-5);
+ }
+}
+
BOOST_AUTO_TEST_SUITE_END();
diff --git a/src/mlpack/tests/sdp_primal_dual_test.cpp b/src/mlpack/tests/sdp_primal_dual_test.cpp
index a039c87..3f00cd4 100644
--- a/src/mlpack/tests/sdp_primal_dual_test.cpp
+++ b/src/mlpack/tests/sdp_primal_dual_test.cpp
@@ -154,7 +154,7 @@ BOOST_AUTO_TEST_CASE(SmallMaxCutFeasibleSdp)
ydense0.set_size(0);
X0.eye(g.NumVertices(), g.NumVertices());
- ysparse0 = -1.1 * arma::vec(arma::sum(arma::abs(sdp.SparseC()), 0));
+ ysparse0 = -1.1 * arma::vec(arma::sum(arma::abs(sdp.SparseC()), 0).t());
Z0 = -Diag(ysparse0) + sdp.SparseC();
PrimalDualSolver solver(sdp, X0, ysparse0, ydense0, Z0);
More information about the mlpack-git
mailing list