Merge pull request #1037 from shikharbhardwaj/parallel_sgd

Implementation of parallel SGD

src/mlpack/core/optimizers/lbfgs/test_functions.cpp

@@ -196,6 +196,17 @@ void GeneralizedRosenbrockFunction::Gradient(const arma::mat& coordinates,
   gradient[i + 1] = 200 * (coordinates[i + 1] - std::pow(coordinates[i], 2));
 }
 
+void GeneralizedRosenbrockFunction::Gradient(const arma::mat& coordinates,
+                                             const size_t i,
+                                             arma::sp_mat& gradient) const
+{
+  gradient.set_size(n);
+
+  gradient[i] = 400 * (std::pow(coordinates[i], 3) - coordinates[i] *
+      coordinates[i + 1]) + 2 * (coordinates[i] - 1);
+  gradient[i + 1] = 200 * (coordinates[i + 1] - std::pow(coordinates[i], 2));
+}
+
 const arma::mat& GeneralizedRosenbrockFunction::GetInitialPoint() const
 {
   return initialPoint;

src/mlpack/core/optimizers/lbfgs/test_functions.hpp

@@ -129,6 +129,10 @@ class GeneralizedRosenbrockFunction
                 const size_t i,
                 arma::mat& gradient) const;
 
+  void Gradient(const arma::mat& coordinates,
+                const size_t i,
+                arma::sp_mat& gradient) const;
+
   const arma::mat& GetInitialPoint() const;
 
  private:

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

@@ -0,0 +1,13 @@
+set(SOURCES
+  parallel_sgd.hpp
+  parallel_sgd_impl.hpp
+  sparse_test_function.hpp
+  sparse_test_function_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)

src/mlpack/core/optimizers/parallel_sgd/decay_policies/CMakeLists.txt

@@ -0,0 +1,11 @@
+set(SOURCES
+  constant_step.hpp
+  exponential_backoff.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/parallel_sgd/decay_policies/constant_step.hpp

@@ -0,0 +1,55 @@
+/**
+ * @file constant_step.hpp
+ * @author Shikhar Bhardwaj
+ *
+ * Constant step size policy for parallel Stochastic Gradient Descent.
+ *
+ * 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_PARALLEL_SGD_CONSTANT_STEP_HPP
+#define MLPACK_CORE_OPTIMIZERS_PARALLEL_SGD_CONSTANT_STEP_HPP
+
+#include <mlpack/prereqs.hpp>
+
+namespace mlpack {
+namespace optimization {
+
+/**
+ * Implementation of the ConstantStep stepsize decay policy for parallel SGD.
+ */
+class ConstantStep
+{
+ public:
+  /**
+   * Member initialization constructor.
+   *
+   * The defaults here are not necessarily good for the given problem, so it is
+   * suggested that the values used be tailored to the task at hand.
+   *
+   * @param step The intial stepsize to use.
+   */
+  ConstantStep(const double step = 0.01) : step(step) { /* Nothing to do */ }
+
+  /**
+   * This function is called in each iteration before the gradient update.
+   *
+   * @param numEpoch The iteration number for which the stepsize is to be
+   *    calculated.
+   * @return The step size for the current iteration.
+   */
+  double StepSize(const size_t /* numEpoch */)
+  {
+    return step;
+  }
+ private:
+  //! The initial stepsize, which remains unchanged
+  double step;
+};
+
+} // namespace optimization
+} // namespace mlpack
+
+#endif

src/mlpack/core/optimizers/parallel_sgd/decay_policies/exponential_backoff.hpp

@@ -0,0 +1,92 @@
+/**
+ * @file exponential_backoff.hpp
+ * @author Shikhar Bhardwaj
+ *
+ * Exponential backoff step size decay policy for parallel Stochastic Gradient
+ * Descent.
+ *
+ * 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_PARALLEL_SGD_EXP_BACKOFF_HPP
+#define MLPACK_CORE_OPTIMIZERS_PARALLEL_SGD_EXP_BACKOFF_HPP
+
+#include <mlpack/prereqs.hpp>
+
+namespace mlpack {
+namespace optimization {
+
+/**
+ * Exponential backoff stepsize reduction policy for parallel SGD.
+ *
+ * For more information, see the following.
+ *
+ * @misc{1106.5730,
+ *   Author = {Feng Niu and Benjamin Recht and Christopher Re and Stephen J.
+ *             Wright},
+ *   Title = {HOGWILD!: A Lock-Free Approach to Parallelizing Stochastic
+ *            Gradient Descent},
+ *   Year = {2011},
+ *   Eprint = {arXiv:1106.5730},
+ * }
+ *
+ * This stepsize update scheme gives robust 1/k convergence rates to the
+ * implementation of parallel SGD.
+ */
+class ExponentialBackoff
+{
+ public:
+  /**
+   * Member initializer constructor to construct the exponential backoff policy
+   * with the required parameters.
+   *
+   * @param firstBackoffEpoch The number of updates to run before the first
+   * stepsize backoff.
+   * @param step The initial stepsize(gamma).
+   * @param beta The reduction factor. This should be a value in range (0, 1).
+   */
+  ExponentialBackoff(const size_t firstBackoffEpoch,
+                     const double step,
+                     const double beta) :
+    firstBackoffEpoch(firstBackoffEpoch),
+    cutoffEpoch(firstBackoffEpoch),
+    step(step),
+    beta(beta)
+  { /* Nothing to do. */ }
+
+  /**
+   * Get the step size for the current gradient update.
+   *
+   * @param numEpoch The iteration number of the current update.
+   * @return The stepsize for the current iteration.
+   */
+  double StepSize(const size_t numEpoch)
+  {
+    if (numEpoch >= cutoffEpoch)
+    {
+      step *= beta;
+      cutoffEpoch += firstBackoffEpoch / beta;
+    }
+    return step;
+  }
+
+ private:
+  //! The first iteration at which the stepsize should be reduced.
+  size_t firstBackoffEpoch;
+
+  //! The iteration at which the next decay will be performed.
+  size_t cutoffEpoch;
+
+  //! The initial stepsize.
+  double step;
+
+  //! The reduction factor, should be in range (0, 1).
+  double beta;
+};
+
+} // namespace optimization
+} // namespace mlpack
+
+#endif

src/mlpack/core/optimizers/parallel_sgd/parallel_sgd.hpp

@@ -0,0 +1,151 @@
+/**
+ * @file parallel_sgd.hpp
+ * @author Shikhar Bhardwaj
+ *
+ * Parallel Stochastic Gradient Descent.
+ *
+ * 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_PARALLEL_SGD_HPP
+#define MLPACK_CORE_OPTIMIZERS_PARALLEL_SGD_HPP
+
+#include <mlpack/prereqs.hpp>
+#include <mlpack/core/math/random.hpp>
+#include "decay_policies/constant_step.hpp"
+
+namespace mlpack {
+namespace optimization {
+
+/**
+ * An implementation of parallel stochastic gradient descent using the lock-free
+ * HOGWILD! approach.
+ *
+ * For more information, see the following.
+ * @misc{1106.5730,
+ *   Author = {Feng Niu and Benjamin Recht and Christopher Re and Stephen J.
+ *             Wright},
+ *   Title = {HOGWILD!: A Lock-Free Approach to Parallelizing Stochastic
+ *            Gradient Descent},
+ *   Year = {2011},
+ *   Eprint = {arXiv:1106.5730},
+ * }
+ *
+ * For Parallel SGD to work, a SparseFunctionType template parameter is
+ * required. This class must implement the following functions:
+ *
+ *   size_t NumFunctions();
+ *   double Evaluate(const arma::mat& coordinates, const size_t i);
+ *   void Gradient(const arma::mat& coordinates,
+ *                 const size_t i,
+ *                 arma::sp_mat& gradient);
+ *
+ * In these functions the parameter id refers to which individual function (or
+ * gradient) is being evaluated. In case of a data-dependent function, the id
+ * would refer to the index of the datapoint(or training example).
+ * The data is distributed uniformly among the threads made available to the
+ * program by the OpenMP runtime.
+ *
+ * The Gradient function interface is slightly changed from the
+ * DecomposableFunctionType interface, it takes in a sparse matrix as the
+ * out-param for the gradient, as ParallelSGD is only expected to be relevant in
+ * situations where the computed gradient is sparse.
+ *
+ * @tparam DecayPolicyType Step size update policy used by parallel SGD
+ *     to update the stepsize after each iteration.
+ */
+template <typename DecayPolicyType = ConstantStep>
+class ParallelSGD
+{
+ public:
+  /**
+   * Construct the parallel SGD optimizer to optimize the given function with
+   * the given parameters. One iteration means one batch of datapoints processed
+   * by each thread.
+   *
+   * The defaults here are not necessarily good for the given problem, so it is
+   * suggested that the values used be tailored to the task at hand.
+   *
+   * @param maxIterations Maximum number of iterations allowed (0 means no
+   *     limit).
+   * @param threadShareSize Number of datapoints to be processed in one
+   *     iteration by each thread.
+   * @param tolerance Maximum absolute tolerance to terminate the algorithm.
+   * @param shuffle If true, the function order is shuffled; otherwise, each
+   *     function is visited in linear order.
+   * @param decayPolicy The step size update policy to use.
+  */
+  ParallelSGD(const size_t maxIterations,
+              const size_t threadShareSize,
+              const double tolerance = 1e-5,
+              const bool shuffle = true,
+              const DecayPolicyType& decayPolicy = DecayPolicyType());
+
+  /**
+   * Optimize the given function using the parallel SGD algorithm. The given
+   * starting point will be modified to store the finishing point of the
+   * algorithm, and the value of the loss function at the final point is
+   * returned.
+   *
+   * @tparam SparseFunctionType Type of function to be optimized.
+   * @param function Function to be optimized(minimized).
+   * @param iterate Starting point(will be modified).
+   * @return Objective value at the final point.
+   */
+  template <typename SparseFunctionType>
+  double Optimize(SparseFunctionType& function, arma::mat& iterate);
+
+  //! Get the maximum number of iterations (0 indicates no limits).
+  size_t MaxIterations() const { return maxIterations; }
+  //! Modify the maximum number of iterations (0 indicates no limits).
+  size_t& MaxIterations() { return maxIterations; }
+
+  //! Get the number of datapoints to be processed in one iteration by each
+  //! thread.
+  size_t ThreadShareSize() const { return threadShareSize; }
+  //! Modify the number of datapoints to be processed in one iteration by each
+  //! thread.
+  size_t& ThreadShareSize() { return threadShareSize; }
+
+  //! Get the tolerance for termination.
+  double Tolerance() const { return tolerance; }
+  //! Modify the tolerance for termination.
+  double& Tolerance() { return tolerance; }
+
+  //! Get whether or not the individual functions are shuffled.
+  bool Shuffle() const { return shuffle; }
+  //! Modify whether or not the individual functions are shuffled.
+  bool& Shuffle() { return shuffle; }
+
+  //! Get the step size decay policy.
+  DecayPolicyType& DecayPolicy() const { return decayPolicy; }
+  //! Modify the step size decay policy.
+  DecayPolicyType& DecayPolicy() { return decayPolicy; }
+
+ private:
+  //! The maximum number of allowed iterations.
+  size_t maxIterations;
+
+  //! The number of datapoints to be processed in one iteration by each thread.
+  size_t threadShareSize;
+
+  //! The tolerance for termination.
+  double tolerance;
+
+  //! Controls whether or not the individual functions are shuffled when
+  //! iterating.
+  bool shuffle;
+
+  //! The step size decay policy.
+  DecayPolicyType decayPolicy;
+};
+
+} // namespace optimization
+} // namespace mlpack
+
+// Include implementation.
+#include "parallel_sgd_impl.hpp"
+
+#endif

src/mlpack/core/optimizers/parallel_sgd/parallel_sgd_impl.hpp

@@ -0,0 +1,130 @@
+/**
+ * @file parallel_sgd_impl.hpp
+ * @author Shikhar Bhardwaj
+ *
+ * Implementation of Parallel Stochastic Gradient Descent.
+ *
+ * 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_PARALLEL_SGD_IMPL_HPP
+#define MLPACK_CORE_OPTIMIZERS_PARALLEL_SGD_IMPL_HPP
+
+// In case it hasn't been included yet.
+#include "parallel_sgd.hpp"
+
+namespace mlpack {
+namespace optimization {
+
+template <typename DecayPolicyType>
+ParallelSGD<DecayPolicyType>::ParallelSGD(
+    const size_t maxIterations,
+    const size_t threadShareSize,
+    const double tolerance,
+    const bool shuffle,
+    const DecayPolicyType& decayPolicy) :
+    maxIterations(maxIterations),
+    threadShareSize(threadShareSize),
+    tolerance(tolerance),
+    shuffle(shuffle),
+    decayPolicy(decayPolicy)
+{ /* Nothing to do. */ }
+
+template <typename DecayPolicyType>
+template <typename SparseFunctionType>
+double ParallelSGD<DecayPolicyType>::Optimize(
+    SparseFunctionType& function,
+    arma::mat& iterate)
+{
+  double overallObjective = DBL_MAX;
+  double lastObjective;
+
+  // The order in which the functions will be visited.
+  arma::Col<size_t> visitationOrder = arma::linspace<arma::Col<size_t>>(0,
+      (function.NumFunctions() - 1), function.NumFunctions());
+
+  // Iterate till the objective is within tolerance or the maximum number of
+  // allowed iterations is reached. If maxIterations is 0, this will iterate
+  // till convergence.
+  for (size_t i = 1; i != maxIterations; ++i)
+  {
+    // Calculate the overall objective.
+    lastObjective = overallObjective;
+    overallObjective = 0;
+
+    for (size_t j = 0; j < function.NumFunctions(); ++j)
+    {
+      overallObjective += function.Evaluate(iterate, j);
+    }
+
+    // Output current objective function.
+    Log::Info << "Parallel SGD: iteration " << i << ", objective "
+      << overallObjective << "." << std::endl;
+
+    if (std::isnan(overallObjective) || std::isinf(overallObjective))
+    {
+      Log::Warn << "Parallel SGD: converged to " << overallObjective
+        << "; terminating with failure. Try a smaller step size?"
+        << std::endl;
+      return overallObjective;
+    }
+
+    if (std::abs(lastObjective - overallObjective) < tolerance)
+    {
+      Log::Info << "SGD: minimized within tolerance " << tolerance << "; "
+        << "terminating optimization." << std::endl;
+      return overallObjective;
+    }
+
+    // Get the stepsize for this iteration
+    double stepSize = decayPolicy.StepSize(i);
+
+    // Shuffle for uniform sampling of functions by each thread.
+
+    if (shuffle) // Determine order of visitation.
+      std::shuffle(visitationOrder.begin(), visitationOrder.end(),
+          mlpack::math::randGen);
+
+    #pragma omp parallel
+    {
+      // Each processor gets a subset of the instances.
+      // Each subset is of size threadShareSize.
+      size_t threadId = omp_get_thread_num();
+
+      for (size_t j = threadId * threadShareSize;
+          j < (threadId + 1) * threadShareSize && j < visitationOrder.n_elem;
+          ++j)
+      {
+        // Each instance affects only some components of the decision variable.
+        // So the gradient is sparse.
+        arma::sp_mat gradient;
+
+        // Evaluate the sparse gradient.
+        function.Gradient(iterate, visitationOrder[j], gradient);
+
+        // Update the decision variable with non-zero components of the
+        // gradient.
+        for (size_t i = 0; i < gradient.n_cols; ++i)
+        {
+          // Iterate over the non-zero elements.
+          for (arma::sp_mat::iterator cur = gradient.begin_col(i);
+              cur != gradient.end_col(i); ++cur)
+          {
+            #pragma omp atomic
+            iterate(cur.row(), i) -= stepSize * (*cur);
+          }
+        }
+      }
+    }
+  }
+  Log::Info << "\n Parallel SGD terminated with objective : "
+    << overallObjective << std::endl;
+  return overallObjective;
+}
+
+} // namespace optimization
+} // namespace mlpack
+
+#endif

src/mlpack/core/optimizers/parallel_sgd/sparse_test_function.hpp

@@ -0,0 +1,63 @@
+/**
+ * @file sparse_test_function.hpp
+ * @author Shikhar Bhardwaj
+ *
+ * Sparse test function for Parallel SGD.
+ *
+ * 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_PARALLEL_SGD_SPARSE_TEST_FUNCTION_HPP
+#define MLPACK_CORE_OPTIMIZERS_PARALLEL_SGD_SPARSE_TEST_FUNCTION_HPP
+
+#include <mlpack/prereqs.hpp>
+
+namespace mlpack {
+namespace optimization {
+namespace test {
+
+// A simple test function. Each dimension has a parabola with a
+// distinct minimum. Each update is guaranteed to be sparse(only a single
+// dimension is updated in the decision variable by each thread). At the end of
+// a reasonable number of iterations, each value in the decision variable should
+// be at the vertex of the parabola in that dimension.
+class SparseTestFunction
+{
+ public:
+  //! Set members in the default constructor.
+  SparseTestFunction();
+
+  //! Return 4 (the number of functions).
+  size_t NumFunctions() const { return 4; }
+
+  //! Get the starting point.
+  arma::mat GetInitialPoint() const { return arma::mat("0; 0; 0; 0;"); }
+
+  //! Evaluate a function.
+  double Evaluate(const arma::mat& coordinates, const size_t i) const;
+
+  //! Evaluate the gradient of a function.
+  void Gradient(const arma::mat& coordinates,
+                const size_t i,
+                arma::sp_mat& gradient) const;
+ private:
+  // Each quadratic polynomial is monic. The intercept and coefficient of the
+  // first order term is stored.
+
+  //! The vector storing the intercepts
+  arma::vec intercepts;
+
+  //! The vector having coefficients of the first order term
+  arma::vec bi;
+};
+
+} // namespace test
+} // namespace optimization
+} // namespace mlpack
+
+// Include implementation
+#include "sparse_test_function_impl.hpp"
+
+#endif

src/mlpack/core/optimizers/parallel_sgd/sparse_test_function_impl.hpp

@@ -0,0 +1,50 @@
+/**
+ * @file sparse_test_function.hpp
+ * @author Shikhar Bhardwaj
+ *
+ * Sparse test function for Parallel SGD.
+ *
+ * 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_PARALLEL_SGD_SPARSE_TEST_FUNCTION_IMPL_HPP
+#define MLPACK_CORE_OPTIMIZERS_PARALLEL_SGD_SPARSE_TEST_FUNCTION_IMPL_HPP
+
+// In case it hasn't been included yet.
+#include "sparse_test_function.hpp"
+
+namespace mlpack {
+namespace optimization {
+namespace test {
+
+//! The default constructor sets the members.
+SparseTestFunction::SparseTestFunction()
+{
+  intercepts = arma::vec("20 12 15 100");
+  bi = arma::vec("-4 -2 -3 -8");
+}
+
+//! Evaluate a function.
+double SparseTestFunction::Evaluate(
+    const arma::mat& coordinates, const size_t i) const
+{
+  return coordinates[i] * coordinates[i] + bi[i] * coordinates[i] +
+    intercepts[i];
+}
+
+  //! Evaluate the gradient of a function.
+void SparseTestFunction::Gradient(const arma::mat& coordinates,
+    const size_t i,
+    arma::sp_mat& gradient) const
+{
+  gradient = arma::sp_mat(coordinates.n_rows, 1);
+  gradient[i] = 2 * coordinates[i] + bi[i];
+}
+
+} // namespace test
+} // namespace optimization
+} // namespace mlpack
+
+#endif

src/mlpack/methods/logistic_regression/logistic_regression_function.hpp

@@ -94,15 +94,18 @@ class LogisticRegressionFunction
    * Evaluate the gradient of the logistic regression log-likelihood function
    * with the given parameters, and with respect to only one point in the
    * dataset.  This is useful for optimizers such as SGD, which require a
-   * separable objective function.
+   * separable objective function.  The type of the gradient parameter is a
+   * template argument to allow the computation of a sparse gradient.
    *
+   * @tparam GradType The type of the gradient out-param.
    * @param parameters Vector of logistic regression parameters.
    * @param i Index of points to use for objective function gradient evaluation.
    * @param gradient Vector to output gradient into.
    */
+  template <typename GradType>
   void Gradient(const arma::mat& parameters,
                 const size_t i,
-                arma::mat& gradient) const;
+                GradType& gradient) const;
 
   //! Return the initial point for the optimization.
   const arma::mat& GetInitialPoint() const { return initialPoint; }

src/mlpack/methods/logistic_regression/logistic_regression_function_impl.hpp

@@ -149,14 +149,15 @@ void LogisticRegressionFunction<MatType>::Gradient(
  * function with respect to individual points.  This is useful for optimizers
  * that use a separable objective function, such as SGD.
  */
-template<typename MatType>
+template <typename MatType>
+template <typename GradType>
 void LogisticRegressionFunction<MatType>::Gradient(
     const arma::mat& parameters,
     const size_t i,
-    arma::mat& gradient) const
+    GradType& gradient) const
 {
   // Calculate the regularization term.
-  arma::mat regularization;
+  GradType regularization;
   regularization = lambda * parameters.col(0).subvec(1, parameters.n_elem - 1)
       / predictors.n_cols;
 

src/mlpack/methods/nca/nca_softmax_error_function.hpp

@@ -14,6 +14,7 @@
 #define MLPACK_METHODS_NCA_NCA_SOFTMAX_ERROR_FUNCTION_HPP
 
 #include <mlpack/prereqs.hpp>
+#include <mlpack/core/metrics/lmetric.hpp>
 
 namespace mlpack {
 namespace nca {
@@ -92,15 +93,18 @@ class SoftmaxErrorFunction
    * matrix on only one point of the dataset.  This is the separable
    * implementation, where the objective function is decomposed into the sum of
    * many objective functions, and here, only one of those constituent objective
-   * functions is returned.
+   * functions is returned. The type of the gradient parameter is a template
+   * argument to allow the computation of a sparse gradient.
    *
+   * @tparam GradType The type of the gradient out-param.
    * @param covariance Covariance matrix of Mahalanobis distance.
    * @param i Index of point to use for objective function.
    * @param gradient Matrix to store the calculated gradient in.
    */
+  template <typename GradType>
   void Gradient(const arma::mat& covariance,
                 const size_t i,
-                arma::mat& gradient);
+                GradType& gradient);
 
   /**
    * Get the initial point.

src/mlpack/methods/nca/nca_softmax_error_function_impl.hpp

@@ -132,10 +132,11 @@ void SoftmaxErrorFunction<MetricType>::Gradient(const arma::mat& coordinates,
 }
 
 //! The separable implementation.
-template<typename MetricType>
+template <typename MetricType>
+template <typename GradType>
 void SoftmaxErrorFunction<MetricType>::Gradient(const arma::mat& coordinates,
                                                 const size_t i,
-                                                arma::mat& gradient)
+                                                GradType& gradient)
 {
   // We will need to calculate p_i before this evaluation is done, so these two
   // variables will hold the information necessary for that.
@@ -144,8 +145,8 @@ void SoftmaxErrorFunction<MetricType>::Gradient(const arma::mat& coordinates,
 
   // The gradient involves two matrix terms which are eventually combined into
   // one.
-  arma::mat firstTerm;
-  arma::mat secondTerm;
+  GradType firstTerm;
+  GradType secondTerm;
 
   firstTerm.zeros(coordinates.n_rows, coordinates.n_cols);
   secondTerm.zeros(coordinates.n_rows, coordinates.n_cols);
@@ -166,7 +167,7 @@ void SoftmaxErrorFunction<MetricType>::Gradient(const arma::mat& coordinates,
     // If the points are in the same class, we must add to the second term of
     // the gradient as well as the numerator of p_i.  We will divide by the
     // denominator of p_ik later.  For x_ik we are not using stretched points.
-    arma::vec x_ik = dataset.col(i) - dataset.col(k);
+    GradType x_ik = dataset.col(i) - dataset.col(k);
     if (labels[i] == labels[k])
     {
       numerator += eval;

src/mlpack/methods/regularized_svd/CMakeLists.txt

@@ -4,7 +4,7 @@ set(SOURCES
   regularized_svd.hpp
   regularized_svd_impl.hpp
   regularized_svd_function.hpp
-  regularized_svd_function.cpp
+  regularized_svd_function_impl.hpp
 )
 
 # Add directory name to sources.

src/mlpack/methods/regularized_svd/regularized_svd_function.hpp

@@ -15,10 +15,19 @@
 
 #include <mlpack/prereqs.hpp>
 #include <mlpack/core/optimizers/sgd/sgd.hpp>
+#include <mlpack/core/optimizers/parallel_sgd/parallel_sgd.hpp>
+#include <mlpack/core/optimizers/parallel_sgd/decay_policies/exponential_backoff.hpp>
 
 namespace mlpack {
 namespace svd {
 
+/**
+ * The data is stored in a matrix of type MatType, so that this class can be
+ * used with both dense and sparse matrix types.
+ *
+ * @tparam MatType The matrix type of the dataset.
+ */
+template <typename MatType = arma::mat>
 class RegularizedSVDFunction
 {
  public:
@@ -31,7 +40,7 @@ class RegularizedSVDFunction
    * @param rank Rank used for matrix factorization.
    * @param lambda Regularization parameter used for optimization.
    */
-  RegularizedSVDFunction(const arma::mat& data,
+  RegularizedSVDFunction(const MatType& data,
                          const size_t rank,
                          const double lambda);
 
@@ -62,6 +71,22 @@ class RegularizedSVDFunction
   void Gradient(const arma::mat& parameters,
                 arma::mat& gradient) const;
 
+  /**
+   * Evaluates the gradient of the cost function over one training example.
+   * This function is useful for optimizers like SGD. The type of the gradient
+   * parameter is a template argument to allow the computation of a sparse
+   * gradient.
+   *
+   * @tparam GradType The type of the gradient out-param.
+   * @param parameters Parameters(user/item matrices) of the decomposition.
+   * @param id The index of the training example.
+   * @param gradient Calculated gradient for the parameters.
+   */
+  template <typename GradType>
+  void Gradient(const arma::mat& parameters,
+                size_t id,
+                GradType& gradient) const;
+
   //! Return the initial point for the optimization.
   const arma::mat& GetInitialPoint() const { return initialPoint; }
 
@@ -85,7 +110,7 @@ class RegularizedSVDFunction
 
  private:
   //! Rating data.
-  const arma::mat& data;
+  const MatType& data;
   //! Initial parameter point.
   arma::mat initialPoint;
   //! Rank used for matrix factorization.
@@ -105,17 +130,26 @@ namespace mlpack {
 namespace optimization {
 
   /**
-   * Template specialization for SGD optimizer. Used because the gradient
-   * affects only a small number of parameters per example, and thus the normal
-   * abstraction does not work as fast as we might like it to.
+   * Template specialization for the SGD and parallel SGD optimizer. Used
+   * because the gradient affects only a small number of parameters per example,
+   * and thus the normal abstraction does not work as fast as we might like it
+   * to.
    */
-  template<>
-  template<>
-  double StandardSGD::Optimize(
-      mlpack::svd::RegularizedSVDFunction& function,
+  template <>
+  template <>
+  inline double StandardSGD::Optimize(
+      mlpack::svd::RegularizedSVDFunction<arma::mat>& function,
+      arma::mat& parameters);
+
+  template <>
+  template <>
+  inline double ParallelSGD<ExponentialBackoff>::Optimize(
+      mlpack::svd::RegularizedSVDFunction<arma::mat>& function,
       arma::mat& parameters);
 
 } // namespace optimization
 } // namespace mlpack
 
+#include "regularized_svd_function_impl.hpp"
+
 #endif

src/mlpack/methods/regularized_svd/regularized_svd_function_impl.hpp

@@ -10,15 +10,19 @@
  * http://www.opensource.org/licenses/BSD-3-Clause for more information.
  */
 
+
+#ifndef MLPACK_METHODS_REGULARIZED_SVD_REGULARIZED_FUNCTION_SVD_IMPL_HPP
+#define MLPACK_METHODS_REGULARIZED_SVD_REGULARIZED_FUNCTION_SVD_IMPL_HPP
+
 #include "regularized_svd_function.hpp"
-#include <mlpack/core/optimizers/sgd/sgd.hpp>
 
 namespace mlpack {
 namespace svd {
 
-RegularizedSVDFunction::RegularizedSVDFunction(const arma::mat& data,
-                                               const size_t rank,
-                                               const double lambda) :
+template <typename MatType>
+RegularizedSVDFunction<MatType>::RegularizedSVDFunction(const MatType& data,
+                                                        const size_t rank,
+                                                        const double lambda) :
     data(data),
     rank(rank),
     lambda(lambda)
@@ -31,7 +35,9 @@ RegularizedSVDFunction::RegularizedSVDFunction(const arma::mat& data,
   initialPoint.randu(rank, numUsers + numItems);
 }
 
-double RegularizedSVDFunction::Evaluate(const arma::mat& parameters) const
+template <typename MatType>
+double RegularizedSVDFunction<MatType>::Evaluate(const arma::mat& parameters)
+const
 {
   // The cost for the optimization is as follows:
   //          f(u, v) = sum((rating(i, j) - u(i).t() * v(j))^2)
@@ -66,8 +72,9 @@ double RegularizedSVDFunction::Evaluate(const arma::mat& parameters) const
   return cost;
 }
 
-double RegularizedSVDFunction::Evaluate(const arma::mat& parameters,
-                                        const size_t i) const
+template <typename MatType>
+double RegularizedSVDFunction<MatType>::Evaluate(const arma::mat& parameters,
+                                                 const size_t i) const
 {
   // Indices for accessing the the correct parameter columns.
   const size_t user = data(0, i);
@@ -88,8 +95,9 @@ double RegularizedSVDFunction::Evaluate(const arma::mat& parameters,
   return (ratingErrorSquared + regularizationError);
 }
 
-void RegularizedSVDFunction::Gradient(const arma::mat& parameters,
-                                      arma::mat& gradient) const
+template <typename MatType>
+void RegularizedSVDFunction<MatType>::Gradient(const arma::mat& parameters,
+                                               arma::mat& gradient) const
 {
   // For an example with rating corresponding to user 'i' and item 'j', the
   // gradients for the parameters is as follows:
@@ -122,6 +130,30 @@ void RegularizedSVDFunction::Gradient(const arma::mat& parameters,
   }
 }
 
+template <typename MatType>
+template <typename GradType>
+void RegularizedSVDFunction<MatType>::Gradient(const arma::mat& parameters,
+                                               size_t id,
+                                               GradType &gradient) const
+{
+  gradient.zeros(rank, numUsers + numItems);
+
+  const size_t user = data(0, id);
+  const size_t item = data(1, id) + numUsers;
+
+  // Prediction error for the example.
+  const double rating = data(2, id);
+  double ratingError = rating - arma::dot(parameters.col(user),
+                                          parameters.col(item));
+
+  // Gradient is non-zero only for the parameter columns corresponding to the
+  // example.
+  gradient.col(user) += 2 * (lambda * parameters.col(user) -
+                             ratingError * parameters.col(item));
+  gradient.col(item) += 2 * (lambda * parameters.col(item) -
+                             ratingError * parameters.col(user));
+}
+
 } // namespace svd
 } // namespace mlpack
 
@@ -129,10 +161,10 @@ void RegularizedSVDFunction::Gradient(const arma::mat& parameters,
 namespace mlpack {
 namespace optimization {
 
-template<>
-template<>
+template <>
+template <>
 double StandardSGD::Optimize(
-    mlpack::svd::RegularizedSVDFunction& function,
+    mlpack::svd::RegularizedSVDFunction<arma::mat>& function,
     arma::mat& parameters)
 {
   // Find the number of functions to use.
@@ -186,5 +218,110 @@ double StandardSGD::Optimize(
   return overallObjective;
 }
 
+
+template <>
+template <>
+inline double ParallelSGD<ExponentialBackoff>::Optimize(
+    mlpack::svd::RegularizedSVDFunction<arma::mat>& function,
+    arma::mat& iterate)
+{
+  double overallObjective = DBL_MAX;
+  double lastObjective;
+
+  // The order in which the functions will be visited.
+  arma::Col<size_t> visitationOrder = arma::linspace<arma::Col<size_t>>(0,
+      (function.NumFunctions() - 1), function.NumFunctions());
+
+  const arma::mat data = function.Dataset();
+
+  // Iterate till the objective is within tolerance or the maximum number of
+  // allowed iterations is reached. If maxIterations is 0, this will iterate
+  // till convergence.
+  for (size_t i = 1; i != maxIterations; ++i)
+  {
+    // Calculate the overall objective.
+    lastObjective = overallObjective;
+    overallObjective = 0;
+
+    #pragma omp parallel for reduction(+:overallObjective)
+    for (omp_size_t j = 0; j < (omp_size_t) function.NumFunctions(); ++j)
+    {
+      overallObjective += function.Evaluate(iterate, j);
+    }
+
+    // Output current objective function.
+    Log::Info << "Parallel SGD: iteration " << i << ", objective "
+      << overallObjective << "." << std::endl;
+
+    if (std::isnan(overallObjective) || std::isinf(overallObjective))
+    {
+      Log::Warn << "Parallel SGD: converged to " << overallObjective
+        << "; terminating with failure. Try a smaller step size?"
+        << std::endl;
+      return overallObjective;
+    }
+
+    if (std::abs(lastObjective - overallObjective) < tolerance)
+    {
+      Log::Info << "SGD: minimized within tolerance " << tolerance << "; "
+        << "terminating optimization." << std::endl;
+      return overallObjective;
+    }
+
+    // Get the stepsize for this iteration
+    double stepSize = decayPolicy.StepSize(i);
+
+    if (shuffle) // Determine order of visitation.
+      std::shuffle(visitationOrder.begin(), visitationOrder.end(),
+          mlpack::math::randGen);
+
+    #pragma omp parallel
+    {
+      // Each processor gets a subset of the instances.
+      // Each subset is of size threadShareSize.
+      size_t threadId = omp_get_thread_num();
+
+      for (size_t j = threadId * threadShareSize;
+          j < (threadId + 1) * threadShareSize && j < visitationOrder.n_elem;
+          ++j)
+      {
+        const size_t numUsers = function.NumUsers();
+
+        // Indices for accessing the the correct parameter columns.
+        const size_t user = data(0, visitationOrder[j]);
+        const size_t item = data(1, visitationOrder[j]) + numUsers;
+
+        // Prediction error for the example.
+        const double rating = data(2, visitationOrder[j]);
+        double ratingError = rating - arma::dot(iterate.col(user),
+            iterate.col(item));
+
+        double lambda = function.Lambda();
+
+        arma::mat userUpdate = stepSize * (lambda * iterate.col(user) -
+            ratingError * iterate.col(item));
+        arma::mat itemUpdate = stepSize * (lambda * iterate.col(item) -
+            ratingError * iterate.col(user));
+
+        // Gradient is non-zero only for the parameter columns corresponding to
+        // the example.
+        for (size_t i = 0; i < iterate.n_rows; ++i)
+        {
+          #pragma omp atomic
+          iterate(i, user) -= userUpdate(i);
+          #pragma omp atomic
+          iterate(i, item) -= itemUpdate(i);
+        }
+      }
+    }
+  }
+  Log::Info << "\n Parallel SGD terminated with objective : "
+    << overallObjective << std::endl;
+
+  return overallObjective;
+}
+
 } // namespace optimization
 } // namespace mlpack
+
+#endif

src/mlpack/methods/regularized_svd/regularized_svd_impl.hpp

@@ -34,7 +34,7 @@ void RegularizedSVD<OptimizerType>::Apply(const arma::mat& data,
                                           arma::mat& v)
 {
   // Make the optimizer object using a RegularizedSVDFunction object.
-  RegularizedSVDFunction rSVDFunc(data, rank, lambda);
+  RegularizedSVDFunction<arma::mat> rSVDFunc(data, rank, lambda);
   mlpack::optimization::StandardSGD optimizer(alpha, iterations * data.n_cols);
 
   // Get optimized parameters.

src/mlpack/methods/sparse_svm/CMakeLists.txt

@@ -0,0 +1,16 @@
+# Define the files we need to compile
+# Anything not in this list will not be compiled into the output library
+# Do not include test programs here
+set(SOURCES
+  sparse_svm_function.hpp
+  sparse_svm_function_impl.hpp
+)
+
+# add directory name to sources
+set(DIR_SRCS)
+foreach(file ${SOURCES})
+    set(DIR_SRCS ${DIR_SRCS} ${CMAKE_CURRENT_SOURCE_DIR}/${file})
+endforeach()
+# append sources (with directory name) to list of all mlpack sources (used at
+# the parent scope)
+set(MLPACK_SRCS ${MLPACK_SRCS} ${DIR_SRCS} PARENT_SCOPE)

src/mlpack/methods/sparse_svm/sparse_svm_function.hpp

@@ -0,0 +1,78 @@
+/**
+ * @file sparse_svm_function.hpp
+ * @author Shikhar Bhardwaj
+ *
+ * Implementation of the hinge loss function for training a sparse SVM with the
+ * parallel SGD 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_METHODS_SPARSE_SVM_SPARSE_SVM_FUNCTION_HPP
+#define MLPACK_METHODS_SPARSE_SVM_SPARSE_SVM_FUNCTION_HPP
+
+#include <mlpack/prereqs.hpp>
+
+class SparseSVMFunction{
+ public:
+  //! Nothing to do for the default constructor.
+  SparseSVMFunction() {}
+
+  //! Member initialization constructor.
+  SparseSVMFunction(const arma::sp_mat& dataset, const arma::vec& labels);
+
+  /**
+   * Evaluate the hinge loss function on the specified datapoint.
+   *
+   * @param parameters The parameters of the SVM.
+   * @param id Index of the datapoint to use for function evaluation.
+   * @return The value of the loss function at the given parameters.
+   */
+  double Evaluate(const arma::mat& parameters, size_t id);
+
+  /**
+   * Evaluate the gradient the gradient of the hinge loss function, following
+   * the SparseFunctionType requirements on the Gradient function.
+   *
+   * @param parameters The parameters of the SVM.
+   * @param id Index of the datapoint to use for the gradient evaluation.
+   * @param gradient Sparse matrix to output the gradient into.
+   */
+  template <typename GradType>
+  void Gradient(const arma::mat& parameters, size_t id, GradType& gradient);
+
+  //! Return the initial point for the optimization.
+  const arma::mat& InitialPoint() const { return initialPoint; }
+  //! Modify the initial point for the optimization.
+  arma::mat& InitialPoint() { return initialPoint; }
+
+  //! Get the dataset.
+  const arma::sp_mat& Dataset() const { return dataset; }
+  //! Modify the dataset.
+  arma::sp_mat& Dataset() { return dataset; }
+
+  //! Get the labels.
+  const arma::vec& Labels() const { return labels; }
+  //! Modify the labels.
+  arma::vec& Labels() { return labels; }
+
+  //! Return the number of functions.
+  size_t NumFunctions();
+
+ private:
+  //! The initial point, from which to start the optimization.
+  arma::mat initialPoint;
+
+  //! The datapoints for training.
+  arma::sp_mat dataset;
+
+  //! The labels, y_i.
+  arma::vec labels;
+};
+
+// Include implementation
+#include "sparse_svm_function_impl.hpp"
+
+#endif // MLPACK_METHODS_SPARSE_SVM_SPARSE_SVM_FUNCTION_HPP

src/mlpack/methods/sparse_svm/sparse_svm_function_impl.hpp

@@ -0,0 +1,47 @@
+/**
+ * @file sparse_svm_function_impl.hpp
+ * @author Shikhar Bhardwaj
+ *
+ * Implementation of the hinge loss function for training a sparse SVM with the
+ * parallel SGD 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_METHODS_SPARSE_SVM_SPARSE_SVM_FUNCTION_IMPL_HPP
+#define MLPACK_METHODS_SPARSE_SVM_SPARSE_SVM_FUNCTION_IMPL_HPP
+
+// In case it hasn't been included yet.
+#include "sparse_svm_function.hpp"
+
+SparseSVMFunction::SparseSVMFunction(
+    const arma::sp_mat& dataset, const arma::vec& labels) :
+  dataset(dataset), labels(labels)
+{ /* Nothing to do */ }
+
+double SparseSVMFunction::Evaluate(const arma::mat& parameters, size_t id)
+{
+  // The hinge loss function.
+  return std::max(0.0, 1 - labels(id) * arma::dot(dataset.col(id), parameters));
+}
+
+template <typename GradType>
+void SparseSVMFunction::Gradient(
+    const arma::mat& parameters, size_t id, GradType& gradient)
+{
+  // Evaluate the gradient of the hinge loss function.
+  double dot = 1 - labels(id) * arma::dot(parameters, dataset.col(id));
+  gradient = (dot < 0) ? GradType(parameters.n_rows, 1) :
+    (-1 * GradType(dataset.col(id) * labels(id)));
+}
+
+size_t SparseSVMFunction::NumFunctions()
+{
+  // The number of points in the dataset is the number of functions, as this
+  // is a data dependent function.
+  return dataset.n_cols;
+}
+
+#endif // MLPACK_METHODS_SPARSE_SVM_SPARSE_SVM_FUNCTION_IMPL_HPP

src/mlpack/tests/CMakeLists.txt

@@ -68,6 +68,7 @@ add_executable(mlpack_test
   nmf_test.cpp
   nystroem_method_test.cpp
   octree_test.cpp
+  parallel_sgd_test.cpp
   pca_test.cpp
   perceptron_test.cpp
   prefixedoutstream_test.cpp

src/mlpack/tests/parallel_sgd_test.cpp

@@ -0,0 +1,119 @@
+/**
+ * @file parallel_sgd_test.cpp
+ * @author Shikhar Bhardwaj
+ *
+ * Test file for Parallel SGD.
+ *
+ * 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.
+ */
+#include <mlpack/core.hpp>
+#include <mlpack/core/optimizers/parallel_sgd/decay_policies/constant_step.hpp>
+#include <mlpack/core/optimizers/parallel_sgd/decay_policies/exponential_backoff.hpp>
+#include <mlpack/core/optimizers/parallel_sgd/sparse_test_function.hpp>
+#include <mlpack/core/optimizers/lbfgs/test_functions.hpp>
+// We need some thorough testing
+#define private public
+#include <mlpack/core/optimizers/parallel_sgd/parallel_sgd.hpp>
+#undef private
+
+#include <boost/test/unit_test.hpp>
+#include "test_tools.hpp"
+
+using namespace std;
+using namespace arma;
+using namespace mlpack;
+using namespace mlpack::optimization;
+using namespace mlpack::optimization::test;
+
+BOOST_AUTO_TEST_SUITE(ParallelSGDTest);
+
+/**
+ * Test the correctness of the Parallel SGD implementation using a specified
+ * sparse test function, with guaranteed disjoint updates between different
+ * threads.
+ */
+BOOST_AUTO_TEST_CASE(SimpleParallelSGDTest)
+{
+  SparseTestFunction f;
+
+  ConstantStep decayPolicy(0.4);
+
+  // The batch size for this test should be chosen according to the threads
+  // available on the system. If the update does not touch each datapoint, the
+  // test will fail.
+
+  size_t threadsAvailable = omp_get_max_threads();
+
+  for (size_t i = threadsAvailable; i > 0; --i)
+  {
+    omp_set_num_threads(i);
+
+    size_t batchSize = std::ceil((float) f.NumFunctions() / i);
+
+    ParallelSGD<ConstantStep> s(10000, batchSize, 1e-5, true, decayPolicy);
+
+    arma::mat coordinates = f.GetInitialPoint();
+    double result = s.Optimize(f, coordinates);
+
+    // The final value of the objective function should be close to the optimal
+    // value, that is the sum of values at the vertices of the parabolas.
+    BOOST_REQUIRE_CLOSE(result, 123.75, 0.01);
+
+    // The co-ordinates should be the vertices of the parabolas.
+    BOOST_REQUIRE_CLOSE(coordinates[0], 2, 0.02);
+    BOOST_REQUIRE_CLOSE(coordinates[1], 1, 0.02);
+    BOOST_REQUIRE_CLOSE(coordinates[2], 1.5, 0.02);
+    BOOST_REQUIRE_CLOSE(coordinates[3], 4, 0.02);
+  }
+}
+
+/**
+ * When run with a single thread, parallel SGD should be identical to normal
+ * SGD.
+ */
+BOOST_AUTO_TEST_CASE(GeneralizedRosenbrockTest)
+{
+  // Loop over several variants.
+  for (size_t i = 10; i < 50; i += 5)
+  {
+    // Create the generalized Rosenbrock function.
+    GeneralizedRosenbrockFunction f(i);
+
+    ConstantStep decayPolicy(0.001);
+
+    ParallelSGD<ConstantStep> s(0, f.NumFunctions(), 1e-12, true, decayPolicy);
+
+    arma::mat coordinates = f.GetInitialPoint();
+
+    omp_set_num_threads(1);
+    double result = s.Optimize(f, coordinates);
+
+    BOOST_REQUIRE_SMALL(result, 1e-8);
+    for (size_t j = 0; j < i; ++j)
+      BOOST_REQUIRE_CLOSE(coordinates[j], (double) 1.0, 0.01);
+  }
+}
+
+/**
+ * Test the correctness of the Exponential backoff stepsize decay policy.
+ */
+BOOST_AUTO_TEST_CASE(ExponentialBackoffDecayTest)
+{
+  ExponentialBackoff decayPolicy(100, 100, 0.9);
+
+  // At the first iteration, stepsize should be unchanged
+  BOOST_REQUIRE_EQUAL(decayPolicy.StepSize(1), 100);
+  // At the 99th iteration, stepsize should be unchanged
+  BOOST_REQUIRE_EQUAL(decayPolicy.StepSize(99), 100);
+  // At the 100th iteration, stepsize should be changed
+  BOOST_REQUIRE_EQUAL(decayPolicy.StepSize(100), 90);
+  // At the 210th iteration, stepsize should be unchanged
+  BOOST_REQUIRE_EQUAL(decayPolicy.StepSize(210), 90);
+  // At the 211th iteration, stepsize should be changed
+  BOOST_REQUIRE_EQUAL(decayPolicy.StepSize(211), 81);
+}
+
+BOOST_AUTO_TEST_SUITE_END();

src/mlpack/tests/regularized_svd_test.cpp

@@ -11,12 +11,15 @@
  */
 #include <mlpack/core.hpp>
 #include <mlpack/methods/regularized_svd/regularized_svd.hpp>
+#include <mlpack/core/optimizers/parallel_sgd/parallel_sgd.hpp>
+#include <mlpack/core/optimizers/parallel_sgd/decay_policies/constant_step.hpp>
 
 #include <boost/test/unit_test.hpp>
 #include "test_tools.hpp"
 
 using namespace mlpack;
 using namespace mlpack::svd;
+using namespace mlpack::optimization;
 
 BOOST_AUTO_TEST_SUITE(RegularizedSVDTest);
 
@@ -41,7 +44,7 @@ BOOST_AUTO_TEST_CASE(RegularizedSVDFunctionRandomEvaluate)
   data(1, numRatings - 1) = numItems - 1;
 
   // Make a RegularizedSVDFunction with zero regularization.
-  RegularizedSVDFunction rSVDFunc(data, rank, 0);
+  RegularizedSVDFunction<arma::mat> rSVDFunc(data, rank, 0);
 
   for (size_t i = 0; i < numTrials; i++)
   {
@@ -89,9 +92,9 @@ BOOST_AUTO_TEST_CASE(RegularizedSVDFunctionRegularizationEvaluate)
 
   // Make three RegularizedSVDFunction objects with different amounts of
   // regularization.
-  RegularizedSVDFunction rSVDFuncNoReg(data, rank, 0);
-  RegularizedSVDFunction rSVDFuncSmallReg(data, rank, 0.5);
-  RegularizedSVDFunction rSVDFuncBigReg(data, rank, 20);
+  RegularizedSVDFunction<arma::mat> rSVDFuncNoReg(data, rank, 0);
+  RegularizedSVDFunction<arma::mat> rSVDFuncSmallReg(data, rank, 0.5);
+  RegularizedSVDFunction<arma::mat> rSVDFuncBigReg(data, rank, 20);
 
   for (size_t i = 0; i < numTrials; i++)
   {
@@ -146,8 +149,8 @@ BOOST_AUTO_TEST_CASE(RegularizedSVDFunctionGradient)
 
   // Make two RegularizedSVDFunction objects, one with regularization and one
   // without.
-  RegularizedSVDFunction rSVDFunc1(data, rank, 0);
-  RegularizedSVDFunction rSVDFunc2(data, rank, 0.5);
+  RegularizedSVDFunction<arma::mat> rSVDFunc1(data, rank, 0);
+  RegularizedSVDFunction<arma::mat> rSVDFunc2(data, rank, 0.5);
 
   // Calculate gradients for both the objects.
   arma::mat gradient1, gradient2;
@@ -225,7 +228,7 @@ BOOST_AUTO_TEST_CASE(RegularizedSVDFunctionOptimize)
   }
 
   // Make the Reg SVD function and the optimizer.
-  RegularizedSVDFunction rSVDFunc(data, rank, lambda);
+  RegularizedSVDFunction<arma::mat> rSVDFunc(data, rank, lambda);
   mlpack::optimization::StandardSGD optimizer(alpha, iterations * numRatings);
 
   // Obtain optimized parameters after training.
@@ -248,4 +251,65 @@ BOOST_AUTO_TEST_CASE(RegularizedSVDFunctionOptimize)
   BOOST_REQUIRE_SMALL(relativeError, 1e-2);
 }
 
+// Test Regularized SVD with parallel SGD.
+BOOST_AUTO_TEST_CASE(RegularizedSVDFunctionOptimizeHOGWILD)
+{
+  // Define useful constants.
+  const size_t numUsers = 50;
+  const size_t numItems = 50;
+  const size_t numRatings = 100;
+  const size_t rank = 10;
+  const double alpha = 0.01;
+  const double lambda = 0.01;
+
+  // Initiate random parameters.
+  arma::mat parameters = arma::randu(rank, numUsers + numItems);
+
+  // Make a random rating dataset.
+  arma::mat data = arma::randu(3, numRatings);
+  data.row(0) = floor(data.row(0) * numUsers);
+  data.row(1) = floor(data.row(1) * numItems);
+
+  // Manually set last row to maximum user and maximum item.
+  data(0, numRatings - 1) = numUsers - 1;
+  data(1, numRatings - 1) = numItems - 1;
+
+  // Make rating entries based on the parameters.
+  for (size_t i = 0; i < numRatings; i++)
+  {
+    data(2, i) = arma::dot(parameters.col(data(0, i)),
+                           parameters.col(numUsers + data(1, i)));
+  }
+
+  // Make the Reg SVD function and the optimizer.
+  RegularizedSVDFunction<arma::mat> rSVDFunc(data, rank, lambda);
+
+  ConstantStep decayPolicy(alpha);
+
+  // Iterate till convergence.
+  // The threadShareSize is chosen such that each function gets optimized.
+  ParallelSGD<ConstantStep> optimizer(0,
+      std::ceil((float) rSVDFunc.NumFunctions() / omp_get_max_threads()), 1e-5,
+      true, decayPolicy);
+
+  // Obtain optimized parameters after training.
+  arma::mat optParameters = arma::randu(rank, numUsers + numItems);
+  optimizer.Optimize(rSVDFunc, optParameters);
+
+  // Get predicted ratings from optimized parameters.
+  arma::mat predictedData(1, numRatings);
+  for (size_t i = 0; i < numRatings; i++)
+  {
+    predictedData(0, i) = arma::dot(optParameters.col(data(0, i)),
+                                    optParameters.col(numUsers + data(1, i)));
+  }
+
+  // Calculate relative error.
+  const double relativeError = arma::norm(data.row(2) - predictedData, "frob") /
+                               arma::norm(data, "frob");
+
+  // Relative error should be small.
+  BOOST_REQUIRE_SMALL(relativeError, 1e-2);
+}
+
 BOOST_AUTO_TEST_SUITE_END();