Classic FrankWolfe and OMP

src/mlpack/core/optimizers/CMakeLists.txt

@@ -3,6 +3,7 @@ set(DIRS
   ada_grad
   adam
   aug_lagrangian
+  fw
   gradient_descent
   lbfgs
   minibatch_sgd

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

@@ -0,0 +1,16 @@
+set(SOURCES
+  frank_wolfe.hpp
+  frank_wolfe_impl.hpp
+  constr_lpball.hpp
+  update_classic.hpp
+  update_span.hpp
+  func_sq.hpp
+  test_func_fw.hpp
+)
+
+set(DIR_SRCS)
+foreach(file ${SOURCES})
+  set(DIR_SRCS ${DIR_SRCS} ${CMAKE_CURRENT_SOURCE_DIR}/${file})
+endforeach()
+
+set(MLPACK_SRCS ${MLPACK_SRCS} ${DIR_SRCS} PARENT_SCOPE)

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

@@ -0,0 +1,107 @@
+/**
+ * @file constr_lpball.hpp
+ * @author Chenzhe Diao
+ *
+ * Lp ball constrained for FrankWolfe algorithm. Used as LinearConstrSolverType.
+ *
+ * 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_CONSTR_LPBALL_HPP
+#define MLPACK_CORE_OPTIMIZERS_FW_CONSTR_LPBALL_HPP
+
+#include <mlpack/prereqs.hpp>
+
+namespace mlpack {
+namespace optimization {
+
+/**
+ * LinearConstrSolver for FrankWolfe algorithm. Constraint domain given in the
+ * form of lp ball. That is, given \f$ v \f$, solve
+ * \f[
+ * s:=arg\min_{s\in D} <s, v>
+ * \f]
+ * when \f$ D \f$ is an lp ball.
+ *
+ * For \f$ p=1 \f$: take (one) \f$ k = arg\max_j |v_j|\f$, then the solution is:
+ * \f[
+ * s_k = -sign(v_k), \qquad s_j = 0, j\neq k.
+ * \f]
+ *
+ * For \f$ 1<p<\infty \f$: the solution is
+ * \f[
+ * s_j = -sign(v_j) |v_j|^{p-1}
+ * \f]
+ *
+ * For \f$ p=\infty \f$: the solution is
+ * \f[
+ * s_j = -sign(v_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
+{
+ public:
+  /**
+   * Construct the solver of constrained problem. The constrained domain should
+   * be unit lp ball for this class.
+   *
+   * @param p The constraint is unit lp ball.
+   */
+  ConstrLpBallSolver(const double p) : p(p)
+  { /* Do nothing. */ }
+
+ /**
+  * Optimizer of Linear Constrained Problem for FrankWolfe.
+  *
+  * @param v Input local gradient.
+  * @param s Output optimal solution in the constrained domain (lp ball).
+  */
+  void Optimize(const arma::mat& v,
+      arma::mat& s)
+  {
+      if (p == std::numeric_limits<double>::infinity())
+      {
+      // l-inf ball
+      s = -sign(v);
+      return;
+      }
+      else if (p > 1.0)
+      {
+      // lp ball with 1<p<inf
+      s = -sign(v) % pow(abs(v), p-1);
+      return;
+      }
+      else if (p == 1.0)
+      {
+      // l1 ball, used in OMP
+      arma::mat tmp = arma::abs(v);
+      arma::uword k = tmp.index_max();  // linear index of matrix
+      s.zeros(v.n_rows, v.n_cols);
+      s(k) = -sign_double(v(k));
+      return;
+      }
+      else
+      {
+      Log::Fatal << "Wrong norm p!" << std::endl;
+      return;
+      }
+  }
+
+ private:
+  //! lp norm, take 1<p<inf,
+  // use std::numeric_limits<double>::infinity() for inf norm.
+  double p;
+
+  //! Signum function for double.
+  double sign_double(const double x) const {return (x > 0) - (x < 0);}
+};
+
+} // namespace optimization
+} // namespace mlpack
+
+#endif

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

@@ -0,0 +1,168 @@
+/**
+ * @file frank_wolfe.hpp
+ * @author Chenzhe Diao
+ *
+ * Frank-Wolfe Algorithm.
+ *
+ * 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_FRANK_WOLFE_HPP
+#define MLPACK_CORE_OPTIMIZERS_FW_FRANK_WOLFE_HPP
+
+
+#include <mlpack/prereqs.hpp>
+#include "update_classic.hpp"
+#include "update_span.hpp"
+#include "constr_lpball.hpp"
+#include "func_sq.hpp"
+
+namespace mlpack {
+namespace optimization {
+
+/**
+ * Frank-Wolfe is a technique to minimize a continuously differentiable convex
+ * function \f$ f \f$ over a compact convex subset \f$ D \f$ of a vector space. 
+ * It is also known as conditional gradient method.
+ *
+ * To find minimum of a function using Frank-Wolfe in each iteration \f$ k \f$:
+ * 1. One optimize the linearized constrained problem, using LinearConstrSolver:
+ * \f[
+ * s_k:= arg\min_{s\in D} <s_k, \nabla f(x_k)>
+ * \f]
+ *
+ * 2. Update \f$ x \f$ using UpdateRule:
+ * \f[
+ * x_{k+1} := (1-\gamma) x_k + \gamma s_k
+ * \f]
+ * for some \f$ \gamma \in (0, 1) \f$, or use Fully-Corrective Variant:
+ * \f[
+ * x_{k+1}:= arg\min_{x\in conv(s_0, \cdots, s_k)} f(x)
+ * \f]
+ *
+ *
+ * The algorithm continues until \f$ k \f$ reaches the maximum number of iterations, 
+ * or when the duality gap is bounded by a certain tolerance \f$ \epsilon \f$. 
+ * That is,
+ *
+ * \f[
+ * g(x):= \max_{s\in D} <x-s, \nabla f(x)> \quad \leq \epsilon,
+ * \f]
+ *
+ * we also know that \f$ g(x) \geq f(x) - f(x^*) \f$, where \f$ x^* \f$ is the optimal 
+ * solution.
+ *
+ * The parameter \f$ \epsilon \f$ is specified by the tolerance parameter to the
+ * constructor.
+ *
+ * For FrankWolfe to work, FunctionType, LinearConstrSolverType and UpdateRuleType 
+ * template parameters are required.
+ * These classes must implement the following functions:
+ *
+ * FunctionType:
+ *
+ *   double Evaluate(const arma::mat& coordinates);
+ *   void Gradient(const arma::mat& coordinates,
+ *                 arma::mat& gradient);
+ *
+ * LinearConstrSolverType:
+ *
+ *   void Optimize(const arma::mat& gradient,
+ *                 arma::mat& s);
+ *
+ * UpdateRuleType:
+ *
+ *   void Update(const arma::mat& old_coords,
+ *               const arma::mat& s,
+ *               arma::mat& new_coords,
+ *               const size_t num_iter);
+ *
+ * @tparam FunctionType Objective function type to be
+ *     minimized.
+ * @tparam LinearConstrSolverType Solver for the linear constrained problem.
+ * @tparam UpdateRuleType Rule to update the solution in each iteration.
+ *
+ */
+template<
+    typename LinearConstrSolverType,
+    typename UpdateRuleType
+>
+class FrankWolfe
+{
+ public:
+  /**
+   * Construct the Frank-Wolfe optimizer with the given function and
+   * parameters. Notice that the constraint domain \f$ D \f$ is input
+   * at the initialization of linear_constr_solver, the function to be  
+   * optimized is stored in update_rule.
+   *
+   * @param linear_constr_solver Solver for linear constrained problem.
+   * @param update_rule Rule for updating solution in each iteration.
+   * @param maxIterations Maximum number of iterations allowed (0 means no
+   *     limit).
+   * @param tolerance Maximum absolute tolerance to terminate algorithm.
+   */
+  FrankWolfe(const LinearConstrSolverType linear_constr_solver,
+        const UpdateRuleType update_rule,
+        const size_t maxIterations = 100000,
+        const double tolerance = 1e-10);
+
+  /**
+   * Optimize the given function using FrankWolfe.  The given starting
+   * point will be modified to store the finishing point of the algorithm, and
+   * the final objective value is returned.
+   *
+   * @tparam function Function to be optimized.
+   * @param iterate Starting point (will be modified).
+   * @return Objective value of the final point.
+   */
+  template<typename FunctionType>
+  double Optimize(FunctionType& function, arma::mat& iterate);
+
+  //! Get the linear constrained solver.
+  LinearConstrSolverType LinearConstrSolver()
+      const { return linear_constr_solver; }
+  //! Modify the linear constrained solver.
+  LinearConstrSolverType& LinearConstrSolver() { return linear_constr_solver; }
+
+  //! Get the update rule.
+  UpdateRuleType UpdateRule() const { return update_rule; }
+  //! Modify the update rule.
+  UpdateRuleType& UpdateRule() { return update_rule; }
+
+  //! Get the maximum number of iterations (0 indicates no limit).
+  size_t MaxIterations() const { return maxIterations; }
+  //! Modify the maximum number of iterations (0 indicates no limit).
+  size_t& MaxIterations() { return maxIterations; }
+
+  //! Get the tolerance for termination.
+  double Tolerance() const { return tolerance; }
+  //! Modify the tolerance for termination.
+  double& Tolerance() { return tolerance; }
+
+ private:
+  //! The solver for constrained linear problem in first step.
+  LinearConstrSolverType linear_constr_solver;
+
+  //! The rule to update, used in the second step.
+  UpdateRuleType update_rule;
+
+  //! The maximum number of allowed iterations.
+  size_t maxIterations;
+
+  //! The tolerance for termination.
+  double tolerance;
+};
+
+//! Orthogonal Matching Pursuit
+using OMP = FrankWolfe<ConstrLpBallSolver, UpdateSpan>;
+
+} // namespace optimization
+} // namespace mlpack
+
+// Include implementation.
+#include "frank_wolfe_impl.hpp"
+
+#endif

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

@@ -0,0 +1,95 @@
+/**
+ * @file frank_wolfe_impl.hpp
+ * @author Chenzhe Diao
+ *
+ * Frank-Wolfe Algorithm.
+ *
+ * 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_FRANK_WOLFE_IMPL_HPP
+#define MLPACK_CORE_OPTIMIZERS_FW_FRANK_WOLFE_IMPL_HPP
+
+// In case it hasn't been included yet.
+#include "frank_wolfe.hpp"
+
+namespace mlpack {
+namespace optimization {
+
+template<
+    typename LinearConstrSolverType,
+    typename UpdateRuleType
+>
+FrankWolfe<LinearConstrSolverType, UpdateRuleType>::FrankWolfe(
+    const LinearConstrSolverType linear_constr_solver,
+    const UpdateRuleType update_rule,
+    const size_t maxIterations,
+    const double tolerance) :
+    linear_constr_solver(linear_constr_solver),
+    update_rule(update_rule),
+    maxIterations(maxIterations),
+    tolerance(tolerance)
+{ /* Nothing to do*/ }
+
+
+//! Optimize the function (minimize).
+template<
+    typename LinearConstrSolverType,
+    typename UpdateRuleType
+>
+template<typename FunctionType>
+double FrankWolfe<LinearConstrSolverType, UpdateRuleType>
+::Optimize(FunctionType& function, arma::mat& iterate)
+{
+    // To keep track of the function value
+    double CurrentObjective = function.Evaluate(iterate);
+    double PreviousObjective = DBL_MAX;
+
+    arma::mat gradient(iterate.n_rows, iterate.n_cols);
+    arma::mat s(iterate.n_rows, iterate.n_cols);
+    arma::mat iterate_new(iterate.n_rows, iterate.n_cols);
+    double gap = 0;
+
+    for (size_t i=1; i != maxIterations; ++i)
+    {
+    // Output current objective function
+    Log::Info << "Iteration " << i << ", objective "
+        << CurrentObjective << "." << std::endl;
+
+    // Reset counter variables.
+    PreviousObjective = CurrentObjective;
+
+    // Calculate the gradient
+    function.Gradient(iterate, gradient);
+
+    // Solve linear constrained problem, solution saved in s.
+    linear_constr_solver.Optimize(gradient, s);
+
+    // Check duality gap for return condition
+    gap = std::fabs(dot(iterate-s, gradient));
+    if (gap < tolerance)
+    {
+        Log::Info << "FrankWolfe: minimized within tolerance "
+        << tolerance << "; " << "terminating optimization." << std::endl;
+        return CurrentObjective;
+    }
+
+
+    // Update solution, save in iterate_new
+    update_rule.Update(function, iterate, s, iterate_new, i);
+
+    iterate = std::move(iterate_new);
+    CurrentObjective = function.Evaluate(iterate);
+    }
+  Log::Info << "Frank Wolfe: maximum iterations (" << maxIterations
+      << ") reached; " << "terminating optimization." << std::endl;
+  return CurrentObjective;
+}
+
+
+} // namespace optimization
+} // namespace mlpack
+
+#endif