Add more general linear operator constraint to sdp

include/ensmallen_bits/sdp/lin_alg.hpp

@@ -15,6 +15,38 @@
 namespace ens {
 namespace math {
 
+/**
+* Matrix representation of a general Linear operator representated 
+* by a function/method.
+*
+* @param input Input function.
+* @param output The matrix representation.
+*/
+template<typename ElemType, typename MatBType>
+inline void convertToMatrix(std::function<ElemType(arma::Mat<ElemType>)> input, 
+                            MatBType& output)
+{
+    const size_t n = output.n_rows;
+
+    arma::Mat<ElemType> tempInputCoordinates;
+    tempInputCoordinates.zeros(n, n);
+
+    size_t idx = 0;
+    for (size_t i = 0; i < n; i++)
+    {
+        if (i > 0) tempInputCoordinates(i - 1, n - 1) = 0;
+        for (size_t j = i; j < n; j++)
+        {
+            if (j > i) tempInputCoordinates(i, j - 1) = 0;
+            tempInputCoordinates(i, j) = 1;
+            if (i == j)
+                output(i, j) = input(tempInputCoordinates);
+            else
+                output(i, j) = output(j, i) = 0.5 * arma::datum::sqrt2 * input(tempInputCoordinates);
+        }
+    }
+}
+
 inline size_t SvecIndex(size_t i, size_t j, size_t n)
 {
   if (i > j)
@@ -109,6 +141,8 @@ inline void Smat(const MatAType& input, MatBType& output)
   }
 }
 
+
+
 /**
  * If A is a symmetric matrix, then SymKronId returns an operator Op such that
  *

include/ensmallen_bits/sdp/lrsdp.hpp

@@ -32,6 +32,21 @@ template <typename SDPType>
 class LRSDP
 {
  public:
+/**
+* Create an LRSDP to be optimized.  The solution will end up being a matrix
+* of size (rows) x (rank).  To construct each constraint and the objective
+* function, use the function SDP() in order to access the SDPType object
+* associated with this optimizer.
+*
+* @param numSparseConstraints Number of sparse constraints in the problem.
+* @param numDenseConstraints Number of dense constraints in the problem.
+* @param initialPoint Initial point of the optimization.
+* @param maxIterations Maximum number of iterations.
+*/
+LRSDP(const size_t numSparseConstraints,
+      const size_t numDenseConstraints,
+      const arma::Mat<typename SDPType::ElemType>& initialPoint,
+      const size_t maxIterations = 1000);
   /**
    * Create an LRSDP to be optimized.  The solution will end up being a matrix
    * of size (rows) x (rank).  To construct each constraint and the objective
@@ -40,11 +55,14 @@ class LRSDP
    *
    * @param numSparseConstraints Number of sparse constraints in the problem.
    * @param numDenseConstraints Number of dense constraints in the problem.
+   * @param numLinearOperatorConstraints Number of general linear operator 
+     constraints in the problem.
    * @param initialPoint Initial point of the optimization.
    * @param maxIterations Maximum number of iterations.
    */
   LRSDP(const size_t numSparseConstraints,
         const size_t numDenseConstraints,
+        const size_t numLinearOperatorConstraints,
         const arma::Mat<typename SDPType::ElemType>& initialPoint,
         const size_t maxIterations = 1000);
 
@@ -53,16 +71,15 @@ class LRSDP
    * given initial point.  Note that the SDP may be modified later by calling
    * SDP() to access the object.
    *
-   * TODO: this is currently not implemented.
    *
    * @param sdp SDP to be solved.
    * @param initialPoint Initial point of the optimization.
    * @param maxIterations Maximum number of iterations.
    *
+   */
   LRSDP(const SDPType& sdp,
         const arma::mat& initialPoint,
         const size_t maxIterations = 1000);
-   */
 
   /**
    * Optimize the LRSDP and return the final objective value.  The given

include/ensmallen_bits/sdp/lrsdp_function.hpp

@@ -48,14 +48,16 @@ class LRSDPFunction
    * constraints. Note n_cols of the initialPoint specifies the rank.
    *
    * Set the A_x, B_x, and C_x  matrices for each constraint using the A_x(),
-   * B_x(), and C_x() functions, for x in {sparse, dense}.
+   * B_x(), and C_x() functions, for x in {sparse, dense, linearOperator}.
    *
    * @param numSparseConstraints
    * @param numDenseConstraints
+   * @param numLinearOperatorConstraints
    * @param initialPoint
    */
   LRSDPFunction(const size_t numSparseConstraints,
                 const size_t numDenseConstraints,
+                const size_t numLinearOperatorConstraints,
                 const arma::Mat<typename SDPType::ElemType>& initialPoint);
 
   /**
@@ -118,6 +120,12 @@ class LRSDPFunction
     return rrt.As<typename std::remove_reference<MatType>::type>();
   }
 
+  // Cache constraint matrices of the general linear operators 
+  // stored as functions.
+  const std::vector<arma::Mat<typename SDPType::ElemType>> LoA() const {
+      return loA;
+  }
+
   //! Modify R*R^T matrix.
   template<typename MatType>
   MatType& RRT()
@@ -137,6 +145,9 @@ class LRSDPFunction
 
   //! Cache R*R^T matrix.
   Any rrt;
+
+  //! Cache linear operators in matrix form
+  std::vector<arma::Mat<typename SDPType::ElemType>> loA;
 };
 
 // Declare specializations in lrsdp_function.cpp.

include/ensmallen_bits/sdp/lrsdp_function_impl.hpp

@@ -23,22 +23,33 @@ LRSDPFunction<SDPType>::LRSDPFunction(
     const SDPType& sdp,
     const arma::Mat<typename SDPType::ElemType>& initialPoint):
     sdp(sdp),
-    initialPoint(initialPoint)
+    initialPoint(initialPoint),
+    loA(sdp.NumLinearOperatorConstraints())
 {
   if (initialPoint.n_rows < initialPoint.n_cols)
   {
     Warn << "LRSDPFunction::LRSDPFunction(): solution matrix will have "
         << "more columns than rows.  It may be more efficient to find the "
         << "transposed solution." << std::endl;
   }
+
+  for (size_t i = 0; i < sdp.NumLinearOperatorConstraints(); i++)
+  {
+      loA[i].zeros(sdp.N(), sdp.N());
+      math::convertToMatrix<typename SDPType::ElemType,
+          typename SDPType::DenseConstraintType>(sdp.LinearOperators()[i], 
+              loA[i]);
+  }
 }
 
 template<typename SDPType>
 LRSDPFunction<SDPType>::LRSDPFunction(
     const size_t numSparseConstraints,
     const size_t numDenseConstraints,
+    const size_t numLinearOperatorConstraints,
     const arma::Mat<typename SDPType::ElemType>& initialPoint):
-    sdp(initialPoint.n_rows, numSparseConstraints, numDenseConstraints),
+    sdp(initialPoint.n_rows, numSparseConstraints, numDenseConstraints, 
+        numLinearOperatorConstraints),
     initialPoint(initialPoint)
 {
   if (initialPoint.n_rows < initialPoint.n_cols)
@@ -91,11 +102,18 @@ typename MatType::elem_type LRSDPFunction<SDPType>::EvaluateConstraint(
     return accu(SDP().SparseA()[index] % rrt.As<MatType>()) -
         SDP().SparseB()[index];
   }
-  const size_t index1 = index - SDP().NumSparseConstraints();
+  else if (index < SDP().NumSparseConstraints() + SDP().NumDenseConstraints()) {
+      const size_t index1 = index - SDP().NumSparseConstraints();
+
+      // For computation optimization we will be taking R^T * A first.
+      return trace((trans(coordinates) * SDP().DenseA()[index1]) * coordinates)
+          - SDP().DenseB()[index1];
+  }
 
-  // For computation optimization we will be taking R^T * A first.
-  return trace((trans(coordinates) * SDP().DenseA()[index1]) * coordinates)
-                 - SDP().DenseB()[index1];
+  const size_t index2 = index - SDP().NumSparseConstraints() 
+      - SDP().NumDenseConstraints();
+  return SDP().LinearOperators()[index2](RRT<arma::Mat<typename SDPType::ElemType>>())
+      - SDP().LinearOperatorsB()[index2];
 }
 
 template<typename SDPType>
@@ -124,7 +142,7 @@ void UpdateRRT(LRSDPFunction<SDPType>& function,
 //! used with an LRSDPFunction.
 template <typename MatrixType, typename VecType, typename MatType>
 static inline void
-UpdateObjective(typename MatType::elem_type& objective,
+UpdateObjectiveFromMatrixConstraints(typename MatType::elem_type& objective,
                 const MatType& rrt,
                 const std::vector<MatrixType>& ais,
                 const VecType& bis,
@@ -143,12 +161,35 @@ UpdateObjective(typename MatType::elem_type& objective,
     objective += (sigma / 2.) * constraint * constraint;
   }
 }
+template <typename ElemType, typename VecType, typename MatType>
+static inline void
+UpdateObjectiveFromLinearOperatorConstraints(
+    typename MatType::elem_type& objective,
+    const MatType& rrt,
+    const std::vector<std::function<ElemType(
+           arma::Mat<ElemType>)>>& ais,
+    const VecType& bis,
+    const arma::vec& lambda,
+    const size_t lambdaOffset,
+    const double sigma)
+{
+    for (size_t i = 0; i < ais.size(); ++i)
+    {
+        // Take the trace subtracted by the b_i.
+        // Here taking R^T * A first is not recommended as we are already
+        // using pre-computed R * R^T. Taking R^T * A first will result in increase
+        // in number of computations.
+        const double constraint = ais[i](rrt) - bis[i];
+        objective -= (lambda[lambdaOffset + i] * constraint);
+        objective += (sigma / 2.) * constraint * constraint;
+    }
+}
 
 //! Utility function for calculating part of the gradient when AugLagrangian is
 //! used with an LRSDPFunction.
 template <typename MatrixType, typename VecType, typename MatType>
 static inline void
-UpdateGradient(MatType& s,
+UpdateGradientFromMatrixConstraints(MatType& s,
                const MatType& rrt,
                const std::vector<MatrixType>& ais,
                const VecType& bis,
@@ -166,6 +207,28 @@ UpdateGradient(MatType& s,
     s -= y * ais[i];
   }
 }
+template <typename ElemType, typename VecType, typename MatType>
+static inline void
+UpdateGradientFromLinearOperatorConstraints(MatType& s,
+    const MatType& rrt,
+    const std::vector<std::function<ElemType(
+        arma::Mat<ElemType>)>>&ais,
+    const VecType& bis,
+    const arma::vec& lambda,
+    const size_t lambdaOffset,
+    const double sigma,
+    const std::vector<arma::Mat<ElemType>> loA)
+{
+    for (size_t i = 0; i < ais.size(); ++i)
+    {
+        // Here taking R^T * A first is not recommended as we are already
+        // using pre-computed R * R^T. Taking R^T * A first will result in increase
+        // in number of computations.
+        const double constraint = ais[i](rrt) - bis[i];
+        const double y = lambda[lambdaOffset + i] - sigma * constraint;
+        s -= y * loA[i];
+    }
+}
 
 template<typename SDPType, typename MatType>
 static inline double
@@ -211,11 +274,18 @@ EvaluateImpl(LRSDPFunction<SDPType>& function,
       trace((trans(coordinates) * function.SDP().C()) * coordinates);
 
   // Now each constraint.
-  UpdateObjective(objective, function.template RRT<MatType>(),
+  UpdateObjectiveFromMatrixConstraints(
+      objective, function.template RRT<MatType>(),
       function.SDP().SparseA(), function.SDP().SparseB(), lambda, 0, sigma);
-  UpdateObjective(objective, function.template RRT<MatType>(),
+  UpdateObjectiveFromMatrixConstraints(
+      objective, function.template RRT<MatType>(),
       function.SDP().DenseA(), function.SDP().DenseB(), lambda,
       function.SDP().NumSparseConstraints(), sigma);
+  UpdateObjectiveFromLinearOperatorConstraints(
+      objective, function.template RRT<arma::Mat<typename SDPType::ElemType>>(),
+      function.SDP().LinearOperators(), function.SDP().LinearOperatorsB(), 
+      lambda, function.SDP().NumSparseConstraints() + 
+      function.SDP().NumDenseConstraints(), sigma);
 
   return objective;
 }
@@ -238,12 +308,17 @@ GradientImpl(const LRSDPFunction<SDPType>& function,
   const MatType& rrt = function.template RRT<MatType>();
   MatType s(function.SDP().C());
 
-  UpdateGradient(
+  UpdateGradientFromMatrixConstraints(
       s, rrt, function.SDP().SparseA(), function.SDP().SparseB(),
       lambda, 0, sigma);
-  UpdateGradient(
+  UpdateGradientFromMatrixConstraints(
       s, rrt, function.SDP().DenseA(), function.SDP().DenseB(),
       lambda, function.SDP().NumSparseConstraints(), sigma);
+  UpdateGradientFromLinearOperatorConstraints(
+      s, function.template RRT<arma::Mat<typename SDPType::ElemType>>(), 
+      function.SDP().LinearOperators(),
+      function.SDP().LinearOperatorsB(), lambda, 
+      function.SDP().NumDenseConstraints(), sigma, function.LoA());
 
   gradient = 2 * s * coordinates;
 }

include/ensmallen_bits/sdp/lrsdp_impl.hpp

@@ -17,12 +17,32 @@
 
 namespace ens {
 
+template<typename SDPType>
+LRSDP<SDPType>::LRSDP(const size_t numSparseConstraints,
+    const size_t numDenseConstraints,
+    const arma::Mat<typename SDPType::ElemType>& initialPoint,
+    const size_t maxIterations) :
+    function(numSparseConstraints, numDenseConstraints,
+        0, initialPoint),
+    maxIterations(maxIterations)
+{ }
+
 template<typename SDPType>
 LRSDP<SDPType>::LRSDP(const size_t numSparseConstraints,
                       const size_t numDenseConstraints,
+                      const size_t numLinearOperatorConstraints,
                       const arma::Mat<typename SDPType::ElemType>& initialPoint,
                       const size_t maxIterations) :
-    function(numSparseConstraints, numDenseConstraints, initialPoint),
+    function(numSparseConstraints, numDenseConstraints, 
+        numLinearOperatorConstraints, initialPoint),
+    maxIterations(maxIterations)
+{ }
+
+template<typename SDPType>
+LRSDP<SDPType>::LRSDP(const SDPType& sdp,
+    const arma::mat& initialPoint,
+    const size_t maxIterations):
+    function(sdp, initialPoint),
     maxIterations(maxIterations)
 { }
 

include/ensmallen_bits/sdp/primal_dual.hpp

@@ -74,6 +74,7 @@ class PrimalDualSolver
    * @param coordinates The initial primal SDP coordinates to optimize.
    * @param ySparse The initial ySparse to optimize.
    * @param yDense The initial yDense to optimize.
+   * @param yDense The initial yLinearOperators to optimize.
    * @param dualCoordinates The initial dual SDP coordinates to optimize.
    * @param callbacks Callback functions.
    */
@@ -82,6 +83,7 @@ class PrimalDualSolver
                                        MatType& coordinates,
                                        MatType& ySparse,
                                        MatType& yDense,
+                                       MatType& yLinearOperators,
                                        MatType& dualCoordinates,
                                        CallbackTypes&&... callbacks);