[mlpack] GSoC 2014 simulated annealing optimizer

Zhihao Lou lzh1984 at gmail.com
Tue Mar 25 17:06:12 EDT 2014


Hi Ryan,

I've got the  SA done. It currently can only handle unbounded search space,
but otherwise it should work. Where should I check if the search space has
hard boundary?

My code used element-wise arma::min(), which only available in Armadillo
after version 3.930. Let me know if that's not acceptable. It is
straightforward to loop through all elements manually.

Attached is my patch, which is generated by running
svn diff > sa.diff
I'm not sure this is the correct way of doing this. Let me know if there's
more I have to do.

Zhihao Lou


On Tue, Mar 25, 2014 at 7:18 AM, Ryan Curtin <gth671b at mail.gatech.edu>wrote:

> On Mon, Mar 24, 2014 at 11:56:35PM -0500, Zhihao Lou wrote:
> > Hi Ryan,
> >
> > Thanks for the help. I'm starting working on the SA optimizer now.
> >
> > What's the project's preferred random number generator for uniform [0,1]
> > scalar? Or should I use arma::randu<vec>(1) ? Or my favorite random
> > generator?
>
> You can find some useful utilities in mlpack::math:
> http://www.mlpack.org/doxygen.php?doc=namespacemlpack_1_1math.html
>
> --
> Ryan Curtin    | "Bye-bye, goofy woman.  I enjoyed repeatedly
> ryan at ratml.org | throwing you to the ground." - Ben Jabituya
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <https://mailman.cc.gatech.edu/pipermail/mlpack/attachments/20140325/25af62aa/attachment.html>
-------------- next part --------------
Index: src/mlpack/tests/CMakeLists.txt
===================================================================
--- src/mlpack/tests/CMakeLists.txt	(revision 16374)
+++ src/mlpack/tests/CMakeLists.txt	(working copy)
@@ -35,6 +35,7 @@
   pca_test.cpp
   radical_test.cpp
   range_search_test.cpp
+  sa_test.cpp
   save_restore_utility_test.cpp
   sgd_test.cpp
   sort_policy_test.cpp
Index: src/mlpack/tests/sa_test.cpp
===================================================================
--- src/mlpack/tests/sa_test.cpp	(revision 0)
+++ src/mlpack/tests/sa_test.cpp	(revision 0)
@@ -0,0 +1,41 @@
+/*
+ * @file sa_test.cpp
+ * @auther Zhihao Lou
+ *
+ * Test file for SA (simulated annealing).
+ */
+#include <mlpack/core.hpp>
+#include <mlpack/core/optimizers/sa/sa.hpp>
+#include <mlpack/core/optimizers/lbfgs/test_functions.hpp>
+
+#include <mlpack/core/metrics/ip_metric.hpp>
+#include <mlpack/core/metrics/lmetric.hpp>
+#include <mlpack/core/metrics/mahalanobis_distance.hpp>
+
+#include <boost/test/unit_test.hpp>
+#include "old_boost_test_definitions.hpp"
+
+using namespace std;
+using namespace arma;
+using namespace mlpack;
+using namespace mlpack::optimization;
+using namespace mlpack::optimization::test;
+using namespace mlpack::metric;
+
+BOOST_AUTO_TEST_SUITE(SATest);
+
+BOOST_AUTO_TEST_CASE(GeneralizedRosenbrockTest)
+{
+  size_t dim = 50;
+  GeneralizedRosenbrockFunction f(dim);
+
+  SA<GeneralizedRosenbrockFunction> sa(f,1e-5, 1000.,1000, 100, 1e-9, 3, 20, 0.3, 0.3, 10000000);
+  arma::mat coordinates = f.GetInitialPoint();
+  double result = sa.Optimize(coordinates);
+
+  BOOST_REQUIRE_SMALL(result, 1e-6);
+  for (size_t j = 0; j < dim; ++j)
+      BOOST_REQUIRE_CLOSE(coordinates[j], (double) 1.0, 1e-2);
+}
+
+BOOST_AUTO_TEST_SUITE_END();
Index: src/mlpack/core/optimizers/sa/sa.hpp
===================================================================
--- src/mlpack/core/optimizers/sa/sa.hpp	(revision 0)
+++ src/mlpack/core/optimizers/sa/sa.hpp	(revision 0)
@@ -0,0 +1,175 @@
+/*
+ * @file sa.hpp
+ * @author Zhihao Lou
+ *
+ * Simulated Annealing (SA)
+ */
+#ifndef __MLPACK_CORE_OPTIMIZERS_SA_SA_HPP
+#define __MLPACK_CORE_OPTIMIZERS_SA_SA_HPP
+
+namespace mlpack {
+namespace optimization {
+/* 
+ * Simulated Annealing is an stochastic optimization algorithm which is able to
+ * deliver near-optimal results quickly without knowing the gradient of the
+ * function being optimized. It has unique hill climbing capability that make
+ * it less vulnerable to local minima. This implementation uses exponential
+ * cooling schedule and feedback move control.
+ *
+ * The exponential cooling schedule cools the temperature T at every step
+ * \f[
+ * T_{n+1}=(1-\lambda)T_{n}
+ * \f]
+ * where \f$ 0<\lambda<1 \f$ is the cooling speed. The smaller \f$ \lambda \f$
+ * is, the slower the cooling speed, and better the final result will be. Some
+ * literature uses \f$ \alpha=(-1\lambda) \f$ instead. In practice, \f$ \alpha \f$
+ * is very close to 1 and will be awkward to input (e.g. alpha=0.999999 vs
+ * lambda=1e-6).
+ *
+ * The algorithm keeps the temperature at initial temperature for initMove
+ * steps to get ride of the dependency of initial condition. After that, it
+ * cools every step until the system considered frozen or maxIterations is 
+ * reached.
+ *
+ * Every step SA only perturb one parameter at a time. The process that SA
+ * perturbed all parameters in a problem is called a sweep. Every moveCtrlSweep
+ * the algorithm does feedback move control to change the average move size
+ * depending on the responsiveness of each parameter. Parameter gain controls
+ * the proportion of the feedback control.
+ *
+ * The system is considered "frozen" when its score failed to change more then
+ * tolerance for consecutive maxToleranceSweep sweeps.
+ *
+ * For SA to work, a function must implement the following methods:
+ *   double Evaluate(const arma::mat& coordinates);
+ *   arma::mat& GetInitialPoint();
+ *
+ * @tparam FunctionType objective function type to be minimized.
+ */
+template<typename FunctionType>
+class SA
+{
+ public:
+  /* 
+   * Construct the SA optimizer with the given function and paramters.
+   *
+   * @param function Function to be minimized.
+   * @param lambda Cooling speed.
+   * @param initT Initial temperature.
+   * @param initMoves Iterations without changing temperature.
+   * @param moveCtrlSweep Sweeps per move control.
+   * @param tolerance Tolerance to consider system frozen.
+   * @param maxToleranceSweep Maximum sweeps below tolerance to consider system frozen.
+   * @param maxMoveCoef Maximum move size.
+   * @param initMoveCoef Initial move size.
+   * @param gain Proportional control in feedback move control.
+   * @param maxIterations Maximum number of iterations allowed (0 indicates no limit).
+   */
+  SA(FunctionType& function,
+     const double lambda = 0.001,
+     const double initT = 10000.,
+     const size_t initMoves = 1000,
+     const size_t moveCtrlSweep = 100,
+     const double tolerance = 1e-5,
+     const size_t maxToleranceSweep = 3,
+     const double maxMoveCoef = 20,
+     const double initMoveCoef = 0.3,
+     const double gain = 0.3,
+     const size_t maxIterations = 1000000);
+  /* 
+   * Optimize the given function using simulated annealing. The given starting
+   * point will be modified to store the finishing point of the algorithm, and
+   * the final objective value is returned.
+   *
+   * @param iterate Starting point (will be modified).
+   * @return Objective value of the final point.
+   */
+  double Optimize(arma::mat& iterate);
+
+  //! Get the instantiated function to be optimized.
+  const FunctionType& Function() const {return function;}
+  //! Modify the instantiated function.
+  FunctionType& Function() {return function;}
+
+  //! Get the cooling speed.
+  double Lambda() const {return lambda;}
+  //! Modify the cooling speed.
+  double& Lambda() {return lambda;}
+
+  //! Get the temperature.
+  double Temperature() const {return T;}
+  //! Modify the temperature.
+  double& Temperature() {return T;}
+
+  //! Get the initial moves.
+  size_t InitMoves() const {return initMoves;}
+  //! Modify the initial moves.
+  size_t& InitMoves() {return initMoves;}
+
+  //! Get sweeps per move control.
+  size_t MoveCtrlSweep() const {return moveCtrlSweep;}
+  //! Modify sweeps per move control.
+  size_t& MoveCtrlSweep() {return moveCtrlSweep;}
+
+  //! Get the tolerance.
+  double Tolerance() const {return tolerance;}
+  //! Modify the tolerance.
+  double& Tolerance() {return tolerance;}
+
+  //! Get the maxToleranceSweep.
+  size_t MaxToleranceSweep() const {return maxToleranceSweep;}
+  //! Modify the maxToleranceSweep.
+  size_t& MaxToleranceSweep() {return maxToleranceSweep;}
+
+  //! Get the gain.
+  double Gain() const {return gain;}
+  //! Modify the gain.
+  double& Gain() {return gain;}
+
+  //! Get the maxIterations.
+  size_t MaxIterations() const {return maxIterations;}
+  //! Modify the maxIterations.
+  size_t& MaxIterations() {return maxIterations;}
+
+  //! Get Maximum move size of each parameter
+  arma::mat MaxMove() const {return maxMove;}
+  //! Modify maximum move size of each parameter
+  arma::mat& MaxMove() {return maxMove;}
+
+  //! Get move size of each parameter
+  arma::mat MoveSize() const {return moveSize;}
+  //! Modify  move size of each parameter
+  arma::mat& MoveSize() {return moveSize;}
+
+  std::string ToString() const;
+ private:
+  FunctionType &function;
+  double lambda;
+  double T;
+  size_t initMoves;
+  size_t moveCtrlSweep;
+  double tolerance;
+  size_t maxToleranceSweep;
+  double gain;
+  size_t maxIterations;
+  arma::mat maxMove;
+  arma::mat moveSize;
+
+
+  // following variables are initialized inside Optimize
+  arma::mat accept;
+  double energy; 
+  size_t idx;
+  size_t nVars;
+  size_t sweepCounter;
+
+  void GenerateMove(arma::mat& iterate);
+  void MoveControl(size_t nMoves);
+};
+
+}; // namespace optimization
+}; // namespace mlpack
+
+#include "sa_impl.hpp"
+
+#endif
Index: src/mlpack/core/optimizers/sa/sa_impl.hpp
===================================================================
--- src/mlpack/core/optimizers/sa/sa_impl.hpp	(revision 0)
+++ src/mlpack/core/optimizers/sa/sa_impl.hpp	(revision 0)
@@ -0,0 +1,172 @@
+/*
+ * @file sa_impl.hpp
+ * @auther Zhihao Lou
+ *
+ * The implementation of the SA optimizer.
+ */
+#ifndef __MLPACK_CORE_OPTIMIZERS_SA_SA_IMPL_HPP
+#define __MLPACK_CORE_OPTIMIZERS_SA_SA_IMPL_HPP
+
+namespace mlpack {
+namespace optimization {
+
+template<typename FunctionType>
+SA<FunctionType>::SA(FunctionType& function,
+                     const double lambda,
+                     const double initT,
+                     const size_t initMoves,
+                     const size_t moveCtrlSweep,
+                     const double tolerance,
+                     const size_t maxToleranceSweep,
+                     const double maxMoveCoef,
+                     const double initMoveCoef,
+                     const double gain,
+                     const size_t maxIterations) : 
+    function(function),
+    lambda(lambda),
+    T(initT),
+    initMoves(initMoves),
+    moveCtrlSweep(moveCtrlSweep),
+    tolerance(tolerance),
+    maxToleranceSweep(maxToleranceSweep),
+    gain(gain),
+    maxIterations(maxIterations)
+{
+  const size_t rows = function.GetInitialPoint().n_rows;
+  const size_t cols = function.GetInitialPoint().n_cols;
+
+  maxMove.set_size(rows, cols);
+  maxMove.fill(maxMoveCoef);
+  moveSize.set_size(rows, cols);
+  moveSize.fill(initMoveCoef);
+  accept.zeros(rows, cols);
+}
+
+//! Optimize the function (minimize).
+template<typename FunctionType>
+double SA<FunctionType>::Optimize(arma::mat &iterate)
+{
+  const size_t rows = function.GetInitialPoint().n_rows;
+  const size_t cols = function.GetInitialPoint().n_cols;
+  size_t i;
+  size_t frozenCount = 0;
+  double alpha = 1 - lambda;
+  energy = function.Evaluate(iterate);
+  size_t old_energy = energy;
+  math::RandomSeed(std::time(NULL));
+
+  nVars = rows * cols;
+  idx = 0;
+  sweepCounter = 0;
+  accept.zeros();
+
+  // Initial Moves to get ride of dependency of initial states
+  for (i = 0; i < initMoves; ++i) 
+  {
+      GenerateMove(iterate);
+  }
+
+  // Iteration and cooling
+  for (i = 0; i != maxIterations; ++i)
+  {
+    old_energy = energy;
+    GenerateMove(iterate);
+    T *= alpha;
+    if (std::abs(energy - old_energy) < tolerance)
+      ++ frozenCount;
+    else
+      frozenCount = 0;
+    if (frozenCount >= maxToleranceSweep * nVars)
+    {
+      Log::Info << "SA: minimized within tolerance " << tolerance
+          << " for " << maxToleranceSweep << " times; terminating "
+          << "optimization." << std::endl;
+      return energy;
+    }
+  }
+  Log::Info << "SA: maximum iterations (" << maxIterations << ") reached; "
+      << "terminating optimization." << std::endl;
+  return energy;
+}
+
+template<typename FunctionType>
+void SA<FunctionType>::GenerateMove(arma::mat& iterate)
+{
+  double prevEnergy = energy;
+  double prevValue = iterate(idx);
+  // uniform [-1,1]
+  double unif = 2.0 * math::Random() - 1.0;
+  // double exponential from moveSize(idx) 
+  double move = -1.0 * moveSize(idx) * std::log(std::abs(unif));
+  if (unif < 0)
+      move = -1.0 * move;
+  iterate(idx) += move;
+  energy = function.Evaluate(iterate);
+  double xi = math::Random();
+  double delta = energy - prevEnergy;
+  double criterion = std::exp(-delta / T);
+  // Metropolis
+  if (delta <= 0. || criterion > xi)
+  {
+    accept(idx) += 1.;
+  } 
+  else
+  {
+    iterate(idx) = prevValue;
+    energy = prevEnergy;
+  }
+  ++ idx;
+  if (idx == nVars)
+  {
+    idx = 0;
+    ++ sweepCounter;
+  }
+  if (sweepCounter == moveCtrlSweep)
+  {
+    MoveControl(moveCtrlSweep);
+    sweepCounter = 0;
+  }
+}
+
+template<typename FunctionType>
+void SA<FunctionType>::MoveControl(size_t nMoves)
+{
+  /*
+   * For feedback move control and the mysterious 0.44 value, see
+   * Jimmy K.-C. Lam and Jean-Marc Delosme. An efficient simulated annealing
+   * schedule: derivation. Technical Report 8816, Yale University, 1988
+   */
+  arma::mat target;
+  target.copy_size(accept);
+  target.fill(0.44);
+  moveSize = arma::log(moveSize);
+  moveSize += gain * (accept / (double) nMoves - target);
+  moveSize = arma::exp(moveSize);
+  moveSize = arma::min(maxMove, moveSize);
+  accept.zeros();
+}
+
+template<typename FunctionType>
+std::string SA<FunctionType>::ToString() const
+{
+  std::ostringstream convert;
+  convert << "SGD [" << this << "]" << std::endl;
+  convert << "  Function:" << std::endl;
+  convert << util::Indent(function.ToString(),2);
+  convert << "  Lambda: " << lambda << std::endl;
+  convert << "  Temperature: " << T << std::endl;
+  convert << "  Initial moves: " << initMoves << std::endl;
+  convert << "  Sweeps per move control: " << moveCtrlSweep << std::endl;
+  convert << "  Tolerance: " << tolerance << std::endl;
+  convert << "  Maximum sweeps below tolerance: " << maxToleranceSweep
+      << std::endl;
+  convert << "  Move control gain: " << gain << std::endl;
+  convert << "  Maximum iterations: " << maxIterations << std::endl;
+  return convert.str();
+}
+
+}; // namespace optimization
+}; // namespace mlpack
+
+
+#endif
Index: src/mlpack/core/optimizers/sa/CMakeLists.txt
===================================================================
--- src/mlpack/core/optimizers/sa/CMakeLists.txt	(revision 0)
+++ src/mlpack/core/optimizers/sa/CMakeLists.txt	(revision 0)
@@ -0,0 +1,11 @@
+set(SOURCES
+  sa.hpp
+  sa_impl.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)
Index: src/mlpack/core/optimizers/CMakeLists.txt
===================================================================
--- src/mlpack/core/optimizers/CMakeLists.txt	(revision 16374)
+++ src/mlpack/core/optimizers/CMakeLists.txt	(working copy)
@@ -2,6 +2,7 @@
   aug_lagrangian
   lbfgs
   lrsdp
+  sa
   sgd
 )
 


More information about the mlpack mailing list