[mlpack-git] master: implement the mehrotra's predictor-corrector rule (159d581)

[gitdub at big.cc.gt.atl.ga.us]

Repository : https://github.com/mlpack/mlpack

On branch  : master
Link       : https://github.com/mlpack/mlpack/compare/bb6e5c56aab07e6449d86021246b52a9e323f3a0...bd6cb33f8d8270b02a84e81e38727679bb6c319a

>---------------------------------------------------------------

commit 159d58177a829ec4bbea6c18102f3501040cd47f
Author: Stephen Tu <tu.stephenl at gmail.com>
Date:   Wed Jan 21 00:23:26 2015 -0800

    implement the mehrotra's predictor-corrector rule
    
    still need to remove the code duplication


>---------------------------------------------------------------

159d58177a829ec4bbea6c18102f3501040cd47f
 src/mlpack/core/optimizers/sdp/primal_dual.cpp | 73 +++++++++++++++++++++-----
 src/mlpack/core/optimizers/sdp/primal_dual.hpp |  2 -
 2 files changed, 61 insertions(+), 14 deletions(-)

diff --git a/src/mlpack/core/optimizers/sdp/primal_dual.cpp b/src/mlpack/core/optimizers/sdp/primal_dual.cpp
index 062b875..abaebc0 100644
--- a/src/mlpack/core/optimizers/sdp/primal_dual.cpp
+++ b/src/mlpack/core/optimizers/sdp/primal_dual.cpp
@@ -28,8 +28,7 @@ PrimalDualSolver::PrimalDualSolver(const SDP& sdp)
     ysparse0(arma::ones<arma::vec>(sdp.NumSparseConstraints())),
     ydense0(arma::ones<arma::vec>(sdp.NumDenseConstraints())),
     Z0(arma::eye<arma::mat>(sdp.N(), sdp.N())),
-    sigma(0.5),
-    tau(0.5),
+    tau(0.99),
     normXzTol(1e-7),
     primalInfeasTol(1e-7),
     dualInfeasTol(1e-7),
@@ -48,8 +47,7 @@ PrimalDualSolver::PrimalDualSolver(const SDP& sdp,
     ysparse0(ysparse0),
     ydense0(ydense0),
     Z0(Z0),
-    sigma(0.5),
-    tau(0.5),
+    tau(0.99),
     normXzTol(1e-7),
     primalInfeasTol(1e-7),
     dualInfeasTol(1e-7),
@@ -191,7 +189,7 @@ PrimalDualSolver::Optimize(arma::mat& X,
   for (size_t iteration = 0; iteration < maxIterations; iteration++)
   {
 
-    const double mu = sigma * arma::dot(sx, sz) / n;
+
 
     if (sdp.NumSparseConstraints())
       rp(arma::span(0, sdp.NumSparseConstraints() - 1)) =
@@ -206,8 +204,7 @@ PrimalDualSolver::Optimize(arma::mat& X,
     if (sdp.HasDenseObjective())
       rd += scdense;
 
-    Rc = mu*arma::eye<arma::mat>(n, n) - 0.5*(X*Z + Z*X);
-    math::Svec(Rc, rc);
+
 
     math::SymKronId(Z, E);
     math::SymKronId(X, F);
@@ -259,6 +256,63 @@ PrimalDualSolver::Optimize(arma::mat& X,
         Adense * Einv_F_AdenseT;
     }
 
+    const double sxdotsz = arma::dot(sx, sz);
+    double sigma = 0.;
+    {
+      Rc = -0.5*(X*Z + Z*X);
+      math::Svec(Rc, rc);
+      math::Smat(F * rd - rc, Frd_rc_Mat);
+      SolveLyapunov(Einv_Frd_rc_Mat, Z, 2. * Frd_rc_Mat);
+      math::Svec(Einv_Frd_rc_Mat, Einv_Frd_rc);
+
+      rhs = rp;
+      if (sdp.NumSparseConstraints())
+        rhs(arma::span(0, sdp.NumSparseConstraints() - 1)) += Asparse * Einv_Frd_rc;
+      if (sdp.NumDenseConstraints())
+        rhs(arma::span(sdp.NumSparseConstraints(), sdp.NumConstraints() - 1)) += Adense * Einv_Frd_rc;
+
+      // 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;
+
+      if (sdp.NumSparseConstraints())
+        dysparse = dy(arma::span(0, sdp.NumSparseConstraints() - 1));
+      if (sdp.NumDenseConstraints())
+        dydense = dy(arma::span(sdp.NumSparseConstraints(), sdp.NumConstraints() - 1));
+
+      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;
+      dsz = rd - Asparse.t() * dysparse - Adense.t() * dydense;
+
+      math::Smat(dsx, dX);
+      math::Smat(dsz, dZ);
+
+      // TODO(stephentu): computing these alphahats should take advantage of
+      // the cholesky decomposition of X and Z which we should have available
+      // when we use more efficient methods above.
+
+      double alphahatX = AlphaHat(X, dX);
+      if (alphahatX < 0.)
+        // dX is PSD
+        alphahatX = 1.;
+
+      double alphahatZ = AlphaHat(Z, dZ);
+      if (alphahatZ < 0.)
+        // dZ is PSD
+        alphahatZ = 1.;
+
+      const double alpha = std::min(1., tau * alphahatX);
+      const double beta = std::min(1., tau * alphahatZ);
+
+      sigma = std::pow(arma::dot(X + alpha * dX, Z + beta * dZ) / sxdotsz, 3);
+    }
+
+    const double mu = sigma * sxdotsz / n;
+
+    Rc = mu*arma::eye<arma::mat>(n, n) - 0.5*(X*Z + Z*X + dX*dZ + dZ*dX);
+    math::Svec(Rc, rc);
     math::Smat(F * rd - rc, Frd_rc_Mat);
     SolveLyapunov(Einv_Frd_rc_Mat, Z, 2. * Frd_rc_Mat);
     math::Svec(Einv_Frd_rc_Mat, Einv_Frd_rc);
@@ -304,11 +358,6 @@ PrimalDualSolver::Optimize(arma::mat& X,
     const double alpha = std::min(1., tau * alphahatX);
     const double beta = std::min(1., tau * alphahatZ);
 
-    // TODO(stephentu): Implement the Mehrotra's predictor-corrector rule.  See
-    // Section 7 of [AHO98]. This will require making the above KKT system
-    // solver more modular (since we'll have to solve another similar KKT
-    // system).
-
     X += alpha * dX;
     math::Svec(X, sx);
     ysparse += beta * dysparse;
diff --git a/src/mlpack/core/optimizers/sdp/primal_dual.hpp b/src/mlpack/core/optimizers/sdp/primal_dual.hpp
index 8c3b863..9480b20 100644
--- a/src/mlpack/core/optimizers/sdp/primal_dual.hpp
+++ b/src/mlpack/core/optimizers/sdp/primal_dual.hpp
@@ -37,7 +37,6 @@ class PrimalDualSolver {
     return Optimize(X, ysparse, ydense, Z);
   }
 
-  double& Sigma() { return sigma; }
   double& Tau() { return tau; }
   double& NormXzTol() { return normXzTol; }
   double& PrimalInfeasTol() { return primalInfeasTol; }
@@ -52,7 +51,6 @@ class PrimalDualSolver {
   arma::vec ydense0;
   arma::mat Z0;
 
-  double sigma;
   double tau;
   double normXzTol;
   double primalInfeasTol;