Atom domain FrankWolfe Algorithm for vector problems

src/mlpack/core/optimizers/CMakeLists.txt

@@ -6,7 +6,9 @@ set(DIRS
   fw
   gradient_descent
   lbfgs
+  line_search
   minibatch_sgd
+  proximal
   rmsprop
   sa
   sdp

src/mlpack/core/optimizers/fw/CMakeLists.txt

@@ -1,9 +1,13 @@
 set(SOURCES
+  atoms.hpp
   frank_wolfe.hpp
   frank_wolfe_impl.hpp
   constr_lpball.hpp
+  constr_structure_group.hpp
   update_classic.hpp
   update_span.hpp
+  update_linesearch.hpp
+  update_full_correction.hpp
   func_sq.hpp
   test_func_fw.hpp
 )

src/mlpack/core/optimizers/fw/atoms.hpp

@@ -0,0 +1,203 @@
+/**
+ * @file atoms.hpp
+ * @author Chenzhe Diao
+ *
+ *
+ * mlpack is free software; you may redistribute it and/or modify it under the
+ * terms of the 3-clause BSD license.  You should have received a copy of the
+ * 3-clause BSD license along with mlpack.  If not, see
+ * http://www.opensource.org/licenses/BSD-3-Clause for more information.
+ */
+#ifndef MLPACK_CORE_OPTIMIZERS_FW_ATOMS_HPP
+#define MLPACK_CORE_OPTIMIZERS_FW_ATOMS_HPP
+
+#include <mlpack/prereqs.hpp>
+#include <mlpack/core/optimizers/proximal/proximal.hpp>
+#include "func_sq.hpp"
+
+namespace mlpack {
+namespace optimization {
+
+/**
+ * Class to hold the information and operations of current atoms in the
+ * soluton space.
+ */
+class Atoms
+{
+ public:
+  Atoms(){ /* Nothing to do. */ }
+
+  /**
+   * Add atom into the solution space.
+   *
+   * @param v new atom to be added.
+   * @param c coefficient of the new atom.
+   */
+  void AddAtom(const arma::vec& v, const double c = 0)
+  {
+    if (currentAtoms.is_empty())
+    {
+      CurrentAtoms() = v;
+      CurrentCoeffs().set_size(1);
+      CurrentCoeffs().fill(c);
+    }
+    else
+    {
+      currentAtoms.insert_cols(0, v);
+      arma::vec cVec(1);
+      cVec(0) = c;
+      currentCoeffs.insert_rows(0, cVec);
+    }
+  }
+
+
+  //! Recover the solution coordinate from the coefficients of current atoms.
+  void RecoverVector(arma::mat& x)
+  {
+    x = currentAtoms * currentCoeffs;
+  }
+
+  /** 
+   * Prune the support, delete previous atoms if they don't contribute much.
+   * See Algorithm 2 of paper:
+   * @code
+   * @article{RaoShaWri:2015Forward--backward,
+   *    Author = {Rao, Nikhil and Shah, Parikshit and Wright, Stephen},
+   *    Journal = {IEEE Transactions on Signal Processing},
+   *    Number = {21},
+   *    Pages = {5798--5811},
+   *    Publisher = {IEEE},
+   *    Title = {Forward--backward greedy algorithms for atomic norm regularization},
+   *    Volume = {63},
+   *    Year = {2015}
+   * }
+   * @endcode
+   *
+   * @param F thresholding number.
+   * @param function function to be optimized.
+   */
+  void PruneSupport(const double F, FuncSq& function)
+  {
+    arma::mat atomSqTerm = function.MatrixA() * currentAtoms;
+    atomSqTerm = sum(square(atomSqTerm), 0);
+    atomSqTerm = 0.5 * atomSqTerm.t() % square(currentCoeffs);
+
+    while (true)
+    {
+      // Solve for current gradient.
+      arma::mat x;
+      RecoverVector(x);
+      arma::mat gradient(size(x));
+      function.Gradient(x, gradient);
+
+      // Find possible atom to be deleted.
+      arma::vec gap = atomSqTerm -
+          currentCoeffs % trans(gradient.t() * currentAtoms);
+      arma::uword ind;
+      gap.min(ind);
+
+      // Try deleting the atom.
+      arma::mat newAtoms = currentAtoms;
+      newAtoms.shed_col(ind);
+      // Recalculate the coefficients.
+      arma::vec newCoeffs =
+          solve(function.MatrixA() * newAtoms, function.Vectorb());
+      // Evaluate the function again.
+      double Fnew = function.Evaluate(newAtoms * newCoeffs);
+
+      if (Fnew > F)
+        // Should not delete the atom.
+        break;
+      else
+      {
+        // Delete the atom from current atoms.
+        currentAtoms = newAtoms;
+        currentCoeffs = newCoeffs;
+        atomSqTerm.shed_row(ind);
+      } // else
+    } // while
+  }
+
+
+  /**
+   * Enhance the solution in the convex hull of current atoms with atom norm
+   * constraint tau. Used in UpdateFullCorrection class for update step.
+   *
+   * Minimize the function in the atom domain defined by current atoms,
+   * where the solution still need to have atom norm (defined by current atoms)
+   * less than or equal to tau. We use projected gradient method to solve it,
+   * see the "Enhancement step" of the following paper:
+   * @code
+   * @article{RaoShaWri:2015Forward--backward,
+   *    Author = {Rao, Nikhil and Shah, Parikshit and Wright, Stephen},
+   *    Journal = {IEEE Transactions on Signal Processing},
+   *    Number = {21},
+   *    Pages = {5798--5811},
+   *    Publisher = {IEEE},
+   *    Title = {Forward--backward greedy algorithms for atomic norm regularization},
+   *    Volume = {63},
+   *    Year = {2015}
+   * }
+   * @endcode
+   *
+   * @param function function to be minimized.
+   * @param tau atom norm constraint.
+   * @param stepSize step size for projected gradient method.
+   * @param maxIteration maximum iteration number.
+   * @param tolerance tolerance for projected gradient method.
+   */
+  template<typename FunctionType>
+  void ProjectedGradientEnhancement(FunctionType& function,
+                                    double tau,
+                                    double stepSize,
+                                    size_t maxIteration = 100,
+                                    double tolerance = 1e-3)
+  {
+    arma::mat x;
+    RecoverVector(x);
+    double value = function.Evaluate(x);
+
+    Proximal proximal(tau);
+    for (size_t iter = 1; iter<maxIteration; iter++)
+    {
+      // Update currentCoeffs with gradient descent method.
+      arma::mat g;
+      function.Gradient(x, g);
+      g = currentAtoms.t() * g;
+      currentCoeffs = currentCoeffs - stepSize * g;
+
+      // Projection of currentCoeffs to satisfy the atom norm constraint.
+      proximal.ProjectToL1Ball(currentCoeffs);
+
+      RecoverVector(x);
+      double valueNew = function.Evaluate(x);
+
+      if ((value - valueNew) < tolerance)
+        break;
+
+      value = valueNew;
+    }
+  }
+
+
+  //! Get the current atom coefficients.
+  const arma::vec& CurrentCoeffs() const { return currentCoeffs; }
+  //! Modify the current atom coefficients.
+  arma::vec& CurrentCoeffs() { return currentCoeffs; }
+
+  //! Get the current atoms.
+  const arma::mat& CurrentAtoms() const { return currentAtoms; }
+  //! Modify the current atoms.
+  arma::mat& CurrentAtoms() { return currentAtoms; }
+
+ private:
+  //! Coefficients of current atoms.
+  arma::vec currentCoeffs;
+
+  //! Current atoms in the solution space.
+  arma::mat currentAtoms;
+}; // class Atoms
+}  // namespace optimization
+}  // namespace mlpack
+
+#endif

src/mlpack/core/optimizers/fw/constr_lpball.hpp

@@ -20,28 +20,37 @@ namespace optimization {
 /**
  * LinearConstrSolver for FrankWolfe algorithm. Constraint domain given in the
  * form of lp ball. That is, given \f$ v \f$, solve
- * \f[
+ * \f$
  * s:=arg\min_{s\in D} <s, v>
+ * \f$
+ * when \f$ D \f$ is a regularized lp ball. That is,
+ * \f[
+ * D = \{ x: (\sum_j|\lambda_j x_j|^p)^{1/p}\leq 1 \}.
  * \f]
- * when \f$ D \f$ is an lp ball.
+ * If \f$ \lambda \f$ is not given in the constructor, default is using all
+ * \f$ \lambda_j = 1 \f$ for all \f$ j \f$.
+ *
+ * In applications such as Orthogonal Matching Pursuit (OMP), \f$ \lambda \f$
+ * could be ideally set to the norm of the elements in the dictionary.
  *
- * For \f$ p=1 \f$: take (one) \f$ k = arg\max_j |v_j|\f$, then the solution is:
+ * For \f$ p=1 \f$: take (one) \f$ k = arg\max_j |v_j/\lambda_j|\f$, then the
+ * solution is:
  * \f[
- * s_k = -sign(v_k), \qquad s_j = 0, j\neq k.
+ * s_k = -sign(v_k)/\lambda_k, \qquad s_j = 0, \quad j\neq k.
  * \f]
  *
  * For \f$ 1<p<\infty \f$: the solution is
  * \f[
- * s_j = -sign(v_j) |v_j|^{p-1}
+ * t_j = -sign(v_j) |v_j/\lambda_j|^{q-1}, \qquad
+ * s_j = \frac{t_j}{||t||_p\cdot\lambda_j}, \quad
+ * 1/p + 1/q = 1.
  * \f]
  *
  * For \f$ p=\infty \f$: the solution is
  * \f[
- * s_j = -sign(v_j)
+ * s_j = -sign(v_j)/\lambda_j
  * \f]
  *
- * where \f$ \alpha \f$ is a parameter which specifies the step size.  \f$ i \f$
- * is chosen according to \f$ j \f$ (the iteration number).
  */
 class ConstrLpBallSolver
 {
@@ -55,6 +64,18 @@ class ConstrLpBallSolver
   ConstrLpBallSolver(const double p) : p(p)
   { /* Do nothing. */ }
 
+  /**
+   * Construct the solver of constrained problem, with regularization parameter
+   * lambda here.
+   *
+   * @param p The constraint is unit lp ball.
+   * @param lambda Regularization parameter.
+   */
+  ConstrLpBallSolver(const double p, const arma::vec lambda) :
+      p(p), regFlag(true), lambda(lambda)
+  { /* Do nothing. */ }
+
+
   /**
    * Optimizer of Linear Constrained Problem for FrankWolfe.
    *
@@ -68,20 +89,39 @@ class ConstrLpBallSolver
     {
       // l-inf ball.
       s = -sign(v);
+      if (regFlag)
+        s = s / lambda;   // element-wise division.
     }
     else if (p > 1.0)
     {
       // lp ball with 1<p<inf.
-      s = -sign(v) % pow(abs(v), p-1);
+      if (regFlag)
+        s = v / lambda;
+      else
+        s = v;
+
+      double q = 1 / (1.0 - 1.0 / p);
+      s = - sign(v) % pow(abs(s), q - 1);  // element-wise multiplication.
+      s = arma::normalise(s, p);
+
+      if (regFlag)
+        s = s / lambda;
     }
     else if (p == 1.0)
     {
       // l1 ball, also used in OMP.
-      arma::mat tmp = arma::abs(v);
+      if (regFlag)
+        s = arma::abs(v / lambda);
+      else
+        s = arma::abs(v);
+
       arma::uword k;
-      tmp.max(k);  // k is the linear index of the largest element.
-      s.zeros(v.n_rows, v.n_cols);
+      s.max(k);  // k is the linear index of the largest element.
+      s.zeros();
       s(k) = - mlpack::math::Sign(v(k));
+
+      if (regFlag)
+        s = s / lambda;
     }
     else
     {
@@ -91,10 +131,31 @@ class ConstrLpBallSolver
     return;
   }
 
+  //! Get the p-norm.
+  double P() const { return p; }
+  //! Modify the p-norm.
+  double& P() { return p;}
+
+  //! Get regularization flag.
+  bool RegFlag() const {return regFlag;}
+  //! Modify regularization flag.
+  bool& RegFlag() {return regFlag;}
+
+  //! Get the regularization parameter.
+  arma::vec Lambda() const {return lambda;}
+  //! Modify the regularization parameter.
+  arma::vec& Lambda() {return lambda;}
+
  private:
   //! lp norm, 1<=p<=inf;
   //! use std::numeric_limits<double>::infinity() for inf norm.
   double p;
+
+  //! Regularization flag.
+  bool regFlag = false;
+
+  //! Regularization parameter.
+  arma::vec lambda;
 };
 
 } // namespace optimization