[mlpack-git] master: WIP: a compiling, but non-functioning primal-dual IP SDP solver (ec63d92)

gitdub at big.cc.gt.atl.ga.us gitdub at big.cc.gt.atl.ga.us
Thu Mar 5 22:12:37 EST 2015


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

On branch  : master
Link       : https://github.com/mlpack/mlpack/compare/904762495c039e345beba14c1142fd719b3bd50e...f94823c800ad6f7266995c700b1b630d5ffdcf40

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

commit ec63d92ccb6fe98143d5632bd0b1fb86fdcc72d1
Author: Stephen Tu <tu.stephenl at gmail.com>
Date:   Thu Jan 15 12:53:09 2015 -0800

    WIP: a compiling, but non-functioning primal-dual IP SDP solver


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

ec63d92ccb6fe98143d5632bd0b1fb86fdcc72d1
 src/mlpack/core/math/lin_alg.cpp                   |  77 ++++++
 src/mlpack/core/math/lin_alg.hpp                   |  32 +++
 src/mlpack/core/optimizers/CMakeLists.txt          |   1 +
 .../core/optimizers/{sa => sdp}/CMakeLists.txt     |   7 +-
 src/mlpack/core/optimizers/sdp/primal_dual.cpp     | 304 +++++++++++++++++++++
 src/mlpack/core/optimizers/sdp/primal_dual.hpp     |  62 +++++
 src/mlpack/core/optimizers/sdp/sdp.cpp             |  29 ++
 src/mlpack/core/optimizers/sdp/sdp.hpp             | 114 ++++++++
 8 files changed, 623 insertions(+), 3 deletions(-)

diff --git a/src/mlpack/core/math/lin_alg.cpp b/src/mlpack/core/math/lin_alg.cpp
index d4b49f3..3e78fc8 100644
--- a/src/mlpack/core/math/lin_alg.cpp
+++ b/src/mlpack/core/math/lin_alg.cpp
@@ -208,3 +208,80 @@ void mlpack::math::RemoveRows(const arma::mat& input,
     }
   }
 }
+
+static const double sq2 = 1.41421356237309504880; // sqrt(2)
+static const double one_over_sq2 = 0.707106781186547524401; // 1/sqrt(2)
+
+void mlpack::math::Svec(const arma::mat& input, arma::vec& output)
+{
+  const size_t n = input.n_rows;
+  const size_t n2bar = n * (n + 1) / 2;
+
+  output.zeros(n2bar);
+
+  size_t idx = 0;
+  for (size_t i = 0; i < n; i++)
+  {
+    for (size_t j = i; j < n; j++)
+    {
+      if (i == j)
+        output(idx++) = input(i, j);
+      else
+        output(idx++) = sq2 * input(i, j);
+    }
+  }
+}
+
+void mlpack::math::Smat(const arma::vec& input, arma::mat& output)
+{
+  const size_t n = static_cast<size_t>(ceil((-1. + sqrt(1. + 8. * input.n_elem))/2.));
+
+  output.zeros(n, n);
+
+  size_t idx = 0;
+  for (size_t i = 0; i < n; i++)
+  {
+    for (size_t j = i; j < n; j++)
+    {
+      if (i == j)
+        output(i, j) = input(idx++);
+      else
+        output(i, j) = output(j, i) = one_over_sq2 * input(idx++);
+    }
+  }
+}
+
+size_t mlpack::math::SvecIndex(size_t i, size_t j, size_t n)
+{
+  if (i > j)
+    std::swap(i, j);
+  return (j-i) + (n*(n+1) - (n-i)*(n-i+1))/2;
+}
+
+
+void mlpack::math::SymKronId(const arma::mat& A, arma::mat& op)
+{
+  // TODO(stephentu): there's probably an easier way to build this operator
+
+  const size_t n = A.n_rows;
+  const size_t n2bar = n * (n + 1) / 2;
+  op.zeros(n2bar, n2bar);
+
+  size_t idx = 0;
+  for (size_t i = 0; i < n; i++)
+  {
+    for (size_t j = 0; j < i; j++, idx++)
+    {
+      for (size_t k = 0; k < n; k++)
+      {
+        op(idx, SvecIndex(k, j, n)) +=
+          ((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;
+      }
+    }
+  }
+}
diff --git a/src/mlpack/core/math/lin_alg.hpp b/src/mlpack/core/math/lin_alg.hpp
index 1988beb..1f76f0c 100644
--- a/src/mlpack/core/math/lin_alg.hpp
+++ b/src/mlpack/core/math/lin_alg.hpp
@@ -76,6 +76,38 @@ void RemoveRows(const arma::mat& input,
                 const std::vector<size_t>& rowsToRemove,
                 arma::mat& output);
 
+/**
+ *
+ * @param input A symmetric matrix
+ * @param output
+ */
+void Svec(const arma::mat& input, arma::vec& output);
+
+/**
+ * The inverse of Svec
+ *
+ * @param input
+ * @param output A symmetric matrix
+ */
+void Smat(const arma::vec& input, arma::mat& output);
+
+/**
+ * Return the index such that A[i,j] == factr(i, j) * svec(A)[pos(i, j)],
+ * where factr(i, j) = sqrt(2) if i != j and 1 otherwise.
+ *
+ * @param i
+ * @param j
+ * @param n
+ */
+size_t SvecIndex(size_t i, size_t j, size_t n);
+
+/**
+ *
+ * @param A
+ * @param op
+ */
+void SymKronId(const arma::mat& A, arma::mat& op);
+
 }; // namespace math
 }; // namespace mlpack
 
diff --git a/src/mlpack/core/optimizers/CMakeLists.txt b/src/mlpack/core/optimizers/CMakeLists.txt
index 1243586..3579a2a 100644
--- a/src/mlpack/core/optimizers/CMakeLists.txt
+++ b/src/mlpack/core/optimizers/CMakeLists.txt
@@ -3,6 +3,7 @@ set(DIRS
   lbfgs
   lrsdp
   sa
+  sdp
   sgd
 )
 
diff --git a/src/mlpack/core/optimizers/sa/CMakeLists.txt b/src/mlpack/core/optimizers/sdp/CMakeLists.txt
similarity index 77%
copy from src/mlpack/core/optimizers/sa/CMakeLists.txt
copy to src/mlpack/core/optimizers/sdp/CMakeLists.txt
index 7ee164a..cba5d07 100644
--- a/src/mlpack/core/optimizers/sa/CMakeLists.txt
+++ b/src/mlpack/core/optimizers/sdp/CMakeLists.txt
@@ -1,7 +1,8 @@
 set(SOURCES
-  sa.hpp
-  sa_impl.hpp
-  exponential_schedule.hpp
+  sdp.hpp
+  sdp.cpp
+  primal_dual.hpp
+  primal_dual.cpp
 )
 
 set(DIR_SRCS)
diff --git a/src/mlpack/core/optimizers/sdp/primal_dual.cpp b/src/mlpack/core/optimizers/sdp/primal_dual.cpp
new file mode 100644
index 0000000..2dc0589
--- /dev/null
+++ b/src/mlpack/core/optimizers/sdp/primal_dual.cpp
@@ -0,0 +1,304 @@
+#include "primal_dual.hpp"
+
+namespace mlpack {
+namespace optimization {
+
+PrimalDualSolver::PrimalDualSolver(const SDP& sdp)
+  : sdp(sdp),
+    X0(arma::eye<arma::mat>(sdp.N(), sdp.N())),
+    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),
+    normXzTol(1e-7),
+    primalInfeasTol(1e-7),
+    dualInfeasTol(1e-7)
+{
+
+}
+
+PrimalDualSolver::PrimalDualSolver(const SDP& sdp,
+                                   const arma::mat& X0,
+                                   const arma::vec& ysparse0,
+                                   const arma::vec& ydense0,
+                                   const arma::mat& Z0)
+  : sdp(sdp),
+    X0(X0),
+    ysparse0(ysparse0),
+    ydense0(ydense0),
+    Z0(Z0),
+    sigma(0.5),
+    tau(0.5),
+    normXzTol(1e-7),
+    primalInfeasTol(1e-7),
+    dualInfeasTol(1e-7)
+{
+  if (X0.n_rows != sdp.N() || X0.n_cols != sdp.N())
+    Log::Fatal << "PrimalDualSolver::PrimalDualSolver(): "
+      << "X0 needs to be square n x n matrix"
+      << std::endl;
+
+  if (ysparse0.n_elem != sdp.NumSparseConstraints())
+    Log::Fatal << "PrimalDualSolver::PrimalDualSolver(): "
+      << "ysparse0 needs to have the same length as the number of sparse constraints"
+      << std::endl;
+
+  if (ydense0.n_elem != sdp.NumDenseConstraints())
+    Log::Fatal << "PrimalDualSolver::PrimalDualSolver(): "
+      << "ydense0 needs to have the same length as the number of dense constraints"
+      << std::endl;
+
+  if (Z0.n_rows != sdp.N() || Z0.n_cols != sdp.N())
+    Log::Fatal << "PrimalDualSolver::PrimalDualSolver(): "
+      << "Z0 needs to be square n x n matrix"
+      << std::endl;
+}
+
+static inline arma::mat
+DenseFromSparse(const arma::sp_mat& input)
+{
+  return arma::mat(input);
+}
+
+static inline double
+AlphaHat(const arma::mat& A, const arma::mat& dA)
+{
+  // note: arma::chol(A) returns an upper triangular matrix (instead of the
+  // usual lower triangular)
+  const arma::mat L = arma::trimatl(arma::chol(A).t());
+  const arma::mat Linv = L.i();
+  const arma::vec evals = arma::eig_sym(-Linv * dA * Linv.t());
+  const double alphahatinv = evals(evals.n_elem - 1);
+  return 1. / alphahatinv;
+}
+
+/**
+ * Solve the following Lyapunov equation (for X)
+ *
+ *   AX + XA = H
+ *
+ * where A, H are symmetric matrices
+ *
+ */
+static inline void
+SolveLyapunov(arma::mat& X, const arma::mat& A, const arma::mat& H)
+{
+  arma::syl(X, A, A, -H);
+}
+
+double PrimalDualSolver::Optimize(arma::mat& X,
+                                  arma::vec& ysparse,
+                                  arma::vec& ydense,
+                                  arma::mat& Z)
+{
+  const size_t n = sdp.N();
+  const size_t n2bar = sdp.N2bar();
+
+  // TODO: implementation does not take adv of sparsity yet
+
+  arma::mat Asparse(sdp.NumSparseConstraints(), n2bar);
+  arma::vec Ai;
+
+  for (size_t i = 0; i < sdp.NumSparseConstraints(); i++)
+  {
+    math::Svec(DenseFromSparse(sdp.SparseA()[i]), Ai);
+    Asparse.row(i) = Ai;
+  }
+
+  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;
+  }
+
+  arma::vec scsparse;
+  if (sdp.HasSparseObjective())
+    math::Svec(DenseFromSparse(sdp.SparseC()), scsparse);
+
+  arma::vec scdense;
+  if (sdp.HasDenseObjective())
+    math::Svec(sdp.DenseC(), scdense);
+
+
+  X = X0;
+  ysparse = ysparse0;
+  ydense = ydense0;
+  Z = Z0;
+
+  arma::vec sx, sz, dy, dysparse, dydense, dsx, dsz;
+  arma::mat dX, dZ;
+
+  math::Svec(X, sx);
+  math::Svec(Z, sz);
+
+  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,
+            Einv_Frd_ATdy_rc_Mat;
+
+  rp.set_size(sdp.NumConstraints());
+  rhs.set_size(sdp.NumConstraints());
+
+  Einv_F_AsparseT.set_size(n2bar, sdp.NumSparseConstraints());
+  Einv_F_AsparseT.set_size(n2bar, sdp.NumDenseConstraints());
+  M.set_size(sdp.NumConstraints(), sdp.NumConstraints());
+
+  for (;;)
+  {
+
+    const double mu = sigma * arma::dot(sx, sz) / n;
+
+    if (sdp.NumSparseConstraints())
+      rp(arma::span(0, sdp.NumSparseConstraints() - 1)) =
+        sdp.SparseB() - Asparse * sx;
+    if (sdp.NumDenseConstraints())
+      rp(arma::span(sdp.NumSparseConstraints(), sdp.NumConstraints() - 1)) =
+          sdp.DenseB() - Adense * sx;
+
+    rd = - sz - Asparse.t() * ysparse - Adense.t() * ydense;
+    if (sdp.HasSparseObjective())
+      rd += scsparse;
+    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);
+
+    for (size_t i = 0; i < sdp.NumSparseConstraints(); i++)
+    {
+      SolveLyapunov(Gk, Z, X * sdp.SparseA()[i] + sdp.SparseA()[i] * X);
+      math::Svec(Gk, gk);
+      Einv_F_AsparseT.col(i) = gk;
+    }
+
+    for (size_t i = 0; i < sdp.NumDenseConstraints(); i++)
+    {
+      SolveLyapunov(Gk, Z, X * sdp.DenseA()[i] + sdp.DenseA()[i] * X);
+      math::Svec(Gk, gk);
+      Einv_F_AdenseT.col(i) = gk;
+    }
+
+    if (sdp.NumSparseConstraints())
+    {
+      M.submat(arma::span(0,
+                          sdp.NumSparseConstraints() - 1),
+               arma::span(0,
+                          sdp.NumSparseConstraints() - 1)) =
+        Asparse * Einv_F_AsparseT;
+      if (sdp.NumDenseConstraints())
+      {
+        M.submat(arma::span(0,
+                            sdp.NumSparseConstraints() - 1),
+                 arma::span(sdp.NumSparseConstraints(),
+                            sdp.NumConstraints() - 1)) =
+          Asparse * Einv_F_AdenseT;
+      }
+    }
+    if (sdp.NumDenseConstraints())
+    {
+      if (sdp.NumSparseConstraints())
+      {
+        M.submat(arma::span(sdp.NumSparseConstraints(),
+                            sdp.NumConstraints() - 1),
+                 arma::span(0,
+                            sdp.NumSparseConstraints() - 1)) =
+          Adense * Einv_F_AsparseT;
+      }
+      M.submat(arma::span(sdp.NumSparseConstraints(),
+                          sdp.NumConstraints() - 1),
+               arma::span(sdp.NumSparseConstraints(),
+                          sdp.NumConstraints() - 1)) =
+        Adense * Einv_F_AdenseT;
+    }
+
+
+    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)
+    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));
+
+    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::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);
+
+    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);
+
+    X += alpha * dX;
+    math::Svec(X, sx);
+    ysparse += beta * dysparse;
+    ydense += beta * dydense;
+    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);
+
+    double primal_obj = 0.;
+    if (sdp.HasSparseObjective())
+      primal_obj += arma::dot(sdp.SparseC(), X);
+    if (sdp.HasDenseObjective())
+      primal_obj += arma::dot(sdp.DenseC(), X);
+
+    const double dual_obj =
+      arma::dot(sdp.SparseB(), ysparse) +
+      arma::dot(sdp.DenseB(), ydense);
+
+    const double duality_gap = dual_obj - primal_obj;
+
+
+    Log::Info
+      << "primal=" << primal_obj << ", "
+      << "dual=" << dual_obj << ", "
+      << "gap=" << duality_gap << ", "
+      << "||XZ||=" << norm_XZ << ", "
+      << "primal_infeas=" << primal_infeas << ", "
+      << "mu=" << mu
+      << std::endl;
+
+    if (norm_XZ <= normXzTol &&
+        primal_infeas <= primalInfeasTol)
+      return primal_obj;
+  }
+}
+
+} // namespace optimization
+} // namespace mlpack
diff --git a/src/mlpack/core/optimizers/sdp/primal_dual.hpp b/src/mlpack/core/optimizers/sdp/primal_dual.hpp
new file mode 100644
index 0000000..9f4d038
--- /dev/null
+++ b/src/mlpack/core/optimizers/sdp/primal_dual.hpp
@@ -0,0 +1,62 @@
+/**
+ * @file primal_dual.hpp
+ * @author Stephen Tu
+ *
+ */
+#ifndef __MLPACK_CORE_OPTIMIZERS_SDP_PRIMAL_DUAL_HPP
+#define __MLPACK_CORE_OPTIMIZERS_SDP_PRIMAL_DUAL_HPP
+
+#include <mlpack/core.hpp>
+#include <mlpack/core/optimizers/sdp/sdp.hpp>
+
+namespace mlpack {
+namespace optimization {
+
+class PrimalDualSolver {
+ public:
+
+  PrimalDualSolver(const SDP& sdp);
+
+  PrimalDualSolver(const SDP& sdp,
+                   const arma::mat& X0,
+                   const arma::vec& ysparse0,
+                   const arma::vec& ydense0,
+                   const arma::mat& Z0);
+
+  double Optimize(arma::mat& X,
+                  arma::vec& ysparse,
+                  arma::vec& ydense,
+                  arma::mat &Z);
+
+  double Optimize(arma::mat& X) {
+    arma::vec ysparse, ydense;
+    arma::mat Z;
+    return Optimize(X, ysparse, ydense, Z);
+  }
+
+  double& Sigma() { return sigma; }
+  double& Tau() { return tau; }
+  double& NormXzTol() { return normXzTol; }
+  double& PrimalInfeasTol() { return primalInfeasTol; }
+  double& DualInfeasTol() { return dualInfeasTol; }
+
+ private:
+  SDP sdp;
+
+  arma::mat X0;
+  arma::vec ysparse0;
+  arma::vec ydense0;
+  arma::mat Z0;
+
+  double sigma;
+  double tau;
+  double normXzTol;
+  double primalInfeasTol;
+  double dualInfeasTol;
+
+};
+
+} // namespace optimization
+} // namespace mlpack
+
+#endif
diff --git a/src/mlpack/core/optimizers/sdp/sdp.cpp b/src/mlpack/core/optimizers/sdp/sdp.cpp
new file mode 100644
index 0000000..e4a2dbd
--- /dev/null
+++ b/src/mlpack/core/optimizers/sdp/sdp.cpp
@@ -0,0 +1,29 @@
+/**
+ * @file sdp.cpp
+ * @author Stephen Tu
+ *
+ */
+
+#include "sdp.hpp"
+
+namespace mlpack {
+namespace optimization {
+
+SDP::SDP(const size_t n,
+         const size_t numSparseConstraints,
+         const size_t numDenseConstraints) :
+    n(n),
+    sparseC(n, n),
+    denseC(n, n),
+    hasModifiedSparseObjective(false),
+    hasModifiedDenseObjective(false),
+    sparseA(numSparseConstraints),
+    sparseB(numSparseConstraints),
+    denseA(numDenseConstraints),
+    denseB(numDenseConstraints)
+{
+  denseC.zeros();
+}
+
+} // namespace optimization
+} // namespace mlpack
diff --git a/src/mlpack/core/optimizers/sdp/sdp.hpp b/src/mlpack/core/optimizers/sdp/sdp.hpp
new file mode 100644
index 0000000..83c2002
--- /dev/null
+++ b/src/mlpack/core/optimizers/sdp/sdp.hpp
@@ -0,0 +1,114 @@
+/**
+ * @file sdp.hpp
+ * @author Stephen Tu
+ *
+ */
+#ifndef __MLPACK_CORE_OPTIMIZERS_SDP_SDP_HPP
+#define __MLPACK_CORE_OPTIMIZERS_SDP_SDP_HPP
+
+#include <mlpack/core.hpp>
+
+namespace mlpack {
+namespace optimization {
+
+/**
+ * Specify an SDP in primal form
+ */
+class SDP
+{
+ public:
+
+  SDP(const size_t n,
+      const size_t numSparseConstraints,
+      const size_t numDenseConstraints);
+
+  size_t N() const { return n; }
+
+  size_t N2bar() const { return n * (n + 1) / 2; }
+
+  size_t NumSparseConstraints() const { return sparseB.n_elem; }
+
+  size_t NumDenseConstraints() const { return denseB.n_elem; }
+
+  size_t NumConstraints() const { return sparseB.n_elem + denseB.n_elem; }
+
+  //! Return the sparse objective function matrix (sparseC).
+  const arma::sp_mat& SparseC() const { return sparseC; }
+
+  //! Modify the sparse objective function matrix (sparseC).
+  arma::sp_mat& SparseC() {
+    hasModifiedSparseObjective = true;
+    return sparseC;
+  }
+
+  //! Return the dense objective function matrix (denseC).
+  const arma::mat& DenseC() const { return denseC; }
+
+  //! Modify the dense objective function matrix (denseC).
+  arma::mat& DenseC() {
+    hasModifiedDenseObjective = true;
+    return denseC;
+  }
+
+  //! Return the vector of sparse A matrices (which correspond to the sparse
+  // constraints).
+  const std::vector<arma::sp_mat>& SparseA() const { return sparseA; }
+
+  //! Modify the veector of sparse A matrices (which correspond to the sparse
+  // constraints).
+  std::vector<arma::sp_mat>& SparseA() { return sparseA; }
+
+  //! Return the vector of dense A matrices (which correspond to the dense
+  // constraints).
+  const std::vector<arma::mat>& DenseA() const { return denseA; }
+
+  //! Modify the veector of dense A matrices (which correspond to the dense
+  // constraints).
+  std::vector<arma::mat>& DenseA() { return denseA; }
+
+  //! Return the vector of sparse B values.
+  const arma::vec& SparseB() const { return sparseB; }
+  //! Modify the vector of sparse B values.
+  arma::vec& SparseB() { return sparseB; }
+
+  //! Return the vector of dense B values.
+  const arma::vec& DenseB() const { return denseB; }
+  //! Modify the vector of dense B values.
+  arma::vec& DenseB() { return denseB; }
+
+  bool HasSparseObjective() const { return hasModifiedSparseObjective; }
+
+  bool HasDenseObjective() const { return hasModifiedDenseObjective; }
+
+ private:
+
+  //! Dimension of the objective variable.
+  size_t n;
+
+  //! Sparse objective function matrix c.
+  arma::sp_mat sparseC;
+
+  //! Dense objective function matrix c.
+  arma::mat denseC;
+
+  //! If false, sparseC is zero
+  bool hasModifiedSparseObjective;
+
+  //! If false, denseC is zero
+  bool hasModifiedDenseObjective;
+
+  //! A_i for each sparse constraint.
+  std::vector<arma::sp_mat> sparseA;
+  //! b_i for each sparse constraint.
+  arma::vec sparseB;
+
+  //! A_i for each dense constraint.
+  std::vector<arma::mat> denseA;
+  //! b_i for each dense constraint.
+  arma::vec denseB;
+};
+
+} // namespace optimization
+} // namespace mlpack
+
+#endif



More information about the mlpack-git mailing list