diff --git a/ChangeLog b/ChangeLog
index bb57404..b6a09b7 100644
--- a/ChangeLog
+++ b/ChangeLog
@@ -1,3 +1,15 @@
+2019-08-09  James Balamuta  <balamut2@illinois.edu>
+
+	* DESCRIPTION (Version, Date): Release 1.16.0
+
+	* NEWS.md: Update for Ensmallen release 1.16.0
+
+	* inst/include/ensmallen_bits: Upgraded to Ensmallen 1.16.0
+	* inst/include/ensmallen.hpp: ditto
+
+	* inst/include/ensmallen_bits/pso/pso.hpp: Fixed return type issue detected
+	  by gcc 9.
+
 2019-05-20  James Balamuta  <balamut2@illinois.edu>
 
 	* DESCRIPTION (Version, Date): Release 1.15.0
diff --git a/DESCRIPTION b/DESCRIPTION
index e19c34b..316cad1 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,6 +1,6 @@
 Package: RcppEnsmallen
 Title: Header-Only C++ Mathematical Optimization Library for 'Armadillo'
-Version: 0.1.15.0.1
+Version: 0.1.16.0.1
 Authors@R: c(
     person("James Joseph", "Balamuta", email = "balamut2@illinois.edu", 
            role = c("aut", "cre", "cph"), 
diff --git a/NEWS.md b/NEWS.md
index a4a7195..761e809 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -1,3 +1,14 @@
+# RcppEnsmallen 0.1.16.0.1
+
+- Upgraded to ensmallen release 1.16.0 "Loud Alarm Clock" (2019-08-09)
+  - Add option to avoid computing exact objective at the end of the optimization
+    ([#109](https://github.com/mlpack/ensmallen/pull/109)).
+  - Fix handling of curvature for BigBatchSGD ([#118](https://github.com/mlpack/ensmallen/pull/118)).
+  - Reduce runtime of tests ([#118](https://github.com/mlpack/ensmallen/pull/118)).
+  - Introduce local-best particle swarm optimization, `LBestPSO`, for
+    unconstrained optimization problems ([#86](https://github.com/mlpack/ensmallen/pull/86)).
+- Fix return type error in `PSO` ([#123](https://github.com/mlpack/ensmallen/pull/123))
+
 # RcppEnsmallen 0.1.15.0.1
 
 - Upgraded to ensmallen release 1.15.0 "Wrong Side Of The Road" (2019-05-14)
diff --git a/cran-comments.md b/cran-comments.md
index de89236..db95e71 100644
--- a/cran-comments.md
+++ b/cran-comments.md
@@ -1,7 +1,7 @@
 ## Test environments
 
-* local OS X install, R 3.6.0
-* ubuntu 14.04 (on travis-ci), R 3.6.0
+* local OS X install, R 3.6.1
+* ubuntu 14.04 (on travis-ci), R 3.6.1
 * win-builder (devel and release)
 
 ## R CMD check results
diff --git a/inst/include/ensmallen.hpp b/inst/include/ensmallen.hpp
index 18341a3..c685b57 100644
--- a/inst/include/ensmallen.hpp
+++ b/inst/include/ensmallen.hpp
@@ -83,6 +83,7 @@
 #include "ensmallen_bits/lbfgs/lbfgs.hpp"
 #include "ensmallen_bits/padam/padam.hpp"
 #include "ensmallen_bits/parallel_sgd/parallel_sgd.hpp"
+#include "ensmallen_bits/pso/pso.hpp"
 #include "ensmallen_bits/rmsprop/rmsprop.hpp"
 
 #include "ensmallen_bits/sa/sa.hpp"
diff --git a/inst/include/ensmallen_bits/ada_delta/ada_delta.hpp b/inst/include/ensmallen_bits/ada_delta/ada_delta.hpp
index f7dcdf8..75ee239 100644
--- a/inst/include/ensmallen_bits/ada_delta/ada_delta.hpp
+++ b/inst/include/ensmallen_bits/ada_delta/ada_delta.hpp
@@ -66,6 +66,8 @@ class AdaDelta
    *        function is visited in linear order.
    * @param resetPolicy If true, parameters are reset before every Optimize
    *        call; otherwise, their values are retained.
+   * @param exactObjective Calculate the exact objective (Default: estimate the
+   *        final objective obtained on the last pass over the data).
    */
   AdaDelta(const double stepSize = 1.0,
            const size_t batchSize = 32,
@@ -74,7 +76,8 @@ class AdaDelta
            const size_t maxIterations = 100000,
            const double tolerance = 1e-5,
            const bool shuffle = true,
-           const bool resetPolicy = true);
+           const bool resetPolicy = true,
+           const bool exactObjective = false);
 
   /**
    * Optimize the given function using AdaDelta. The given starting point will
@@ -128,6 +131,11 @@ class AdaDelta
   //! Modify whether or not the individual functions are shuffled.
   bool& Shuffle() { return optimizer.Shuffle(); }
 
+  //! Get whether or not the actual objective is calculated.
+  bool ExactObjective() const { return optimizer.ExactObjective(); }
+  //! Modify whether or not the actual objective is calculated.
+  bool& ExactObjective() { return optimizer.ExactObjective(); }
+
   //! Get whether or not the update policy parameters
   //! are reset before Optimize call.
   bool ResetPolicy() const { return optimizer.ResetPolicy(); }
diff --git a/inst/include/ensmallen_bits/ada_delta/ada_delta_impl.hpp b/inst/include/ensmallen_bits/ada_delta/ada_delta_impl.hpp
index 829761c..9e641f8 100644
--- a/inst/include/ensmallen_bits/ada_delta/ada_delta_impl.hpp
+++ b/inst/include/ensmallen_bits/ada_delta/ada_delta_impl.hpp
@@ -26,7 +26,8 @@ inline AdaDelta::AdaDelta(const double stepSize,
                           const size_t maxIterations,
                           const double tolerance,
                           const bool shuffle,
-                          const bool resetPolicy) :
+                          const bool resetPolicy,
+                          const bool exactObjective) :
     optimizer(stepSize,
               batchSize,
               maxIterations,
@@ -34,7 +35,8 @@ inline AdaDelta::AdaDelta(const double stepSize,
               shuffle,
               AdaDeltaUpdate(rho, epsilon),
               NoDecay(),
-              resetPolicy)
+              resetPolicy,
+              exactObjective)
 { /* Nothing to do. */ }
 
 } // namespace ens
diff --git a/inst/include/ensmallen_bits/ada_grad/ada_grad.hpp b/inst/include/ensmallen_bits/ada_grad/ada_grad.hpp
index 9f638a0..3890874 100644
--- a/inst/include/ensmallen_bits/ada_grad/ada_grad.hpp
+++ b/inst/include/ensmallen_bits/ada_grad/ada_grad.hpp
@@ -64,6 +64,8 @@ class AdaGrad
    *        function is visited in linear order.
    * @param resetPolicy If true, parameters are reset before every Optimize
    *        call; otherwise, their values are retained.
+   * @param exactObjective Calculate the exact objective (Default: estimate the
+   *        final objective obtained on the last pass over the data).
    */
   AdaGrad(const double stepSize = 0.01,
           const size_t batchSize = 32,
@@ -71,7 +73,8 @@ class AdaGrad
           const size_t maxIterations = 100000,
           const double tolerance = 1e-5,
           const bool shuffle = true,
-          const bool resetPolicy = true);
+          const bool resetPolicy = true,
+          const bool exactObjective = false);
 
   /**
    * Optimize the given function using AdaGrad. The given starting point will
@@ -119,6 +122,11 @@ class AdaGrad
   //! Modify whether or not the individual functions are shuffled.
   bool& Shuffle() { return optimizer.Shuffle(); }
 
+  //! Get whether or not the actual objective is calculated.
+  bool ExactObjective() const { return optimizer.ExactObjective(); }
+  //! Modify whether or not the actual objective is calculated.
+  bool& ExactObjective() { return optimizer.ExactObjective(); }
+
   //! Get whether or not the update policy parameters
   //! are reset before Optimize call.
   bool ResetPolicy() const { return optimizer.ResetPolicy(); }
diff --git a/inst/include/ensmallen_bits/ada_grad/ada_grad_impl.hpp b/inst/include/ensmallen_bits/ada_grad/ada_grad_impl.hpp
index d4f9870..f3b723f 100644
--- a/inst/include/ensmallen_bits/ada_grad/ada_grad_impl.hpp
+++ b/inst/include/ensmallen_bits/ada_grad/ada_grad_impl.hpp
@@ -23,7 +23,8 @@ inline AdaGrad::AdaGrad(const double stepSize,
                         const size_t maxIterations,
                         const double tolerance,
                         const bool shuffle,
-                        const bool resetPolicy) :
+                        const bool resetPolicy,
+                        const bool exactObjective) :
     optimizer(stepSize,
               batchSize,
               maxIterations,
@@ -31,7 +32,8 @@ inline AdaGrad::AdaGrad(const double stepSize,
               shuffle,
               AdaGradUpdate(epsilon),
               NoDecay(),
-              resetPolicy)
+              resetPolicy,
+              exactObjective)
 { /* Nothing to do. */ }
 
 } // namespace ens
diff --git a/inst/include/ensmallen_bits/adam/adam.hpp b/inst/include/ensmallen_bits/adam/adam.hpp
index 9a485d7..6085eb7 100644
--- a/inst/include/ensmallen_bits/adam/adam.hpp
+++ b/inst/include/ensmallen_bits/adam/adam.hpp
@@ -88,6 +88,8 @@ class AdamType
    *        function is visited in linear order.
    * @param resetPolicy If true, parameters are reset before every Optimize
    *        call; otherwise, their values are retained.
+   * @param exactObjective Calculate the exact objective (Default: estimate the
+   *        final objective obtained on the last pass over the data).
    */
   AdamType(const double stepSize = 0.001,
            const size_t batchSize = 32,
@@ -97,7 +99,8 @@ class AdamType
            const size_t maxIterations = 100000,
            const double tolerance = 1e-5,
            const bool shuffle = true,
-           const bool resetPolicy = true);
+           const bool resetPolicy = true,
+           const bool exactObjective = false);
 
   /**
    * Optimize the given function using Adam. The given starting point will be
@@ -155,6 +158,11 @@ class AdamType
   //! Modify whether or not the individual functions are shuffled.
   bool& Shuffle() { return optimizer.Shuffle(); }
 
+  //! Get whether or not the actual objective is calculated.
+  bool ExactObjective() const { return optimizer.ExactObjective(); }
+  //! Modify whether or not the actual objective is calculated.
+  bool& ExactObjective() { return optimizer.ExactObjective(); }
+
   //! Get whether or not the update policy parameters
   //! are reset before Optimize call.
   bool ResetPolicy() const { return optimizer.ResetPolicy(); }
diff --git a/inst/include/ensmallen_bits/adam/adam_impl.hpp b/inst/include/ensmallen_bits/adam/adam_impl.hpp
index a0daf83..7758f72 100644
--- a/inst/include/ensmallen_bits/adam/adam_impl.hpp
+++ b/inst/include/ensmallen_bits/adam/adam_impl.hpp
@@ -30,7 +30,8 @@ AdamType<UpdateRule>::AdamType(
     const size_t maxIterations,
     const double tolerance,
     const bool shuffle,
-    const bool resetPolicy) :
+    const bool resetPolicy,
+    const bool exactObjective) :
     optimizer(stepSize,
               batchSize,
               maxIterations,
@@ -38,7 +39,8 @@ AdamType<UpdateRule>::AdamType(
               shuffle,
               UpdateRule(epsilon, beta1, beta2),
               NoDecay(),
-              resetPolicy)
+              resetPolicy,
+              exactObjective)
 { /* Nothing to do. */ }
 
 } // namespace ens
diff --git a/inst/include/ensmallen_bits/bigbatch_sgd/adaptive_stepsize.hpp b/inst/include/ensmallen_bits/bigbatch_sgd/adaptive_stepsize.hpp
index c189837..0844c4d 100644
--- a/inst/include/ensmallen_bits/bigbatch_sgd/adaptive_stepsize.hpp
+++ b/inst/include/ensmallen_bits/bigbatch_sgd/adaptive_stepsize.hpp
@@ -126,20 +126,25 @@ class AdaptiveStepsize
 
     // Update sample variance & norm of the gradient.
     sampleVariance = vB;
-    gradientNorm = std::pow(arma::norm(gradient / backtrackingBatchSize, 2), 2.0);
-
-    // Compute curvature.
-    double v = arma::trace(arma::trans(iterate - iteratePrev) *
-        (gradient - gradPrevIterate)) /
-        std::pow(arma::norm(iterate - iteratePrev, 2), 2.0);
+    gradientNorm = std::pow(arma::norm(gradient / backtrackingBatchSize, 2),
+        2.0);
+
+    // Compute curvature.  If it can't be computed (typically due to floating
+    // point representation issues), call it 0.  If the curvature is 0, the step
+    // size will not decay.
+    const double vNum = arma::trace(arma::trans(iterate - iteratePrev) *
+        (gradient - gradPrevIterate));
+    const double vDenom = std::pow(arma::norm(iterate - iteratePrev, 2), 2.0);
+    const double vTmp = vNum / vDenom;
+    const double v = std::isfinite(vTmp) ? vTmp : 0.0;
 
     // Update previous iterate.
     iteratePrev = iterate;
 
     // TODO: Develop an absolute strategy to deal with stepSizeDecay updates in
-    // case we arrive at local minima. See #1469 for more details.
+    // case we arrive at local minima. See mlpack#1469 for more details.
     double stepSizeDecay = 0;
-    if (gradientNorm && sampleVariance && batchSize)
+    if (gradientNorm && sampleVariance && batchSize && v)
     {
       if (batchSize < function.NumFunctions())
       {
diff --git a/inst/include/ensmallen_bits/bigbatch_sgd/bigbatch_sgd.hpp b/inst/include/ensmallen_bits/bigbatch_sgd/bigbatch_sgd.hpp
index bdd51db..5b12230 100644
--- a/inst/include/ensmallen_bits/bigbatch_sgd/bigbatch_sgd.hpp
+++ b/inst/include/ensmallen_bits/bigbatch_sgd/bigbatch_sgd.hpp
@@ -91,13 +91,16 @@ class BigBatchSGD
    * @param tolerance Maximum absolute tolerance to terminate algorithm.
    * @param shuffle If true, the batch order is shuffled; otherwise, each
    *        batch is visited in linear order.
+   * @param exactObjective Calculate the exact objective (Default: estimate the
+   *        final objective obtained on the last pass over the data).
    */
   BigBatchSGD(const size_t batchSize = 1000,
               const double stepSize = 0.01,
               const double batchDelta = 0.1,
               const size_t maxIterations = 100000,
               const double tolerance = 1e-5,
-              const bool shuffle = true);
+              const bool shuffle = true,
+              const bool exactObjective = false);
   /**
    * Optimize the given function using big-batch SGD.  The given starting point
    * will be modified to store the finishing point of the algorithm, and the
@@ -147,6 +150,11 @@ class BigBatchSGD
   //! Modify the update policy.
   UpdatePolicyType& UpdatePolicy() { return updatePolicy; }
 
+  //! Get whether or not the actual objective is calculated.
+  bool ExactObjective() const { return exactObjective; }
+  //! Modify whether or not the actual objective is calculated.
+  bool& ExactObjective() { return exactObjective; }
+
  private:
   //! The size of the current batch.
   size_t batchSize;
@@ -167,6 +175,9 @@ class BigBatchSGD
   //! iterating.
   bool shuffle;
 
+  //! Controls whether or not the actual Objective value is calculated.
+  bool exactObjective;
+
   //! The update policy used to update the parameters in each iteration.
   UpdatePolicyType updatePolicy;
 };
diff --git a/inst/include/ensmallen_bits/bigbatch_sgd/bigbatch_sgd_impl.hpp b/inst/include/ensmallen_bits/bigbatch_sgd/bigbatch_sgd_impl.hpp
index 99cb49a..5626232 100644
--- a/inst/include/ensmallen_bits/bigbatch_sgd/bigbatch_sgd_impl.hpp
+++ b/inst/include/ensmallen_bits/bigbatch_sgd/bigbatch_sgd_impl.hpp
@@ -26,13 +26,15 @@ BigBatchSGD<UpdatePolicyType>::BigBatchSGD(
     const double batchDelta,
     const size_t maxIterations,
     const double tolerance,
-    const bool shuffle) :
+    const bool shuffle,
+    const bool exactObjective) :
     batchSize(batchSize),
     stepSize(stepSize),
     batchDelta(batchDelta),
     maxIterations(maxIterations),
     tolerance(tolerance),
     shuffle(shuffle),
+    exactObjective(exactObjective),
     updatePolicy(UpdatePolicyType())
 { /* Nothing to do. */ }
 
@@ -186,12 +188,15 @@ double BigBatchSGD<UpdatePolicyType>::Optimize(
   Info << "Big-batch SGD: maximum iterations (" << maxIterations << ") "
       << "reached; terminating optimization." << std::endl;
 
-  // Calculate final objective.
-  overallObjective = 0;
-  for (size_t i = 0; i < numFunctions; i += batchSize)
+  // Calculate final objective if exactObjective is set to true.
+  if (exactObjective)
   {
-    const size_t effectiveBatchSize = std::min(batchSize, numFunctions - i);
-    overallObjective += f.Evaluate(iterate, i, effectiveBatchSize);
+    overallObjective = 0;
+    for (size_t i = 0; i < numFunctions; i += batchSize)
+    {
+      const size_t effectiveBatchSize = std::min(batchSize, numFunctions - i);
+      overallObjective += f.Evaluate(iterate, i, effectiveBatchSize);
+    }
   }
   return overallObjective;
 }
diff --git a/inst/include/ensmallen_bits/ens_version.hpp b/inst/include/ensmallen_bits/ens_version.hpp
index 915de88..d816d28 100644
--- a/inst/include/ensmallen_bits/ens_version.hpp
+++ b/inst/include/ensmallen_bits/ens_version.hpp
@@ -15,13 +15,13 @@
 #define ENS_VERSION_MAJOR 1
 // The minor version is two digits so regular numerical comparisons of versions
 // work right.  The first minor version of a release is always 10.
-#define ENS_VERSION_MINOR 15
+#define ENS_VERSION_MINOR 16
 #define ENS_VERSION_PATCH 0
 // If this is a release candidate, it will be reflected in the version name
 // (i.e. the version name will be "RC1", "RC2", etc.).  Otherwise the version
 // name will typically be a seemingly arbitrary set of words that does not
 // contain the capitalized string "RC".
-#define ENS_VERSION_NAME "Wrong Side Of The Road"
+#define ENS_VERSION_NAME "Loud Alarm Clock"
 
 namespace ens {
 
diff --git a/inst/include/ensmallen_bits/eve/eve.hpp b/inst/include/ensmallen_bits/eve/eve.hpp
index 3d8d60b..437c964 100644
--- a/inst/include/ensmallen_bits/eve/eve.hpp
+++ b/inst/include/ensmallen_bits/eve/eve.hpp
@@ -72,6 +72,8 @@ class Eve
    * @param tolerance Maximum absolute tolerance to terminate algorithm.
    * @param shuffle If true, the function order is shuffled; otherwise, each
    *        function is visited in linear order.
+   * @param exactObjective Calculate the exact objective (Default: estimate the
+   *        final objective obtained on the last pass over the data).
    */
   Eve(const double stepSize = 0.001,
       const size_t batchSize = 32,
@@ -82,7 +84,8 @@ class Eve
       const double clip = 10,
       const size_t maxIterations = 100000,
       const double tolerance = 1e-5,
-      const bool shuffle = true);
+      const bool shuffle = true,
+      const bool exactObjective = false);
 
   /**
    * Optimize the given function using stochastic gradient descent.  The given
@@ -147,6 +150,11 @@ class Eve
   //! Modify whether or not the individual functions are shuffled.
   bool& Shuffle() { return shuffle; }
 
+  //! Get whether or not the actual objective is calculated.
+  bool ExactObjective() const { return exactObjective; }
+  //! Modify whether or not the actual objective is calculated.
+  bool& ExactObjective() { return exactObjective; }
+
  private:
   //! The step size for each example.
   double stepSize;
@@ -178,6 +186,9 @@ class Eve
   //! Controls whether or not the individual functions are shuffled when
   //! iterating.
   bool shuffle;
+
+  //! Controls whether or not the actual Objective value is calculated.
+  bool exactObjective;
 };
 
 } // namespace ens
diff --git a/inst/include/ensmallen_bits/eve/eve_impl.hpp b/inst/include/ensmallen_bits/eve/eve_impl.hpp
index f355fa2..93126e5 100644
--- a/inst/include/ensmallen_bits/eve/eve_impl.hpp
+++ b/inst/include/ensmallen_bits/eve/eve_impl.hpp
@@ -29,7 +29,8 @@ inline Eve::Eve(const double stepSize,
                 const double clip,
                 const size_t maxIterations,
                 const double tolerance,
-                const bool shuffle) :
+                const bool shuffle,
+                const bool exactObjective) :
     stepSize(stepSize),
     batchSize(batchSize),
     beta1(beta1),
@@ -39,7 +40,8 @@ inline Eve::Eve(const double stepSize,
     clip(clip),
     maxIterations(maxIterations),
     tolerance(tolerance),
-    shuffle(shuffle)
+    shuffle(shuffle),
+    exactObjective(exactObjective)
 { /* Nothing to do. */ }
 
 //! Optimize the function (minimize).
@@ -154,12 +156,15 @@ double Eve::Optimize(
   Info << "Eve: maximum iterations (" << maxIterations << ") reached; "
       << "terminating optimization." << std::endl;
 
-  // Calculate final objective.
-  overallObjective = 0;
-  for (size_t i = 0; i < numFunctions; i += batchSize)
+  // Calculate final objective if exactObjective is set to true.
+  if (exactObjective)
   {
-    const size_t effectiveBatchSize = std::min(batchSize, numFunctions - i);
-    overallObjective += f.Evaluate(iterate, i, effectiveBatchSize);
+    overallObjective = 0;
+    for (size_t i = 0; i < numFunctions; i += batchSize)
+    {
+      const size_t effectiveBatchSize = std::min(batchSize, numFunctions - i);
+      overallObjective += f.Evaluate(iterate, i, effectiveBatchSize);
+    }
   }
   return overallObjective;
 }
diff --git a/inst/include/ensmallen_bits/ftml/ftml.hpp b/inst/include/ensmallen_bits/ftml/ftml.hpp
index 87dd3b1..85442c2 100644
--- a/inst/include/ensmallen_bits/ftml/ftml.hpp
+++ b/inst/include/ensmallen_bits/ftml/ftml.hpp
@@ -66,6 +66,8 @@ class FTML
    *        function is visited in linear order.
    * @param resetPolicy If true, parameters are reset before every Optimize
    *        call; otherwise, their values are retained.
+   * @param exactObjective Calculate the exact objective (Default: estimate the
+   *        final objective obtained on the last pass over the data).
    */
   FTML(const double stepSize = 0.001,
        const size_t batchSize = 32,
@@ -75,7 +77,8 @@ class FTML
        const size_t maxIterations = 100000,
        const double tolerance = 1e-5,
        const bool shuffle = true,
-       const bool resetPolicy = true);
+       const bool resetPolicy = true,
+       const bool exactObjective = false);
 
   /**
    * Optimize the given function using FTML. The given starting point will
@@ -133,6 +136,11 @@ class FTML
   //! Modify whether or not the individual functions are shuffled.
   bool& Shuffle() { return optimizer.Shuffle(); }
 
+  //! Get whether or not the actual objective is calculated.
+  bool ExactObjective() const { return optimizer.ExactObjective(); }
+  //! Modify whether or not the actual objective is calculated.
+  bool& ExactObjective() { return optimizer.ExactObjective(); }
+
   //! Get whether or not the update policy parameters
   //! are reset before Optimize call.
   bool ResetPolicy() const { return optimizer.ResetPolicy(); }
diff --git a/inst/include/ensmallen_bits/ftml/ftml_impl.hpp b/inst/include/ensmallen_bits/ftml/ftml_impl.hpp
index 47701fb..05f0576 100644
--- a/inst/include/ensmallen_bits/ftml/ftml_impl.hpp
+++ b/inst/include/ensmallen_bits/ftml/ftml_impl.hpp
@@ -25,7 +25,8 @@ inline FTML::FTML(const double stepSize,
                   const size_t maxIterations,
                   const double tolerance,
                   const bool shuffle,
-                  const bool resetPolicy) :
+                  const bool resetPolicy,
+                  const bool exactObjective) :
     optimizer(stepSize,
               batchSize,
               maxIterations,
@@ -33,7 +34,8 @@ inline FTML::FTML(const double stepSize,
               shuffle,
               FTMLUpdate(epsilon, beta1, beta2),
               NoDecay(),
-              resetPolicy)
+              resetPolicy,
+              exactObjective)
 { /* Nothing to do. */ }
 
 } // namespace ens
diff --git a/inst/include/ensmallen_bits/fw/proximal/proximal_impl.hpp b/inst/include/ensmallen_bits/fw/proximal/proximal_impl.hpp
index f08e8e2..ca109c1 100644
--- a/inst/include/ensmallen_bits/fw/proximal/proximal_impl.hpp
+++ b/inst/include/ensmallen_bits/fw/proximal/proximal_impl.hpp
@@ -44,7 +44,7 @@ inline void Proximal::ProjectToL1Ball(arma::vec& v, double tau)
   arma::vec simplexSum = arma::cumsum(simplexSol);
 
   double nu = 0;
-  size_t rho;
+  size_t rho = 0;
   for (size_t j = 1; j <= simplexSol.n_rows; j++)
   {
     rho = simplexSol.n_rows - j;
diff --git a/inst/include/ensmallen_bits/katyusha/katyusha.hpp b/inst/include/ensmallen_bits/katyusha/katyusha.hpp
index d245db1..8cabfbf 100644
--- a/inst/include/ensmallen_bits/katyusha/katyusha.hpp
+++ b/inst/include/ensmallen_bits/katyusha/katyusha.hpp
@@ -63,6 +63,8 @@ class KatyushaType
    * @param tolerance Maximum absolute tolerance to terminate algorithm.
    * @param shuffle If true, the function order is shuffled; otherwise, each
    *    function is visited in linear order.
+   * @param exactObjective Calculate the exact objective (Default: estimate the
+   *        final objective obtained on the last pass over the data).
    */
   KatyushaType(const double convexity = 1.0,
                const double lipschitz = 10.0,
@@ -70,7 +72,8 @@ class KatyushaType
                const size_t maxIterations = 1000,
                const size_t innerIterations = 0,
                const double tolerance = 1e-5,
-               const bool shuffle = true);
+               const bool shuffle = true,
+               const bool exactObjective = false);
 
   /**
    * Optimize the given function using Katyusha. The given starting point will
@@ -120,6 +123,11 @@ class KatyushaType
   //! Modify whether or not the individual functions are shuffled.
   bool& Shuffle() { return shuffle; }
 
+  //! Get whether or not the actual objective is calculated.
+  bool ExactObjective() const { return exactObjective; }
+  //! Modify whether or not the actual objective is calculated.
+  bool& ExactObjective() { return exactObjective; }
+
  private:
   //! The convexity regularization term.
   double convexity;
@@ -142,6 +150,9 @@ class KatyushaType
   //! Controls whether or not the individual functions are shuffled when
   //! iterating.
   bool shuffle;
+
+  //! Controls whether or not the actual Objective value is calculated.
+  bool exactObjective;
 };
 
 // Convenience typedefs.
diff --git a/inst/include/ensmallen_bits/katyusha/katyusha_impl.hpp b/inst/include/ensmallen_bits/katyusha/katyusha_impl.hpp
index 1e9f9c7..b7bd840 100644
--- a/inst/include/ensmallen_bits/katyusha/katyusha_impl.hpp
+++ b/inst/include/ensmallen_bits/katyusha/katyusha_impl.hpp
@@ -27,14 +27,16 @@ KatyushaType<Proximal>::KatyushaType(
     const size_t maxIterations,
     const size_t innerIterations,
     const double tolerance,
-    const bool shuffle) :
+    const bool shuffle,
+    const bool exactObjective) :
     convexity(convexity),
     lipschitz(lipschitz),
     batchSize(batchSize),
     maxIterations(maxIterations),
     innerIterations(innerIterations),
     tolerance(tolerance),
-    shuffle(shuffle)
+    shuffle(shuffle),
+    exactObjective(exactObjective)
 { /* Nothing to do. */ }
 
 //! Optimize the function (minimize).
@@ -197,12 +199,15 @@ double KatyushaType<Proximal>::Optimize(
   Info << "Katyusha: maximum iterations (" << maxIterations << ") reached"
       << "; terminating optimization." << std::endl;
 
-  // Calculate final objective.
-  overallObjective = 0;
-  for (size_t i = 0; i < numFunctions; i += batchSize)
+  // Calculate final objective if exactObjective is set to true.
+  if (exactObjective)
   {
-    const size_t effectiveBatchSize = std::min(batchSize, numFunctions - i);
-    overallObjective += function.Evaluate(iterate, i, effectiveBatchSize);
+    overallObjective = 0;
+    for (size_t i = 0; i < numFunctions; i += batchSize)
+    {
+      const size_t effectiveBatchSize = std::min(batchSize, numFunctions - i);
+      overallObjective += function.Evaluate(iterate, i, effectiveBatchSize);
+    }
   }
   return overallObjective;
 }
diff --git a/inst/include/ensmallen_bits/padam/padam.hpp b/inst/include/ensmallen_bits/padam/padam.hpp
index c3c24ae..21f39cd 100644
--- a/inst/include/ensmallen_bits/padam/padam.hpp
+++ b/inst/include/ensmallen_bits/padam/padam.hpp
@@ -64,6 +64,8 @@ class Padam
    *        function is visited in linear order.
    * @param resetPolicy If true, parameters are reset before every Optimize
    *        call; otherwise, their values are retained.
+   * @param exactObjective Calculate the exact objective (Default: estimate the
+   *        final objective obtained on the last pass over the data).
    */
   Padam(const double stepSize = 0.001,
         const size_t batchSize = 32,
@@ -74,7 +76,8 @@ class Padam
         const size_t maxIterations = 100000,
         const double tolerance = 1e-5,
         const bool shuffle = true,
-        const bool resetPolicy = true) :
+        const bool resetPolicy = true,
+        const bool exactObjective = false) :
       optimizer(stepSize,
                 batchSize,
                 maxIterations,
@@ -82,7 +85,8 @@ class Padam
                 shuffle,
                 PadamUpdate(epsilon, beta1, beta2, partial),
                 NoDecay(),
-                resetPolicy)
+                resetPolicy,
+                exactObjective)
   { /* Nothing to do here. */ }
 
   /**
@@ -146,6 +150,11 @@ class Padam
   //! Modify whether or not the individual functions are shuffled.
   bool& Shuffle() { return optimizer.Shuffle(); }
 
+  //! Get whether or not the actual objective is calculated.
+  bool ExactObjective() const { return optimizer.ExactObjective(); }
+  //! Modify whether or not the actual objective is calculated.
+  bool& ExactObjective() { return optimizer.ExactObjective(); }
+
   //! Get whether or not the update policy parameters
   //! are reset before Optimize call.
   bool ResetPolicy() const { return optimizer.ResetPolicy(); }
diff --git a/inst/include/ensmallen_bits/problems/ackley_function.hpp b/inst/include/ensmallen_bits/problems/ackley_function.hpp
new file mode 100644
index 0000000..2c2b98b
--- /dev/null
+++ b/inst/include/ensmallen_bits/problems/ackley_function.hpp
@@ -0,0 +1,127 @@
+/**
+ * @file ackley_function.hpp
+ * @author Suryoday Basak
+ *
+ * Definition of the Ackley function.
+ *
+ * ensmallen 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 ensmallen.  If not, see
+ * http://www.opensource.org/licenses/BSD-3-Clause for more information.
+ */
+#ifndef ENSMALLEN_PROBLEMS_ACKLEY_FUNCTION_HPP
+#define ENSMALLEN_PROBLEMS_ACKLEY_FUNCTION_HPP
+
+namespace ens {
+namespace test {
+
+/**
+ * The Ackley function, defined by
+ *
+ * \f[
+ * f(x_1,x_2) = -20 * e^(-0.2 * sqrt(0.5 * (x_1^2 + x_2^2))) -
+ *    e * (0.5(cos(2 * pi * x_1) + cos(2 * pi * x_2))) + e + 20
+ * \f]
+ *
+ * This should optimize to f(x) = 0, at x = [0, 0].
+ *
+ * For more information, please refer to:
+ *
+ * @code
+ * @book{Ackley1987,
+ *   doi       = {10.1007/978-1-4613-1997-9},
+ *   url       = {https://doi.org/10.1007/978-1-4613-1997-9},
+ *   year      = {1987},
+ *   publisher = {Springer {US}},
+ *   author    = {David H. Ackley},
+ *   title     = {A Connectionist Machine for Genetic Hillclimbing}
+ * }
+ * @endcode
+ */
+class AckleyFunction
+{
+ public:
+  /**
+   * Initialize the AckleyFunction.
+   *
+   * @param c Multiplicative constant with a default value of 2 * pi.
+   * @param epsilon Coefficient to avoid division by zero (numerical stability).
+   */
+  AckleyFunction(const double c = 2 * arma::datum::pi,
+                 const double epsilon = 1e-8);
+
+  /**
+   * Shuffle the order of function visitation. This may be called by the
+   * optimizer.
+   */
+  void Shuffle();
+
+  //! Return 1 (the number of functions).
+  size_t NumFunctions() const { return 1; }
+
+  //! Get the starting point.
+  arma::mat GetInitialPoint() const { return arma::mat("-5.0; 5.0"); }
+
+  /**
+   * Evaluate a function for a particular batch-size.
+   *
+   * @param coordinates The function coordinates.
+   * @param begin The first function.
+   * @param batchSize Number of points to process.
+   */
+  double Evaluate(const arma::mat& coordinates,
+                  const size_t begin,
+                  const size_t batchSize) const;
+
+  /**
+   * Evaluate a function with the given coordinates.
+   *
+   * @param coordinates The function coordinates.
+   */
+  double Evaluate(const arma::mat& coordinates) const;
+
+  /**
+   * Evaluate the gradient of a function for a particular batch-size.
+   *
+   * @param coordinates The function coordinates.
+   * @param begin The first function.
+   * @param gradient The function gradient.
+   * @param batchSize Number of points to process.
+   */
+  void Gradient(const arma::mat& coordinates,
+                const size_t begin,
+                arma::mat& gradient,
+                const size_t batchSize) const;
+
+  /**
+   * Evaluate the gradient of a function with the given coordinates.
+   *
+   * @param coordinates The function coordinates.
+   * @param gradient The function gradient.
+   */
+  void Gradient(const arma::mat& coordinates, arma::mat& gradient);
+
+  //! Get the value used for c.
+  double MultiplicativeConstant() const { return c; }
+  //! Modify the value used for c.
+  double& MultiplicativeConstant() { return c; }
+
+  //! Get the value used for numerical stability.
+  double Epsilon() const { return epsilon; }
+  //! Modify the value used for numerical stability.
+  double& Epsilon() { return epsilon; }
+
+ private:
+  //! The value of the multiplicative constant.
+  double c;
+  //! The value used for numerical stability.
+  double epsilon;
+};
+
+} // namespace test
+} // namespace ens
+
+// Include implementation.
+#include "ackley_function_impl.hpp"
+
+#endif // ENSMALLEN_PROBLEMS_ACKLEY_FUNCTION_HPP
diff --git a/inst/include/ensmallen_bits/problems/ackley_function_impl.hpp b/inst/include/ensmallen_bits/problems/ackley_function_impl.hpp
new file mode 100644
index 0000000..327fb29
--- /dev/null
+++ b/inst/include/ensmallen_bits/problems/ackley_function_impl.hpp
@@ -0,0 +1,76 @@
+/**
+ * @file ackley_function_impl.hpp
+ * @author Suryoday Basak
+ *
+ * Implementation of the Ackley function.
+ *
+ * ensmallen 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 ensmallen.  If not, see
+ * http://www.opensource.org/licenses/BSD-3-Clause for more information.
+ */
+#ifndef ENSMALLEN_PROBLEMS_ACKLEY_FUNCTION_IMPL_HPP
+#define ENSMALLEN_PROBLEMS_ACKLEY_FUNCTION_IMPL_HPP
+
+// In case it hasn't been included yet.
+#include "ackley_function.hpp"
+using namespace std;
+
+namespace ens {
+namespace test {
+
+inline AckleyFunction::AckleyFunction(const double c, const double epsilon) :
+    c(c), epsilon(epsilon)
+{ /* Nothing to do here */}
+
+inline void AckleyFunction::Shuffle() { /* Nothing to do here */ }
+
+inline double AckleyFunction::Evaluate(const arma::mat& coordinates,
+                                       const size_t /* begin */,
+                                       const size_t /* batchSize */) const
+{
+  // For convenience; we assume these temporaries will be optimized out.
+  const double x1 = coordinates(0);
+  const double x2 = coordinates(1);
+
+  const double objective = -20 * exp(-0.2 * sqrt(0.5 * (x1 * x1 + x2 * x2))) -
+      exp(0.5 * (cos(c * x1) + cos(c * x2))) + exp(1) + 20;
+
+  return objective;
+}
+
+inline double AckleyFunction::Evaluate(const arma::mat& coordinates) const
+{
+  return Evaluate(coordinates, 0, NumFunctions());
+}
+
+inline void AckleyFunction::Gradient(const arma::mat& coordinates,
+                                     const size_t /* begin */,
+                                     arma::mat& gradient,
+                                     const size_t /* batchSize */) const
+{
+  // For convenience; we assume these temporaries will be optimized out.
+  const double x1 = coordinates(0);
+  const double x2 = coordinates(1);
+
+  // Aliases for different terms in the expression of the gradient.
+  const double t0 = sqrt(0.5 * (x1 * x1 + x2 * x2));
+  const double t1 = 2.0 * exp(- 0.2 * t0) / (t0 + epsilon);
+  const double t2 = 0.5 * c * exp(0.5 * (cos(c * x1) +
+      cos(c * x2)));
+
+  gradient.set_size(2, 1);
+  gradient(0) = (x1 * t1) + (t2 * sin(c * x1));
+  gradient(1) = (x2 * t1) + (t2 * sin(c * x2));
+}
+
+inline void AckleyFunction::Gradient(const arma::mat& coordinates,
+                                     arma::mat& gradient)
+{
+  Gradient(coordinates, 0, gradient, 1);
+}
+
+} // namespace test
+} // namespace ens
+
+#endif
diff --git a/inst/include/ensmallen_bits/problems/beale_function.hpp b/inst/include/ensmallen_bits/problems/beale_function.hpp
new file mode 100644
index 0000000..294ae84
--- /dev/null
+++ b/inst/include/ensmallen_bits/problems/beale_function.hpp
@@ -0,0 +1,104 @@
+/**
+ * @file beale_function.hpp
+ * @author Suryoday Basak
+ *
+ * Definition of the Beale function.
+ *
+ * ensmallen 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 ensmallen.  If not, see
+ * http://www.opensource.org/licenses/BSD-3-Clause for more information.
+ */
+#ifndef ENSMALLEN_PROBLEMS_BEALE_FUNCTION_HPP
+#define ENSMALLEN_PROBLEMS_BEALE_FUNCTION_HPP
+
+namespace ens {
+namespace test {
+
+/**
+ * The Beale function, defined by
+ *
+ * \f[
+ * f(x_1,x_2) = (1.5 - x_1 + x_1 * x_2)^2 +
+ *              (2.25 - x_1 + x_1 * x_2^2)^2 +
+ *              (2.625 - x_1 + x_1 * x_2^3)^2
+ * \f]
+ *
+ * This should optimize to f(x) = 0, at x = [3, 0.5].
+ *
+ * For more information, please refer to:
+ *
+ * @code
+ * @misc{1307.5838,
+ *   Author = {Masoumeh Vali},
+ *   Title  = {Rotational Mutation Genetic Algorithm on optimizationProblems},
+ *   Year   = {2013},
+ *   Eprint = {arXiv:1307.5838},
+ * }
+ * @endcode
+ */
+class BealeFunction
+{
+ public:
+  //! Initialize the BealeFunction.
+  BealeFunction();
+
+  /**
+   * Shuffle the order of function visitation. This may be called by the
+   * optimizer.
+   */
+  void Shuffle();
+
+  //! Return 1 (the number of functions).
+  size_t NumFunctions() const { return 1; }
+
+  //! Get the starting point.
+  arma::mat GetInitialPoint() const { return arma::mat("-4.5; 4.5"); }
+
+  /**
+   * Evaluate a function for a particular batch-size.
+   *
+   * @param coordinates The function coordinates.
+   * @param begin The first function.
+   * @param batchSize Number of points to process.
+   */
+  double Evaluate(const arma::mat& coordinates,
+                  const size_t begin,
+                  const size_t batchSize) const;
+
+  /**
+   * Evaluate a function with the given coordinates.
+   *
+   * @param coordinates The function coordinates.
+   */
+  double Evaluate(const arma::mat& coordinates) const;
+
+  /**
+   * Evaluate the gradient of a function for a particular batch-size.
+   *
+   * @param coordinates The function coordinates.
+   * @param begin The first function.
+   * @param gradient The function gradient.
+   * @param batchSize Number of points to process.
+   */
+  void Gradient(const arma::mat& coordinates,
+                const size_t begin,
+                arma::mat& gradient,
+                const size_t batchSize) const;
+
+  /**
+   * Evaluate the gradient of a function with the given coordinates.
+   *
+   * @param coordinates The function coordinates.
+   * @param gradient The function gradient.
+   */
+  void Gradient(const arma::mat& coordinates, arma::mat& gradient);
+};
+
+} // namespace test
+} // namespace ens
+
+// Include implementation.
+#include "beale_function_impl.hpp"
+
+#endif // ENSMALLEN_PROBLEMS_BEALE_FUNCTION_HPP
diff --git a/inst/include/ensmallen_bits/problems/beale_function_impl.hpp b/inst/include/ensmallen_bits/problems/beale_function_impl.hpp
new file mode 100644
index 0000000..451e6e7
--- /dev/null
+++ b/inst/include/ensmallen_bits/problems/beale_function_impl.hpp
@@ -0,0 +1,76 @@
+/**
+ * @file beale_function_impl.hpp
+ * @author Suryoday Basak
+ *
+ * Implementation of the Beale function.
+ *
+ * ensmallen 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 ensmallen.  If not, see
+ * http://www.opensource.org/licenses/BSD-3-Clause for more information.
+ */
+#ifndef ENSMALLEN_PROBLEMS_BEALE_FUNCTION_IMPL_HPP
+#define ENSMALLEN_PROBLEMS_BEALE_FUNCTION_IMPL_HPP
+
+// In case it hasn't been included yet.
+#include "beale_function.hpp"
+using namespace std;
+
+namespace ens {
+namespace test {
+
+inline BealeFunction::BealeFunction() { /* Nothing to do here */ }
+
+inline void BealeFunction::Shuffle() { /* Nothing to do here */ }
+
+inline double BealeFunction::Evaluate(const arma::mat& coordinates,
+                                      const size_t /* begin */,
+                                      const size_t /* batchSize */) const
+{
+  // For convenience; we assume these temporaries will be optimized out.
+  const double x1 = coordinates(0);
+  const double x2 = coordinates(1);
+
+  const double objective = pow(1.5 - x1 + x1 * x2, 2) +
+      pow(2.25 - x1 + x1 * x2 * x2, 2) +  pow(2.625 - x1 + x1 * pow(x2, 3), 2);
+
+  return objective;
+}
+
+inline double BealeFunction::Evaluate(const arma::mat& coordinates) const
+{
+  return Evaluate(coordinates, 0, NumFunctions());
+}
+
+inline void BealeFunction::Gradient(const arma::mat& coordinates,
+                                    const size_t /* begin */,
+                                    arma::mat& gradient,
+                                    const size_t /* batchSize */) const
+{
+  // For convenience; we assume these temporaries will be optimized out.
+  const double x1 = coordinates(0);
+  const double x2 = coordinates(1);
+
+  // Aliases for different terms in the expression of the gradient
+  const double x2Sq = x2 * x2;
+  const double x2Cub = pow(x2, 3);
+
+  gradient.set_size(2, 1);
+  gradient(0) = ((2 * x2 - 2) * (x1 * x2 - x1 + 1.5)) +
+      ((2 * x2Sq - 2) * (x1 * x2Sq - x1 + 2.25)) +
+      ((2 * x2Cub - 2) * (x1 * x2Cub - x1 + 2.625));
+  gradient(1) = (6 * x1 * x2Sq * (x1 * x2Cub - x1 + 2.625)) +
+      (4 * x1 * x2 * (x1 * x2Sq - x1 + 2.25)) +
+      (2 * x1 * (x1 * x2 - x1 + 1.5));
+}
+
+inline void BealeFunction::Gradient(const arma::mat& coordinates,
+                                    arma::mat& gradient)
+{
+  Gradient(coordinates, 0, gradient, 1);
+}
+
+} // namespace test
+} // namespace ens
+
+#endif
diff --git a/inst/include/ensmallen_bits/problems/booth_function.hpp b/inst/include/ensmallen_bits/problems/booth_function.hpp
index 43eeb90..5d360b4 100644
--- a/inst/include/ensmallen_bits/problems/booth_function.hpp
+++ b/inst/include/ensmallen_bits/problems/booth_function.hpp
@@ -44,9 +44,9 @@ class BoothFunction
   BoothFunction();
 
   /**
-  * Shuffle the order of function visitation. This may be called by the
-  * optimizer.
-  */
+   * Shuffle the order of function visitation. This may be called by the
+   * optimizer.
+   */
   void Shuffle();
 
   //! Return 1 (the number of functions).
@@ -55,7 +55,7 @@ class BoothFunction
   //! Get the starting point.
   arma::mat GetInitialPoint() const { return arma::mat("-9; -9"); }
 
-  /*
+  /**
    * Evaluate a function for a particular batch-size.
    *
    * @param coordinates The function coordinates.
@@ -66,14 +66,14 @@ class BoothFunction
                   const size_t begin,
                   const size_t batchSize) const;
 
-  /*
+  /**
    * Evaluate a function with the given coordinates.
    *
    * @param coordinates The function coordinates.
    */
   double Evaluate(const arma::mat& coordinates) const;
 
-  /*
+  /**
    * Evaluate the gradient of a function for a particular batch-size.
    *
    * @param coordinates The function coordinates.
@@ -86,7 +86,7 @@ class BoothFunction
                 arma::mat& gradient,
                 const size_t batchSize) const;
 
-  /*
+  /**
    * Evaluate the gradient of a function with the given coordinates.
    *
    * @param coordinates The function coordinates.
diff --git a/inst/include/ensmallen_bits/problems/bukin_function.hpp b/inst/include/ensmallen_bits/problems/bukin_function.hpp
index 86b5309..4bc949c 100644
--- a/inst/include/ensmallen_bits/problems/bukin_function.hpp
+++ b/inst/include/ensmallen_bits/problems/bukin_function.hpp
@@ -41,7 +41,7 @@ namespace test {
 class BukinFunction
 {
  public:
-  /*
+  /**
    * Initialize the BukinFunction.
    *
    * @param epsilon Coefficient to avoid division by zero (numerical stability).
@@ -49,9 +49,9 @@ class BukinFunction
   BukinFunction(const double epsilon = 1e-8);
 
   /**
-  * Shuffle the order of function visitation. This may be called by the
-  * optimizer.
-  */
+   * Shuffle the order of function visitation. This may be called by the
+   * optimizer.
+   */
   void Shuffle();
 
   //! Return 1 (the number of functions).
@@ -60,7 +60,7 @@ class BukinFunction
   //! Get the starting point.
   arma::mat GetInitialPoint() const { return arma::mat("-10; -2.0"); }
 
-  /*
+  /**
    * Evaluate a function for a particular batch-size.
    *
    * @param coordinates The function coordinates.
@@ -71,14 +71,14 @@ class BukinFunction
                   const size_t begin,
                   const size_t batchSize) const;
 
-  /*
+  /**
    * Evaluate a function with the given coordinates.
    *
    * @param coordinates The function coordinates.
    */
   double Evaluate(const arma::mat& coordinates) const;
 
-  /*
+  /**
    * Evaluate the gradient of a function for a particular batch-size.
    *
    * @param coordinates The function coordinates.
@@ -91,7 +91,7 @@ class BukinFunction
                 arma::mat& gradient,
                 const size_t batchSize) const;
 
-  /*
+  /**
    * Evaluate the gradient of a function with the given coordinates.
    *
    * @param coordinates The function coordinates.
diff --git a/inst/include/ensmallen_bits/problems/colville_function.hpp b/inst/include/ensmallen_bits/problems/colville_function.hpp
index c390ae3..910ed2d 100644
--- a/inst/include/ensmallen_bits/problems/colville_function.hpp
+++ b/inst/include/ensmallen_bits/problems/colville_function.hpp
@@ -45,9 +45,9 @@ class ColvilleFunction
   ColvilleFunction();
 
   /**
-  * Shuffle the order of function visitation. This may be called by the
-  * optimizer.
-  */
+   * Shuffle the order of function visitation. This may be called by the
+   * optimizer.
+   */
   void Shuffle();
 
   //! Return 1 (the number of functions).
@@ -56,7 +56,7 @@ class ColvilleFunction
   //! Get the starting point.
   arma::mat GetInitialPoint() const { return arma::mat("-5; 3; 1; -9"); }
 
-  /*
+  /**
    * Evaluate a function for a particular batch-size.
    *
    * @param coordinates The function coordinates.
@@ -67,14 +67,14 @@ class ColvilleFunction
                   const size_t begin,
                   const size_t batchSize) const;
 
-  /*
+  /**
    * Evaluate a function with the given coordinates.
    *
    * @param coordinates The function coordinates.
    */
   double Evaluate(const arma::mat& coordinates) const;
 
-  /*
+  /**
    * Evaluate the gradient of a function for a particular batch-size.
    *
    * @param coordinates The function coordinates.
@@ -87,7 +87,7 @@ class ColvilleFunction
                 arma::mat& gradient,
                 const size_t batchSize) const;
 
-  /*
+  /**
    * Evaluate the gradient of a function with the given coordinates.
    *
    * @param coordinates The function coordinates.
diff --git a/inst/include/ensmallen_bits/problems/cross_in_tray_function.hpp b/inst/include/ensmallen_bits/problems/cross_in_tray_function.hpp
new file mode 100644
index 0000000..2450c13
--- /dev/null
+++ b/inst/include/ensmallen_bits/problems/cross_in_tray_function.hpp
@@ -0,0 +1,88 @@
+/**
+ * @file cross_in_tray_function.hpp
+ * @author Suryoday Basak
+ *
+ * Definition of the Cross-in-Tray function.
+ *
+ * ensmallen 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 ensmallen.  If not, see
+ * http://www.opensource.org/licenses/BSD-3-Clause for more information.
+ */
+#ifndef ENSMALLEN_PROBLEMS_CROSS_IN_TRAY_FUNCTION_HPP
+#define ENSMALLEN_PROBLEMS_CROSS_IN_TRAY_FUNCTION_HPP
+
+namespace ens {
+namespace test {
+
+/**
+ * The Cross-in-Tray function, defined by
+ *
+ * \f[
+ * f(x1, x2) = - 0.0001 * (|sin(x1) * sin(x2) *
+ *               exp(|100 - (sqrt(x1^2 + x2^2) / pi)|)| + 1)^0.1
+ * \f]
+ *
+ * This should optimize to f(x1, x2) = -2.06261, at
+ *                                (x1, x2) = [1.34941, -1.34941], or
+ *                                (x1, x2) = [1.34941, 1.34941],  or
+ *                                (x1, x2) = [-1.34941, 1.34941], or
+ *                                (x1, x2) = [-1.34941, -1.34941]
+ *
+ * For more information, please refer to:
+ *
+ * @code
+ * @article{1308.4008,
+ *   Author = {Momin Jamil and Xin-She Yang},
+ *   Title  = {A Literature Survey of Benchmark Functions For Global
+ *             Optimization Problems},
+ *   Year   = {2013},
+ *   Eprint = {arXiv:1308.4008},
+ *   Doi    = {10.1504/IJMMNO.2013.055204},
+ * }
+ * @endcode
+ */
+class CrossInTrayFunction
+{
+ public:
+  //! Initialize the CrossInTrayFunction.
+  CrossInTrayFunction();
+
+  /**
+   * Shuffle the order of function visitation. This may be called by the
+   * optimizer.
+   */
+  void Shuffle();
+
+  //! Return 1 (the number of functions).
+  size_t NumFunctions() const { return 1; }
+
+  //! Get the starting point.
+  arma::mat GetInitialPoint() const { return arma::mat("0; 0"); }
+
+  /*
+   * Evaluate a function for a particular batch-size.
+   *
+   * @param coordinates The function coordinates.
+   * @param begin The first function.
+   * @param batchSize Number of points to process.
+   */
+  double Evaluate(const arma::mat& coordinates,
+                  const size_t begin,
+                  const size_t batchSize) const;
+
+  /*
+   * Evaluate a function with the given coordinates.
+   *
+   * @param coordinates The function coordinates.
+   */
+  double Evaluate(const arma::mat& coordinates) const;
+};
+
+} // namespace test
+} // namespace ens
+
+// Include implementation.
+#include "cross_in_tray_function_impl.hpp"
+
+#endif // ENSMALLEN_PROBLEMS_CROSS_IN_TRAY_FUNCTION_HPP
diff --git a/inst/include/ensmallen_bits/problems/cross_in_tray_function_impl.hpp b/inst/include/ensmallen_bits/problems/cross_in_tray_function_impl.hpp
new file mode 100644
index 0000000..4280f33
--- /dev/null
+++ b/inst/include/ensmallen_bits/problems/cross_in_tray_function_impl.hpp
@@ -0,0 +1,48 @@
+/**
+ * @file cross_in_tray_function_impl.hpp
+ * @author Suryoday Basak
+ *
+ * Implementation of the Cross-in-Tray function.
+ *
+ * ensmallen 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 ensmallen.  If not, see
+ * http://www.opensource.org/licenses/BSD-3-Clause for more information.
+ */
+#ifndef ENSMALLEN_PROBLEMS_CROSS_IN_TRAY_FUNCTION_IMPL_HPP
+#define ENSMALLEN_PROBLEMS_CROSS_IN_TRAY_FUNCTION_IMPL_HPP
+
+// In case it hasn't been included yet.
+#include "cross_in_tray_function.hpp"
+using namespace std;
+
+namespace ens {
+namespace test {
+
+inline CrossInTrayFunction::CrossInTrayFunction() { /* Nothing to do here */ }
+
+inline void CrossInTrayFunction::Shuffle() { /* Nothing to do here */ }
+
+inline double CrossInTrayFunction::Evaluate(const arma::mat& coordinates,
+                                            const size_t /* begin */,
+                                            const size_t /* batchSize */) const
+{
+  // For convenience; we assume these temporaries will be optimized out.
+  const double x1 = coordinates(0);
+  const double x2 = coordinates(1);
+
+  const double objective = -0.0001 * pow(abs(sin(x1) * sin(x2) *
+      exp(abs(100 - (sqrt(pow(x1, 2) + pow(x2, 2)) /
+      arma::datum::pi))) + 1), 0.1);
+  return objective;
+}
+
+inline double CrossInTrayFunction::Evaluate(const arma::mat& coordinates) const
+{
+  return Evaluate(coordinates, 0, NumFunctions());
+}
+
+} // namespace test
+} // namespace ens
+
+#endif
diff --git a/inst/include/ensmallen_bits/problems/drop_wave_function.hpp b/inst/include/ensmallen_bits/problems/drop_wave_function.hpp
index 3c00876..602b29e 100644
--- a/inst/include/ensmallen_bits/problems/drop_wave_function.hpp
+++ b/inst/include/ensmallen_bits/problems/drop_wave_function.hpp
@@ -44,9 +44,9 @@ class DropWaveFunction
   DropWaveFunction();
 
   /**
-  * Shuffle the order of function visitation. This may be called by the
-  * optimizer.
-  */
+   * Shuffle the order of function visitation. This may be called by the
+   * optimizer.
+   */
   void Shuffle();
 
   //! Return 1 (the number of functions).
@@ -55,7 +55,7 @@ class DropWaveFunction
   //! Get the starting point.
   arma::mat GetInitialPoint() const { return arma::mat("0.5; 0.5"); }
 
-  /*
+  /**
    * Evaluate a function for a particular batch-size.
    *
    * @param coordinates The function coordinates.
@@ -66,14 +66,14 @@ class DropWaveFunction
                   const size_t begin,
                   const size_t batchSize) const;
 
-  /*
+  /**
    * Evaluate a function with the given coordinates.
    *
    * @param coordinates The function coordinates.
    */
   double Evaluate(const arma::mat& coordinates) const;
 
-  /*
+  /**
    * Evaluate the gradient of a function for a particular batch-size.
    *
    * @param coordinates The function coordinates.
@@ -86,7 +86,7 @@ class DropWaveFunction
                 arma::mat& gradient,
                 const size_t batchSize) const;
 
-  /*
+  /**
    * Evaluate the gradient of a function with the given coordinates.
    *
    * @param coordinates The function coordinates.
diff --git a/inst/include/ensmallen_bits/problems/easom_function.hpp b/inst/include/ensmallen_bits/problems/easom_function.hpp
index c5c7014..dfafd3d 100644
--- a/inst/include/ensmallen_bits/problems/easom_function.hpp
+++ b/inst/include/ensmallen_bits/problems/easom_function.hpp
@@ -44,9 +44,9 @@ class EasomFunction
   EasomFunction();
 
   /**
-  * Shuffle the order of function visitation. This may be called by the
-  * optimizer.
-  */
+   * Shuffle the order of function visitation. This may be called by the
+   * optimizer.
+   */
   void Shuffle();
 
   //! Return 1 (the number of functions).
@@ -55,7 +55,7 @@ class EasomFunction
   //! Get the starting point.
   arma::mat GetInitialPoint() const { return arma::mat("-90.0; 90.0"); }
 
-  /*
+  /**
    * Evaluate a function for a particular batch-size.
    *
    * @param coordinates The function coordinates.
@@ -66,14 +66,14 @@ class EasomFunction
                   const size_t begin,
                   const size_t batchSize) const;
 
-  /*
+  /**
    * Evaluate a function with the given coordinates.
    *
    * @param coordinates The function coordinates.
    */
   double Evaluate(const arma::mat& coordinates) const;
 
-  /*
+  /**
    * Evaluate the gradient of a function for a particular batch-size.
    *
    * @param coordinates The function coordinates.
@@ -86,7 +86,7 @@ class EasomFunction
                 arma::mat& gradient,
                 const size_t batchSize) const;
 
-  /*
+  /**
    * Evaluate the gradient of a function with the given coordinates.
    *
    * @param coordinates The function coordinates.
diff --git a/inst/include/ensmallen_bits/problems/eggholder_function.hpp b/inst/include/ensmallen_bits/problems/eggholder_function.hpp
index 97e2679..d643554 100644
--- a/inst/include/ensmallen_bits/problems/eggholder_function.hpp
+++ b/inst/include/ensmallen_bits/problems/eggholder_function.hpp
@@ -45,9 +45,9 @@ class EggholderFunction
   EggholderFunction();
 
   /**
-  * Shuffle the order of function visitation. This may be called by the
-  * optimizer.
-  */
+   * Shuffle the order of function visitation. This may be called by the
+   * optimizer.
+   */
   void Shuffle();
 
   //! Return 1 (the number of functions).
@@ -56,7 +56,7 @@ class EggholderFunction
   //! Get the starting point.
   arma::mat GetInitialPoint() const { return arma::mat("-333; -333"); }
 
-  /*
+  /**
    * Evaluate a function for a particular batch-size.
    *
    * @param coordinates The function coordinates.
@@ -67,14 +67,14 @@ class EggholderFunction
                   const size_t begin,
                   const size_t batchSize) const;
 
-  /*
+  /**
    * Evaluate a function with the given coordinates.
    *
    * @param coordinates The function coordinates.
    */
   double Evaluate(const arma::mat& coordinates) const;
 
-  /*
+  /**
    * Evaluate the gradient of a function for a particular batch-size.
    *
    * @param coordinates The function coordinates.
@@ -87,7 +87,7 @@ class EggholderFunction
                 arma::mat& gradient,
                 const size_t batchSize) const;
 
-  /*
+  /**
    * Evaluate the gradient of a function with the given coordinates.
    *
    * @param coordinates The function coordinates.
diff --git a/inst/include/ensmallen_bits/problems/fw_test_function.hpp b/inst/include/ensmallen_bits/problems/fw_test_function.hpp
index fe976f3..0de4cfb 100644
--- a/inst/include/ensmallen_bits/problems/fw_test_function.hpp
+++ b/inst/include/ensmallen_bits/problems/fw_test_function.hpp
@@ -24,8 +24,7 @@ namespace ens {
 class TestFuncFW
 {
  public:
-  TestFuncFW()
-  {/* Nothing to do. */}
+  TestFuncFW() {/* Nothing to do. */}
 
   /**
    * Evaluation of the function.
diff --git a/inst/include/ensmallen_bits/problems/generalized_rosenbrock_function.hpp b/inst/include/ensmallen_bits/problems/generalized_rosenbrock_function.hpp
index 4052f10..d34e2eb 100644
--- a/inst/include/ensmallen_bits/problems/generalized_rosenbrock_function.hpp
+++ b/inst/include/ensmallen_bits/problems/generalized_rosenbrock_function.hpp
@@ -53,10 +53,10 @@ class GeneralizedRosenbrockFunction
   GeneralizedRosenbrockFunction(const size_t n);
 
   /**
-  * Shuffle the order of function visitation. This may be called by the
-  * optimizer.
-  */
-  void Shuffle();
+   * Shuffle the order of function visitation. This may be called by the
+   * optimizer.
+   */
+ void Shuffle();
 
   //! Return 1 (the number of functions).
   size_t NumFunctions() const { return n - 1; }
@@ -64,7 +64,7 @@ class GeneralizedRosenbrockFunction
   //! Get the starting point.
   const arma::mat& GetInitialPoint() const { return initialPoint;}
 
-  /*
+  /**
    * Evaluate a function for a particular batch-size.
    *
    * @param coordinates The function coordinates.
@@ -75,14 +75,14 @@ class GeneralizedRosenbrockFunction
                   const size_t begin,
                   const size_t batchSize = 1) const;
 
-  /*
+  /**
    * Evaluate a function with the given coordinates.
    *
    * @param coordinates The function coordinates.
    */
   double Evaluate(const arma::mat& coordinates) const;
 
-  /*
+  /**
    * Evaluate the gradient of a function for a particular batch-size.
    *
    * @param coordinates The function coordinates.
@@ -95,7 +95,7 @@ class GeneralizedRosenbrockFunction
                 arma::mat& gradient,
                 const size_t batchSize = 1) const;
 
-  /*
+  /**
    * Evaluate the gradient of a function for a particular batch-size.
    *
    * @param coordinates The function coordinates.
@@ -107,7 +107,7 @@ class GeneralizedRosenbrockFunction
                 arma::sp_mat& gradient,
                 const size_t count) const;
 
-  /*
+  /**
    * Evaluate the gradient of a function with the given coordinates.
    *
    * @param coordinates The function coordinates.
diff --git a/inst/include/ensmallen_bits/problems/goldstein_price_function.hpp b/inst/include/ensmallen_bits/problems/goldstein_price_function.hpp
new file mode 100644
index 0000000..1bb43e4
--- /dev/null
+++ b/inst/include/ensmallen_bits/problems/goldstein_price_function.hpp
@@ -0,0 +1,113 @@
+/**
+ * @file goldstein_price_function.hpp
+ * @author Suryoday Basak
+ *
+ * Definition of the Goldstein-Price function.
+ *
+ * ensmallen 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 ensmallen.  If not, see
+ * http://www.opensource.org/licenses/BSD-3-Clause for more information.
+ */
+#ifndef ENSMALLEN_PROBLEMS_GOLDSTEIN_PRICE_FUNCTION_HPP
+#define ENSMALLEN_PROBLEMS_GOLDSTEIN_PRICE_FUNCTION_HPP
+
+namespace ens {
+namespace test {
+
+/**
+ * The Goldstein-Price function, defined by
+ * \f[
+ * f(x_1, x_2) = (1 + (x_1 + x_2 + 1)^2 * (19 - 14 * x_1 + 3 * x_1^2 - 14 *
+ *               x_2 + 6 * x_1 * x_2 + 3 * x_2^2)) *
+ *               (30 + (2 * x_1 - 3 * x_2)^2 * (18 - 32 * x_1 + 12 * x^2 +
+ *               48 * x_2 - 36 * x_1 * x_2 + 27 * x_2^2))
+ * \f]
+ *
+ * This should optimize to f(x) = 3, at x = [0, -1].
+ *
+ * For more information, please refer to:
+ *
+ * @code
+ * @article{Picheny:2013:BKI:2579829.2579986,
+ *   author     = {Picheny, Victor and Wagner, Tobias and Ginsbourger, David},
+ *   title      = {A Benchmark of Kriging-based Infill Criteria for Noisy
+ *                 Optimization},
+ *   journal    = {Struct. Multidiscip. Optim.},
+ *   issue_date = {September 2013},
+ *   volume     = {48},
+ *   number     = {3},
+ *   month      = sep,
+ *   year       = {2013},
+ *   issn       = {1615-147X},
+ *   pages      = {607--626},
+ *   numpages   = {20},
+ *   doi        = {10.1007/s00158-013-0919-4},
+ * }
+ * @endcode
+ */
+class GoldsteinPriceFunction
+{
+ public:
+  //! Initialize the GoldsteinPriceFunction.
+  GoldsteinPriceFunction();
+
+  /**
+   * Shuffle the order of function visitation. This may be called by the
+   * optimizer.
+   */
+  void Shuffle();
+
+  //! Return 1 (the number of functions).
+  size_t NumFunctions() const { return 1; }
+
+  //! Get the starting point.
+  arma::mat GetInitialPoint() const { return arma::mat("-2; 2"); }
+
+  /**
+   * Evaluate a function for a particular batch-size.
+   *
+   * @param coordinates The function coordinates.
+   * @param begin The first function.
+   * @param batchSize Number of points to process.
+   */
+  double Evaluate(const arma::mat& coordinates,
+                  const size_t begin,
+                  const size_t batchSize) const;
+
+  /**
+   * Evaluate a function with the given coordinates.
+   *
+   * @param coordinates The function coordinates.
+   */
+  double Evaluate(const arma::mat& coordinates) const;
+
+  /**
+   * Evaluate the gradient of a function for a particular batch-size.
+   *
+   * @param coordinates The function coordinates.
+   * @param begin The first function.
+   * @param gradient The function gradient.
+   * @param batchSize Number of points to process.
+   */
+  void Gradient(const arma::mat& coordinates,
+                const size_t begin,
+                arma::mat& gradient,
+                const size_t batchSize) const;
+
+  /**
+   * Evaluate the gradient of a function with the given coordinates.
+   *
+   * @param coordinates The function coordinates.
+   * @param gradient The function gradient.
+   */
+  void Gradient(const arma::mat& coordinates, arma::mat& gradient);
+};
+
+} // namespace test
+} // namespace ens
+
+// Include implementation.
+#include "goldstein_price_function_impl.hpp"
+
+#endif // ENSMALLEN_PROBLEMS_GOLDSTEIN_PRICE_FUNCTION_HPP
diff --git a/inst/include/ensmallen_bits/problems/goldstein_price_function_impl.hpp b/inst/include/ensmallen_bits/problems/goldstein_price_function_impl.hpp
new file mode 100644
index 0000000..18158fe
--- /dev/null
+++ b/inst/include/ensmallen_bits/problems/goldstein_price_function_impl.hpp
@@ -0,0 +1,93 @@
+/**
+ * @file goldstein_price_function_impl.hpp
+ * @author Suryoday Basak
+ *
+ * Implementation of the Goldstein-Price function.
+ *
+ * ensmallen 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 ensmallen.  If not, see
+ * http://www.opensource.org/licenses/BSD-3-Clause for more information.
+ */
+#ifndef ENSMALLEN_PROBLEMS_GOLDSTEIN_PRICE_FUNCTION_IMPL_HPP
+#define ENSMALLEN_PROBLEMS_GOLDSTEIN_PRICE_FUNCTION_IMPL_HPP
+
+using namespace std;
+
+// In case it hasn't been included yet.
+#include "goldstein_price_function.hpp"
+
+namespace ens {
+namespace test {
+
+inline GoldsteinPriceFunction::GoldsteinPriceFunction()
+{ /* Nothing to do here */ }
+
+inline void GoldsteinPriceFunction::Shuffle() { /* Nothing to do here */ }
+
+inline double GoldsteinPriceFunction::Evaluate(
+    const arma::mat& coordinates,
+    const size_t /* begin */,
+    const size_t /* batchSize */) const
+{
+  // For convenience; we assume these temporaries will be optimized out.
+  const double x1 = coordinates(0);
+  const double x2 = coordinates(1);
+
+  const double x1Sq = pow(x1, 2);
+  const double x2Sq = pow(x2, 2);
+  const double x1x2 = x1 * x2;
+  const double objective = (1 + pow(x1 + x2 + 1, 2) * (19 - 14 * x1 + 3 * x1Sq -
+      14 * x2 + 6 * x1x2 + 3 * x2Sq)) * (30 + pow(2 * x1 - 3 * x2, 2) *
+      (18 - 32 * x1 + 12 * x1Sq + 48 * x2 - 36 * x1x2 + 27 * x2Sq));
+
+  return objective;
+}
+
+inline double GoldsteinPriceFunction::
+Evaluate(const arma::mat& coordinates) const
+{
+  return Evaluate(coordinates, 0, NumFunctions());
+}
+
+inline void GoldsteinPriceFunction::Gradient(const arma::mat& coordinates,
+                                             const size_t /* begin */,
+                                             arma::mat& gradient,
+                                             const size_t /* batchSize */) const
+{
+  // For convenience; we assume these temporaries will be optimized out.
+  const double x1 = coordinates(0);
+  const double x2 = coordinates(1);
+  gradient.set_size(2, 1);
+
+  gradient(0) = (pow(2 * x1 - 3 * x2, 2) * (24 * x1 - 36 * x2 - 32) +
+      (8 * x1 - 12 * x2) * (12 * x1 * x1 - 36 * x1 * x2 - 32 * x1 +
+      27 * x2 * x2 + 48 * x2 + 18)) * (pow(x1 + x2 + 1, 2) *
+      (3 * x1 * x1 + 6 * x1 * x2 - 14 * x1 + 3 * x2 * x2 - 14 * x2 +
+      19) + 1) + (pow(2 * x1 - 3 * x2, 2) * (12 * x1 * x1 - 36 *
+      x1 * x2 - 32 * x1 + 27 * x2 * x2 + 48 * x2 + 18) + 30) *
+      (pow(x1 + x2 + 1, 2) * (6 * x1 + 6 * x2 - 14) + (2 * x1 +
+      2 * x2 + 2) * (3 * x1 * x1 + 6 * x1 * x2 - 14 * x1 + 3 * x2 *
+      x2 - 14 * x2 + 19));
+
+  gradient(1) = ((- 12 * x1 + 18 * x2) * (12 * x1 * x1 - 36 * x1 * x2 - 32 *
+      x1 + 27 * x2 * x2 + 48 * x2 + 18) + pow(2 * x1 - 3 * x2, 2) *
+      (-36 * x1 + 54 * x2 + 48)) * (pow(x1 + x2 + 1, 2) * (3 * x1 *
+      x1 + 6 * x1 * x2 - 14 * x1 + 3 * x2 * x2 - 14 * x2 + 19) + 1) +
+      (pow(2 * x1 - 3 * x2, 2) * (12 * x1 * x1 - 36 * x1 * x2 -
+      32 * x1 + 27 * x2 * x2 + 48 * x2 + 18) + 30) *
+      (pow(x1 + x2 + 1, 2) * (6 * x1 + 6 * x2 - 14) + (2 * x1 +
+      2 * x2 + 2) * (3 * x1 * x1 + 6 * x1 * x2 - 14 * x1 + 3 * x2 *
+      x2 - 14 * x2 + 19));
+}
+
+inline void GoldsteinPriceFunction::Gradient(const arma::mat& coordinates,
+                                             arma::mat& gradient)
+{
+  Gradient(coordinates, 0, gradient, 1);
+}
+
+} // namespace test
+} // namespace ens
+
+#endif
diff --git a/inst/include/ensmallen_bits/problems/himmelblau_function.hpp b/inst/include/ensmallen_bits/problems/himmelblau_function.hpp
new file mode 100644
index 0000000..d3b63d8
--- /dev/null
+++ b/inst/include/ensmallen_bits/problems/himmelblau_function.hpp
@@ -0,0 +1,108 @@
+/**
+ * @file himmelblau_function.hpp
+ * @author Suryoday Basak
+ *
+ * Definition of the Himmelblau function.
+ *
+ * ensmallen 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 ensmallen.  If not, see
+ * http://www.opensource.org/licenses/BSD-3-Clause for more information.
+ */
+#ifndef ENSMALLEN_PROBLEMS_HIMMELBLAU_FUNCTION_HPP
+#define ENSMALLEN_PROBLEMS_HIMMELBLAU_FUNCTION_HPP
+
+namespace ens {
+namespace test {
+
+/**
+ * The Himmelblau function, defined by
+ *
+ * \f[
+ * f(x_1,x_2) = (x_1^2 + y - 11)^2 + (x_1 + x_2^2 - 7)^2
+ * \f]
+ *
+ * This should optimize to f(x) = 0, at x = [3.0,  2.0], or
+ *          x = [-2.805118, 3.131312], or
+ *          x = [-3.779310, -3.283186], or
+ *          x = [3.584428, -1.848126].
+ *
+ * For more information, please refer to:
+ *
+ * @code
+ * @book{davidmautnerhimmelblau1972,
+ *   Author      = {David Mautner Himmelblau},
+ *   title       = {Applied Nonlinear Programming},
+ *   description = {Applied Nonlinear Programming (Book, 1972)},
+ *   publisher   = {McGraw-Hill},
+ *   year        = {1972},
+ *   month       = {jun},
+ *   isbn        = {0070289212},
+ * }
+ * @endcode
+ */
+class HimmelblauFunction
+{
+ public:
+  //! Initialize the HimmelblauFunction.
+  HimmelblauFunction();
+
+  /**
+   * Shuffle the order of function visitation. This may be called by the
+   * optimizer.
+   */
+  void Shuffle();
+
+  //! Return 1 (the number of functions).
+  size_t NumFunctions() const { return 1; }
+
+  //! Get the starting point.
+  arma::mat GetInitialPoint() const { return arma::mat("5; -5"); }
+
+  /**
+   * Evaluate a function for a particular batch-size.
+   *
+   * @param coordinates The function coordinates.
+   * @param begin The first function.
+   * @param batchSize Number of points to process.
+   */
+  double Evaluate(const arma::mat& coordinates,
+                  const size_t begin,
+                  const size_t batchSize) const;
+
+  /**
+   * Evaluate a function with the given coordinates.
+   *
+   * @param coordinates The function coordinates.
+   */
+  double Evaluate(const arma::mat& coordinates) const;
+
+  /**
+   * Evaluate the gradient of a function for a particular batch-size.
+   *
+   * @param coordinates The function coordinates.
+   * @param begin The first function.
+   * @param gradient The function gradient.
+   * @param batchSize Number of points to process.
+   */
+  void Gradient(const arma::mat& coordinates,
+                const size_t begin,
+                arma::mat& gradient,
+                const size_t batchSize) const;
+
+  /**
+   * Evaluate the gradient of a function with the given coordinates.
+   *
+   * @param coordinates The function coordinates.
+   * @param gradient The function gradient.
+   */
+  void Gradient(const arma::mat& coordinates, arma::mat& gradient);
+};
+
+} // namespace test
+} // namespace ens
+
+// Include implementation.
+#include "himmelblau_function_impl.hpp"
+
+#endif // ENSMALLEN_PROBLEMS_HIMMELBLAU_FUNCTION_HPP
diff --git a/inst/include/ensmallen_bits/problems/himmelblau_function_impl.hpp b/inst/include/ensmallen_bits/problems/himmelblau_function_impl.hpp
new file mode 100644
index 0000000..c2af2df
--- /dev/null
+++ b/inst/include/ensmallen_bits/problems/himmelblau_function_impl.hpp
@@ -0,0 +1,71 @@
+/**
+ * @file himmelblau_function_impl.hpp
+ * @author Suryoday Basak
+ *
+ * Implementation of the Beale function.
+ *
+ * ensmallen 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 ensmallen.  If not, see
+ * http://www.opensource.org/licenses/BSD-3-Clause for more information.
+ */
+#ifndef ENSMALLEN_PROBLEMS_HIMMELBLAU_FUNCTION_IMPL_HPP
+#define ENSMALLEN_PROBLEMS_HIMMELBLAU_FUNCTION_IMPL_HPP
+
+// In case it hasn't been included yet.
+#include "himmelblau_function.hpp"
+using namespace std;
+
+namespace ens {
+namespace test {
+
+inline HimmelblauFunction::HimmelblauFunction() { /* Nothing to do here */ }
+
+inline void HimmelblauFunction::Shuffle() { /* Nothing to do here */ }
+
+inline double HimmelblauFunction::Evaluate(const arma::mat& coordinates,
+                                           const size_t /* begin */,
+                                           const size_t /* batchSize */) const
+{
+  // For convenience; we assume these temporaries will be optimized out.
+  const double x1 = coordinates(0);
+  const double x2 = coordinates(1);
+
+  const double objective = pow(x1 * x1 + x2  - 11 , 2) +
+                           pow(x1 + x2 * x2 - 7, 2);
+  return objective;
+}
+
+inline double HimmelblauFunction::Evaluate(const arma::mat& coordinates) const
+{
+  return Evaluate(coordinates, 0, NumFunctions());
+}
+
+inline void HimmelblauFunction::Gradient(const arma::mat& coordinates,
+                                         const size_t /* begin */,
+                                         arma::mat& gradient,
+                                         const size_t /* batchSize */) const
+{
+  // For convenience; we assume these temporaries will be optimized out.
+  const double x1 = coordinates(0);
+  const double x2 = coordinates(1);
+
+  // Aliases for different terms in the expression of the gradient
+  const double x1Sq = x1 * x1;
+  const double x2Sq = x2 * x2;
+
+  gradient.set_size(2, 1);
+  gradient(0) = (4 * x1 * (x1Sq + x2 - 11)) + (2 * (x1 + x2Sq - 7));
+  gradient(1) = (2 * (x1Sq + x2 - 11)) + (4 * x2 * (x1 + x2Sq - 7));
+}
+
+inline void HimmelblauFunction::Gradient(const arma::mat& coordinates,
+                                         arma::mat& gradient)
+{
+  Gradient(coordinates, 0, gradient, 1);
+}
+
+} // namespace test
+} // namespace ens
+
+#endif
diff --git a/inst/include/ensmallen_bits/problems/holder_table_function.hpp b/inst/include/ensmallen_bits/problems/holder_table_function.hpp
new file mode 100644
index 0000000..8db849c
--- /dev/null
+++ b/inst/include/ensmallen_bits/problems/holder_table_function.hpp
@@ -0,0 +1,88 @@
+/**
+ * @file holder_table_function.hpp
+ * @author Suryoday Basak
+ *
+ * Definition of the Holder table function.
+ *
+ * ensmallen 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 ensmallen.  If not, see
+ * http://www.opensource.org/licenses/BSD-3-Clause for more information.
+ */
+#ifndef ENSMALLEN_PROBLEMS_HOLDER_TABLE_FUNCTION_HPP
+#define ENSMALLEN_PROBLEMS_HOLDER_TABLE_FUNCTION_HPP
+
+namespace ens {
+namespace test {
+
+/**
+ * The Holder table function, defined by
+ *
+ * \f[
+ * f(x1, x2) = - |sin(x1) * cos(x2) * exp(|1 - (sqrt(x1^2 + x2^2) / pi)|)|
+ * \f]
+ *
+ * This should optimize to f(x1, x2) = -19.2085, at
+ *                                     (x1, x2) = [-8.05502, -9.66459], or
+ *                                     (x1, x2) = [8.05502, -9.66459],  or
+ *                                     (x1, x2) = [-8.05502, 9.66459],  or
+ *                                     (x1, x2) = [8.05502, 9.66459]
+ *
+ * For more information, please refer to:
+ *
+ * @code
+ * @article{Mishra2006,
+ *   doi = {10.2139/ssrn.926132},
+ *   year = {2006},
+ *   publisher = {Elsevier {BV}},
+ *   author = {S. K. Mishra},
+ *   title = {Some New Test Functions for Global Optimization and
+ *            Performance of Repulsive Particle Swarm Method},
+ *   journal = {{SSRN} Electronic Journal}
+ * }
+ * @endcode
+ */
+class HolderTableFunction
+{
+ public:
+  //! Initialize the HolderTableFunction.
+  HolderTableFunction();
+
+  /**
+   * Shuffle the order of function visitation. This may be called by the
+   * optimizer.
+   */
+  void Shuffle();
+
+  //! Return 1 (the number of functions).
+  size_t NumFunctions() const { return 1; }
+
+  //! Get the starting point.
+  arma::mat GetInitialPoint() const { return arma::mat("7; 7"); }
+
+  /**
+   * Evaluate a function for a particular batch-size.
+   *
+   * @param coordinates The function coordinates.
+   * @param begin The first function.
+   * @param batchSize Number of points to process.
+   */
+  double Evaluate(const arma::mat& coordinates,
+      const size_t begin,
+      const size_t batchSize) const;
+
+  /**
+   * Evaluate a function with the given coordinates.
+   *
+   * @param coordinates The function coordinates.
+   */
+  double Evaluate(const arma::mat& coordinates) const;
+};
+
+} // namespace test
+} // namespace ens
+
+// Include implementation.
+#include "holder_table_function_impl.hpp"
+
+#endif // ENSMALLEN_PROBLEMS_HOLDER_TABLE_FUNCTION_HPP
diff --git a/inst/include/ensmallen_bits/problems/holder_table_function_impl.hpp b/inst/include/ensmallen_bits/problems/holder_table_function_impl.hpp
new file mode 100644
index 0000000..6fb0ece
--- /dev/null
+++ b/inst/include/ensmallen_bits/problems/holder_table_function_impl.hpp
@@ -0,0 +1,48 @@
+/**
+ * @file holder_table_function_impl.hpp
+ * @author Suryoday Basak
+ *
+ * Implementation of the Holder table function.
+ *
+ * ensmallen 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 ensmallen.  If not, see
+ * http://www.opensource.org/licenses/BSD-3-Clause for more information.
+ */
+#ifndef ENSMALLEN_PROBLEMS_HOLDER_TABLE_FUNCTION_IMPL_HPP
+#define ENSMALLEN_PROBLEMS_HOLDER_TABLE_FUNCTION_IMPL_HPP
+
+// In case it hasn't been included yet.
+#include "holder_table_function.hpp"
+using namespace std;
+
+namespace ens {
+namespace test {
+
+inline HolderTableFunction::HolderTableFunction() { /* Nothing to do here */ }
+
+inline void HolderTableFunction::Shuffle() { /* Nothing to do here */ }
+
+inline double HolderTableFunction::Evaluate(const arma::mat& coordinates,
+                                            const size_t /* begin */,
+                                            const size_t /* batchSize */) const
+{
+  // For convenience; we assume these temporaries will be optimized out.
+  const double x1 = coordinates(0);
+  const double x2 = coordinates(1);
+
+  const double objective = - abs(sin(x1) * cos(x2) * exp(abs(1 -
+      (sqrt(x1 * x1 + x2 * x2) / arma::datum::pi))));
+
+  return objective;
+}
+
+inline double HolderTableFunction::Evaluate(const arma::mat& coordinates) const
+{
+  return Evaluate(coordinates, 0, NumFunctions());
+}
+
+} // namespace test
+} // namespace ens
+
+#endif
diff --git a/inst/include/ensmallen_bits/problems/levy_function_n13.hpp b/inst/include/ensmallen_bits/problems/levy_function_n13.hpp
new file mode 100644
index 0000000..708a45c
--- /dev/null
+++ b/inst/include/ensmallen_bits/problems/levy_function_n13.hpp
@@ -0,0 +1,100 @@
+/**
+ * @file levy_function_n13.hpp
+ * @author Suryoday Basak
+ *
+ * Definition of the Levy function N.13.
+ *
+ * ensmallen 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 ensmallen.  If not, see
+ * http://www.opensource.org/licenses/BSD-3-Clause for more information.
+ */
+#ifndef ENSMALLEN_PROBLEMS_LEVY_FUNCTION_N13_HPP
+#define ENSMALLEN_PROBLEMS_LEVY_FUNCTION_N13_HPP
+
+namespace ens {
+namespace test {
+
+/**
+ * The Levy function N.13, defined by
+ *
+ * \f[
+ * f(x_1, x_2) = sin(3 * pi * x_1)^2 + (x_1 - 1)^2 *
+ *               (1 + (sin(3 * pi * x_2)^2)) +
+ *               (x_2 - 1)^2 * (1 + sin(2 * pi * x_2)^2)
+ * \f]
+ * This should optimize to f(x_1, x_2) = 0, at x = [1, 1].
+ *
+ * For more information, please refer to:
+ *
+ * @code
+ * @misc{LevyFunction,
+ *   URL = {http://www.sfu.ca/~ssurjano/levy13.html},
+ * }
+ * @endcode
+ */
+class LevyFunctionN13
+{
+ public:
+  //! Initialize the BealeFunction.
+  LevyFunctionN13();
+
+  /**
+   * Shuffle the order of function visitation. This may be called by the
+   * optimizer.
+   */
+  void Shuffle();
+
+  //! Return 1 (the number of functions).
+  size_t NumFunctions() const { return 1; }
+
+  //! Get the starting point.
+  arma::mat GetInitialPoint() const { return arma::mat("-10; 10"); }
+
+  /**
+   * Evaluate a function for a particular batch-size.
+   *
+   * @param coordinates The function coordinates.
+   * @param begin The first function.
+   * @param batchSize Number of points to process.
+   */
+  double Evaluate(const arma::mat& coordinates,
+                  const size_t begin,
+                  const size_t batchSize) const;
+
+  /**
+   * Evaluate a function with the given coordinates.
+   *
+   * @param coordinates The function coordinates.
+   */
+  double Evaluate(const arma::mat& coordinates) const;
+
+  /**
+   * Evaluate the gradient of a function for a particular batch-size.
+   *
+   * @param coordinates The function coordinates.
+   * @param begin The first function.
+   * @param gradient The function gradient.
+   * @param batchSize Number of points to process.
+   */
+  void Gradient(const arma::mat& coordinates,
+                const size_t begin,
+                arma::mat& gradient,
+                const size_t batchSize) const;
+
+  /**
+   * Evaluate the gradient of a function with the given coordinates.
+   *
+   * @param coordinates The function coordinates.
+   * @param gradient The function gradient.
+   */
+  void Gradient(const arma::mat& coordinates, arma::mat& gradient);
+};
+
+} // namespace test
+} // namespace ens
+
+// Include implementation.
+#include "levy_function_n13_impl.hpp"
+
+#endif // ENSMALLEN_PROBLEMS_LEVY_FUNCTION_N13_HPP
diff --git a/inst/include/ensmallen_bits/problems/levy_function_n13_impl.hpp b/inst/include/ensmallen_bits/problems/levy_function_n13_impl.hpp
new file mode 100644
index 0000000..6e57fc0
--- /dev/null
+++ b/inst/include/ensmallen_bits/problems/levy_function_n13_impl.hpp
@@ -0,0 +1,76 @@
+/**
+ * @file levy_function_n13_impl.hpp
+ * @author Suryoday Basak
+ *
+ * Implementation of the Levy function N.13.
+ *
+ * ensmallen 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 ensmallen.  If not, see
+ * http://www.opensource.org/licenses/BSD-3-Clause for more information.
+ */
+#ifndef ENSMALLEN_PROBLEMS_LEVY_FUNCTION_N13_IMPL_HPP
+#define ENSMALLEN_PROBLEMS_LEVY_FUNCTION_N13_IMPL_HPP
+
+// In case it hasn't been included yet.
+#include "levy_function_n13.hpp"
+using namespace std;
+
+namespace ens {
+namespace test {
+
+inline LevyFunctionN13::LevyFunctionN13() { /* Nothing to do here */ }
+
+inline void LevyFunctionN13::Shuffle() { /* Nothing to do here */ }
+
+inline double LevyFunctionN13::Evaluate(const arma::mat& coordinates,
+                                        const size_t /* begin */,
+                                        const size_t /* batchSize */) const
+{
+  // For convenience; we assume these temporaries will be optimized out.
+  const double x1 = coordinates(0);
+  const double x2 = coordinates(1);
+
+  const double objective = pow(sin(3 * arma::datum::pi * x1), 2) +
+      (pow(x1 - 1, 2) * (1 + pow(sin(3 * arma::datum::pi * x2), 2))) +
+      (pow(x2 - 1, 2) * (1 + pow(sin(2 * arma::datum::pi * x2), 2)));
+
+  return objective;
+}
+
+inline double LevyFunctionN13::Evaluate(const arma::mat& coordinates) const
+{
+  return Evaluate(coordinates, 0, NumFunctions());
+}
+
+inline void LevyFunctionN13::Gradient(const arma::mat& coordinates,
+                                      const size_t /* begin */,
+                                      arma::mat& gradient,
+                                      const size_t /* batchSize */) const
+{
+  // For convenience; we assume these temporaries will be optimized out.
+  const double x1 = coordinates(0);
+  const double x2 = coordinates(1);
+  gradient.set_size(2, 1);
+
+  gradient(0) = (2 * x1 - 2) * (pow(sin(3 * arma::datum::pi * x2), 2) + 1) +
+    6 * arma::datum::pi * sin(3 * arma::datum::pi * x1) *
+    cos(3 * arma::datum::pi * x1);
+
+  gradient(1) = 6 * arma::datum::pi * pow(x1 - 1, 2) * sin(3 *
+      arma::datum::pi * x2) * cos(3 * arma::datum::pi * x2) +
+      4 * arma::datum::pi * pow(x2 - 1, 2) * sin(2 *
+      arma::datum::pi * x2) * cos(2 * arma::datum::pi * x2) +
+      (2 * x2 - 2) * (pow(sin(2 * arma::datum::pi * x2), 2) + 1);
+}
+
+inline void LevyFunctionN13::Gradient(const arma::mat& coordinates,
+                                      arma::mat& gradient)
+{
+  Gradient(coordinates, 0, gradient, 1);
+}
+
+} // namespace test
+} // namespace ens
+
+#endif
diff --git a/inst/include/ensmallen_bits/problems/logistic_regression_function.hpp b/inst/include/ensmallen_bits/problems/logistic_regression_function.hpp
index cd5ae09..6621a6c 100644
--- a/inst/include/ensmallen_bits/problems/logistic_regression_function.hpp
+++ b/inst/include/ensmallen_bits/problems/logistic_regression_function.hpp
@@ -51,8 +51,8 @@ class LogisticRegressionFunction
   const arma::Row<size_t>& Responses() const { return responses; }
 
   /**
-  * Shuffle the order of function visitation.  This may be called by the optimizer.
-  */
+   * Shuffle the order of function visitation.  This may be called by the optimizer.
+   */
   void Shuffle();
 
   /**
diff --git a/inst/include/ensmallen_bits/problems/matyas_function.hpp b/inst/include/ensmallen_bits/problems/matyas_function.hpp
index 9675956..2f47e7d 100644
--- a/inst/include/ensmallen_bits/problems/matyas_function.hpp
+++ b/inst/include/ensmallen_bits/problems/matyas_function.hpp
@@ -44,9 +44,9 @@ class MatyasFunction
   MatyasFunction();
 
   /**
-  * Shuffle the order of function visitation. This may be called by the
-  * optimizer.
-  */
+   * Shuffle the order of function visitation. This may be called by the
+   * optimizer.
+   */
   void Shuffle();
 
   //! Return 1 (the number of functions).
@@ -55,7 +55,7 @@ class MatyasFunction
   //! Get the starting point.
   arma::mat GetInitialPoint() const { return arma::mat("-3; 3"); }
 
-  /*
+  /**
    * Evaluate a function for a particular batch-size.
    *
    * @param coordinates The function coordinates.
@@ -66,14 +66,14 @@ class MatyasFunction
                   const size_t begin,
                   const size_t batchSize) const;
 
-  /*
+  /**
    * Evaluate a function with the given coordinates.
    *
    * @param coordinates The function coordinates.
    */
   double Evaluate(const arma::mat& coordinates) const;
 
-  /*
+  /**
    * Evaluate the gradient of a function for a particular batch-size.
    *
    * @param coordinates The function coordinates.
@@ -86,7 +86,7 @@ class MatyasFunction
                 arma::mat& gradient,
                 const size_t batchSize) const;
 
-  /*
+  /**
    * Evaluate the gradient of a function with the given coordinates.
    *
    * @param coordinates The function coordinates.
diff --git a/inst/include/ensmallen_bits/problems/matyas_function_impl.hpp b/inst/include/ensmallen_bits/problems/matyas_function_impl.hpp
index 1e60f01..525184d 100644
--- a/inst/include/ensmallen_bits/problems/matyas_function_impl.hpp
+++ b/inst/include/ensmallen_bits/problems/matyas_function_impl.hpp
@@ -31,7 +31,7 @@ inline double MatyasFunction::Evaluate(const arma::mat& coordinates,
   const double x2 = coordinates(1);
 
   const double objective = 0.26 * (pow(x1, 2) + std::pow(x2, 2)) -
-    0.48 * x1 * x2;
+      0.48 * x1 * x2;
 
   return objective;
 }
diff --git a/inst/include/ensmallen_bits/problems/mc_cormick_function.hpp b/inst/include/ensmallen_bits/problems/mc_cormick_function.hpp
index edf0146..b3d4a9a 100644
--- a/inst/include/ensmallen_bits/problems/mc_cormick_function.hpp
+++ b/inst/include/ensmallen_bits/problems/mc_cormick_function.hpp
@@ -44,9 +44,9 @@ class McCormickFunction
   McCormickFunction();
 
   /**
-  * Shuffle the order of function visitation. This may be called by the
-  * optimizer.
-  */
+   * Shuffle the order of function visitation. This may be called by the
+   * optimizer.
+   */
   void Shuffle();
 
   //! Return 1 (the number of functions).
@@ -55,7 +55,7 @@ class McCormickFunction
   //! Get the starting point.
   arma::mat GetInitialPoint() const { return arma::mat("-2; 4"); }
 
-  /*
+  /**
    * Evaluate a function for a particular batch-size.
    *
    * @param coordinates The function coordinates.
@@ -66,14 +66,14 @@ class McCormickFunction
                   const size_t begin,
                   const size_t batchSize) const;
 
-  /*
+  /**
    * Evaluate a function with the given coordinates.
    *
    * @param coordinates The function coordinates.
    */
   double Evaluate(const arma::mat& coordinates) const;
 
-  /*
+  /**
    * Evaluate the gradient of a function for a particular batch-size.
    *
    * @param coordinates The function coordinates.
@@ -86,7 +86,7 @@ class McCormickFunction
                 arma::mat& gradient,
                 const size_t batchSize) const;
 
-  /*
+  /**
    * Evaluate the gradient of a function with the given coordinates.
    *
    * @param coordinates The function coordinates.
diff --git a/inst/include/ensmallen_bits/problems/problems.hpp b/inst/include/ensmallen_bits/problems/problems.hpp
index 88f9e58..9504a38 100644
--- a/inst/include/ensmallen_bits/problems/problems.hpp
+++ b/inst/include/ensmallen_bits/problems/problems.hpp
@@ -7,28 +7,38 @@
 #ifndef ENSMALLEN_PROBLEMS_PROBLEMS_HPP
 #define ENSMALLEN_PROBLEMS_PROBLEMS_HPP
 
+#include "ackley_function.hpp"
 #include "aug_lagrangian_test_functions.hpp"
+#include "beale_function.hpp"
 #include "booth_function.hpp"
 #include "bukin_function.hpp"
 #include "colville_function.hpp"
+#include "cross_in_tray_function.hpp"
 #include "drop_wave_function.hpp"
 #include "easom_function.hpp"
 #include "eggholder_function.hpp"
 #include "fw_test_function.hpp"
 #include "generalized_rosenbrock_function.hpp"
+#include "goldstein_price_function.hpp"
 #include "gradient_descent_test_function.hpp"
+#include "himmelblau_function.hpp"
+#include "holder_table_function.hpp"
+#include "levy_function_n13.hpp"
 #include "logistic_regression_function.hpp"
 #include "matyas_function.hpp"
 #include "mc_cormick_function.hpp"
 #include "rastrigin_function.hpp"
 #include "rosenbrock_function.hpp"
 #include "rosenbrock_wood_function.hpp"
+#include "schaffer_function_n2.hpp"
+#include "schaffer_function_n4.hpp"
 #include "schwefel_function.hpp"
 #include "sgd_test_function.hpp"
 #include "softmax_regression_function.hpp"
 #include "sparse_test_function.hpp"
 #include "sphere_function.hpp"
 #include "styblinski_tang_function.hpp"
+#include "three_hump_camel_function.hpp"
 #include "wood_function.hpp"
 
 #endif
diff --git a/inst/include/ensmallen_bits/problems/rastrigin_function.hpp b/inst/include/ensmallen_bits/problems/rastrigin_function.hpp
index 6d78459..d0b5a58 100644
--- a/inst/include/ensmallen_bits/problems/rastrigin_function.hpp
+++ b/inst/include/ensmallen_bits/problems/rastrigin_function.hpp
@@ -47,10 +47,10 @@ class RastriginFunction
   RastriginFunction(const size_t n);
 
   /**
-  * Shuffle the order of function visitation. This may be called by the
-  * optimizer.
-  */
-  void Shuffle();
+   * Shuffle the order of function visitation. This may be called by the
+   * optimizer.
+   */
+ void Shuffle();
 
   //! Return 1 (the number of functions).
   size_t NumFunctions() const { return n; }
@@ -58,7 +58,7 @@ class RastriginFunction
   //! Get the starting point.
   arma::mat GetInitialPoint() const { return initialPoint; }
 
-  /*
+  /**
    * Evaluate a function for a particular batch-size.
    *
    * @param coordinates The function coordinates.
@@ -69,14 +69,14 @@ class RastriginFunction
                   const size_t begin,
                   const size_t batchSize) const;
 
-  /*
+  /**
    * Evaluate a function with the given coordinates.
    *
    * @param coordinates The function coordinates.
    */
   double Evaluate(const arma::mat& coordinates) const;
 
-  /*
+  /**
    * Evaluate the gradient of a function for a particular batch-size.
    *
    * @param coordinates The function coordinates.
@@ -89,7 +89,7 @@ class RastriginFunction
                 arma::mat& gradient,
                 const size_t batchSize) const;
 
-  /*
+  /**
    * Evaluate the gradient of a function with the given coordinates.
    *
    * @param coordinates The function coordinates.
diff --git a/inst/include/ensmallen_bits/problems/rosenbrock_function.hpp b/inst/include/ensmallen_bits/problems/rosenbrock_function.hpp
index 12dfc8d..a294384 100644
--- a/inst/include/ensmallen_bits/problems/rosenbrock_function.hpp
+++ b/inst/include/ensmallen_bits/problems/rosenbrock_function.hpp
@@ -47,9 +47,9 @@ class RosenbrockFunction
   RosenbrockFunction();
 
   /**
-  * Shuffle the order of function visitation. This may be called by the
-  * optimizer.
-  */
+   * Shuffle the order of function visitation. This may be called by the
+   * optimizer.
+   */
   void Shuffle();
 
   //! Return 1 (the number of functions).
@@ -58,7 +58,7 @@ class RosenbrockFunction
   //! Get the starting point.
   arma::mat GetInitialPoint() const { return arma::mat("-1.2; 1"); }
 
-  /*
+  /**
    * Evaluate a function for a particular batch-size.
    *
    * @param coordinates The function coordinates.
@@ -69,14 +69,14 @@ class RosenbrockFunction
                   const size_t begin,
                   const size_t batchSize) const;
 
-  /*
+  /**
    * Evaluate a function with the given coordinates.
    *
    * @param coordinates The function coordinates.
    */
   double Evaluate(const arma::mat& coordinates) const;
 
-  /*
+  /**
    * Evaluate the gradient of a function for a particular batch-size.
    *
    * @param coordinates The function coordinates.
@@ -89,7 +89,7 @@ class RosenbrockFunction
                 arma::mat& gradient,
                 const size_t batchSize) const;
 
-  /*
+  /**
    * Evaluate the gradient of a function with the given coordinates.
    *
    * @param coordinates The function coordinates.
diff --git a/inst/include/ensmallen_bits/problems/rosenbrock_wood_function.hpp b/inst/include/ensmallen_bits/problems/rosenbrock_wood_function.hpp
index e61429c..a07070b 100644
--- a/inst/include/ensmallen_bits/problems/rosenbrock_wood_function.hpp
+++ b/inst/include/ensmallen_bits/problems/rosenbrock_wood_function.hpp
@@ -31,9 +31,9 @@ class RosenbrockWoodFunction
   RosenbrockWoodFunction();
 
   /**
-  * Shuffle the order of function visitation. This may be called by the
-  * optimizer.
-  */
+   * Shuffle the order of function visitation. This may be called by the
+   * optimizer.
+   */
   void Shuffle();
 
   //! Return 1 (the number of functions).
@@ -42,7 +42,7 @@ class RosenbrockWoodFunction
   //! Get the starting point.
   const arma::mat& GetInitialPoint() const { return initialPoint; }
 
-  /*
+  /**
    * Evaluate a function for a particular batch-size.
    *
    * @param coordinates The function coordinates.
@@ -53,14 +53,14 @@ class RosenbrockWoodFunction
                   const size_t begin,
                   const size_t batchSize) const;
 
-  /*
+  /**
    * Evaluate a function with the given coordinates.
    *
    * @param coordinates The function coordinates.
    */
   double Evaluate(const arma::mat& coordinates) const;
 
-  /*
+  /**
    * Evaluate the gradient of a function for a particular batch-size.
    *
    * @param coordinates The function coordinates.
@@ -73,7 +73,7 @@ class RosenbrockWoodFunction
                 arma::mat& gradient,
                 const size_t batchSize) const;
 
-  /*
+  /**
    * Evaluate the gradient of a function with the given coordinates.
    *
    * @param coordinates The function coordinates.
diff --git a/inst/include/ensmallen_bits/problems/schaffer_function_n2.hpp b/inst/include/ensmallen_bits/problems/schaffer_function_n2.hpp
new file mode 100644
index 0000000..0ef8616
--- /dev/null
+++ b/inst/include/ensmallen_bits/problems/schaffer_function_n2.hpp
@@ -0,0 +1,105 @@
+/**
+ * @file schaffer_function_n2.hpp
+ * @author Suryoday Basak
+ *
+ * Definition of Schaffer function N.2.
+ *
+ * ensmallen 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 ensmallen.  If not, see
+ * http://www.opensource.org/licenses/BSD-3-Clause for more information.
+ */
+#ifndef ENSMALLEN_PROBLEMS_SCHAFFER_FUNCTION_N2_HPP
+#define ENSMALLEN_PROBLEMS_SCHAFFER_FUNCTION_N2_HPP
+
+namespace ens {
+namespace test {
+
+/**
+ * The Schaffer function N.2, defined by
+ *
+ * \f[
+ * f(x1, x2) = 0.5 + ((sin^2(x1^2 - x2^2) - 0.5) /
+ *             (1 + 0.001 * (x1^2 + x2^2))^2)
+ * \f]
+ *
+ * This should optimize to f(x1, x2) = 0, at (x1, x2) = [0, 0].
+ *
+ * For more information, please refer to:
+ *
+ * @code
+ * @article{1308.4008,
+ *   Author = {Momin Jamil and Xin-She Yang},
+ *   Title  = {A Literature Survey of Benchmark Functions For Global
+ *             Optimization Problems},
+ *   Year   = {2013},
+ *   Eprint = {arXiv:1308.4008},
+ *   Doi    = {10.1504/IJMMNO.2013.055204},
+ * }
+ * @endcode
+ */
+class SchafferFunctionN2
+{
+ public:
+  //! Initialize the SchafferFunctionN4.
+  SchafferFunctionN2();
+
+  /**
+   * Shuffle the order of function visitation. This may be called by the
+   * optimizer.
+   */
+  void Shuffle();
+
+  //! Return 1 (the number of functions).
+  size_t NumFunctions() const { return 1; }
+
+  //! Get the starting point.
+  arma::mat GetInitialPoint() const { return arma::mat("-100; 100"); }
+
+  /**
+   * Evaluate a function for a particular batch-size.
+   *
+   * @param coordinates The function coordinates.
+   * @param begin The first function.
+   * @param batchSize Number of points to process.
+   */
+  double Evaluate(const arma::mat& coordinates,
+                  const size_t begin,
+                  const size_t batchSize) const;
+
+  /**
+   * Evaluate a function with the given coordinates.
+   *
+   * @param coordinates The function coordinates.
+   */
+  double Evaluate(const arma::mat& coordinates) const;
+
+  /**
+   * Evaluate the gradient of a function for a particular batch-size.
+   *
+   * @param coordinates The function coordinates.
+   * @param begin The first function.
+   * @param gradient The function gradient.
+   * @param batchSize Number of points to process.
+   */
+  void Gradient(const arma::mat& coordinates,
+                const size_t begin,
+                arma::mat& gradient,
+                const size_t batchSize) const;
+
+  /**
+   * Evaluate the gradient of a function with the given coordinates.
+   *
+   * @param coordinates The function coordinates.
+   * @param gradient The function gradient.
+   */
+  void Gradient(const arma::mat& coordinates, arma::mat& gradient);
+};
+
+} // namespace test
+} // namespace ens
+
+// Include implementation.
+#include "schaffer_function_n2_impl.hpp"
+
+#endif // ENSMALLEN_PROBLEMS_SCHAFFER_FUNCTION_N2_HPP
diff --git a/inst/include/ensmallen_bits/problems/schaffer_function_n2_impl.hpp b/inst/include/ensmallen_bits/problems/schaffer_function_n2_impl.hpp
new file mode 100644
index 0000000..c63ae2f
--- /dev/null
+++ b/inst/include/ensmallen_bits/problems/schaffer_function_n2_impl.hpp
@@ -0,0 +1,79 @@
+/**
+ * @file schaffer_function_n2_impl.hpp
+ * @author Suryoday Basak
+ *
+ * Implementation of Schaffer function N.2.
+ *
+ * ensmallen 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 ensmallen.  If not, see
+ * http://www.opensource.org/licenses/BSD-3-Clause for more information.
+ */
+#ifndef ENSMALLEN_PROBLEMS_SCHAFFER_FUNCTION_N2_IMPL_HPP
+#define ENSMALLEN_PROBLEMS_SCHAFFER_FUNCTION_N2_IMPL_HPP
+
+// In case it hasn't been included yet.
+#include "schaffer_function_n2.hpp"
+using namespace std;
+
+namespace ens {
+namespace test {
+
+inline SchafferFunctionN2::SchafferFunctionN2() { /* Nothing to do here */ }
+
+inline void SchafferFunctionN2::Shuffle() { /* Nothing to do here */ }
+
+inline double SchafferFunctionN2::Evaluate(const arma::mat& coordinates,
+                                           const size_t /* begin */,
+                                           const size_t /* batchSize */) const
+{
+  // For convenience; we assume these temporaries will be optimized out.
+  const double x1 = coordinates(0);
+  const double x2 = coordinates(1);
+
+  const double objective = 0.5 + (pow(sin(pow(x1, 2) - pow(x2, 2)), 2) - 0.5) /
+      pow(1 + 0.001 * (pow(x1, 2) + pow(x2, 2)), 2);
+
+  return objective;
+}
+
+inline double SchafferFunctionN2::Evaluate(const arma::mat& coordinates) const
+{
+  return Evaluate(coordinates, 0, NumFunctions());
+}
+
+inline void SchafferFunctionN2::Gradient(const arma::mat& coordinates,
+                                         const size_t /* begin */,
+                                         arma::mat& gradient,
+                                         const size_t /* batchSize */) const
+{
+  // For convenience; we assume these temporaries will be optimized out.
+  const double x1 = coordinates(0);
+  const double x2 = coordinates(1);
+
+  // Aliases for different terms in the expression of the gradient
+  const double x1Sq = x1 * x1;
+  const double x2Sq = x2 * x2;
+  const double sum1 = x1Sq - x2Sq;
+  const double sinSum1 = sin(sum1);
+  const double sum2 = 0.001 * (x1Sq + x2Sq) + 1;
+  const double trigExpression = 4 * sinSum1 * cos(sum1);
+  const double numerator1 = - 0.004 * (pow(sinSum1, 2) - 0.5);
+  const double expr1 = numerator1 / pow(sum2, 3);
+  const double expr2 = trigExpression / pow(sum2, 2);
+
+  gradient.set_size(2, 1);
+  gradient(0) = x1 * (expr1 + expr2);
+  gradient(1) = x2 * (expr1 - expr2);
+}
+
+inline void SchafferFunctionN2::Gradient(const arma::mat& coordinates,
+                                         arma::mat& gradient)
+{
+  Gradient(coordinates, 0, gradient, 1);
+}
+
+} // namespace test
+} // namespace ens
+
+#endif
diff --git a/inst/include/ensmallen_bits/problems/schaffer_function_n4.hpp b/inst/include/ensmallen_bits/problems/schaffer_function_n4.hpp
new file mode 100644
index 0000000..a91009a
--- /dev/null
+++ b/inst/include/ensmallen_bits/problems/schaffer_function_n4.hpp
@@ -0,0 +1,80 @@
+/**
+ * @file schaffer_function_n4.hpp
+ * @author Suryoday Basak
+ *
+ * Definition of Schaffer function N.4.
+ *
+ * ensmallen 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 ensmallen.  If not, see
+ * http://www.opensource.org/licenses/BSD-3-Clause for more information.
+ */
+#ifndef ENSMALLEN_PROBLEMS_SCHAFFER_FUNCTION_N4_HPP
+#define ENSMALLEN_PROBLEMS_SCHAFFER_FUNCTION_N4_HPP
+
+namespace ens {
+namespace test {
+
+/**
+ * The Schaffer function N.4, defined by
+ *
+ * \f[
+ * f(x1, x2) = 0.5 + (cos^2(sin(|x1^2 - x2^2|)) - 0.5) /
+ *     (1 + 0.001 * (x1^2 + x2^2))^2
+ * \f]
+ *
+ * This should optimize to f(x1, x2) = 0.292579, at (x1, x2) = [0, 1.25313], or
+ *                                                  (x1, x2) = [0, -1.25313].
+ *
+ * For more information, please refer to:
+ *
+ * @code
+ * @misc{LevyFunction,
+ *   URL = {http://benchmarkfcns.xyz/benchmarkfcns/schaffern4fcn.html}
+ * }
+ * @endcode
+ */
+class SchafferFunctionN4
+{
+ public:
+  //! Initialize the SchafferFunctionN4.
+  SchafferFunctionN4();
+
+  /**
+   * Shuffle the order of function visitation. This may be called by the
+   * optimizer.
+   */
+  void Shuffle();
+
+  //! Return 1 (the number of functions).
+  size_t NumFunctions() const { return 1; }
+
+  //! Get the starting point.
+  arma::mat GetInitialPoint() const { return arma::mat("-5; 5"); }
+
+  /**
+   * Evaluate a function for a particular batch-size.
+   *
+   * @param coordinates The function coordinates.
+   * @param begin The first function.
+   * @param batchSize Number of points to process.
+   */
+  double Evaluate(const arma::mat& coordinates,
+                  const size_t begin,
+                  const size_t batchSize) const;
+
+  /**
+   * Evaluate a function with the given coordinates.
+   *
+   * @param coordinates The function coordinates.
+   */
+  double Evaluate(const arma::mat& coordinates) const;
+};
+
+} // namespace test
+} // namespace ens
+
+// Include implementation.
+#include "schaffer_function_n4_impl.hpp"
+
+#endif // ENSMALLEN_PROBLEMS_SCHAFFER_FUNCTION_N4_HPP
diff --git a/inst/include/ensmallen_bits/problems/schaffer_function_n4_impl.hpp b/inst/include/ensmallen_bits/problems/schaffer_function_n4_impl.hpp
new file mode 100644
index 0000000..d614b26
--- /dev/null
+++ b/inst/include/ensmallen_bits/problems/schaffer_function_n4_impl.hpp
@@ -0,0 +1,49 @@
+/**
+ * @file schaffer_function_n4_impl.hpp
+ * @author Suryoday Basak
+ *
+ * Implementation of Schaffer function N.4.
+ *
+ * ensmallen 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 ensmallen.  If not, see
+ * http://www.opensource.org/licenses/BSD-3-Clause for more information.
+ */
+#ifndef ENSMALLEN_PROBLEMS_SCHAFFER_FUNCTION_N4_IMPL_HPP
+#define ENSMALLEN_PROBLEMS_SCHAFFER_FUNCTION_N4_IMPL_HPP
+
+// In case it hasn't been included yet.
+#include "schaffer_function_n4.hpp"
+using namespace std;
+
+namespace ens {
+namespace test {
+
+inline SchafferFunctionN4::SchafferFunctionN4() { /* Nothing to do here */ }
+
+inline void SchafferFunctionN4::Shuffle() { /* Nothing to do here */ }
+
+inline double SchafferFunctionN4::Evaluate(const arma::mat& coordinates,
+                                           const size_t /* begin */,
+                                           const size_t /* batchSize */) const
+{
+  // For convenience; we assume these temporaries will be optimized out.
+  const double x1 = coordinates(0);
+  const double x2 = coordinates(1);
+
+  const double objective = 0.5 + (pow(cos(sin(abs(pow(x1, 2) -
+      pow(x2, 2)))), 2) - 0.5) / pow(1 + 0.001 * (pow(x1, 2) +
+      pow(x2, 2)), 2);
+
+  return objective;
+}
+
+inline double SchafferFunctionN4::Evaluate(const arma::mat& coordinates) const
+{
+  return Evaluate(coordinates, 0, NumFunctions());
+}
+
+} // namespace test
+} // namespace ens
+
+#endif
diff --git a/inst/include/ensmallen_bits/problems/schwefel_function.hpp b/inst/include/ensmallen_bits/problems/schwefel_function.hpp
index 22da47f..f4068be 100644
--- a/inst/include/ensmallen_bits/problems/schwefel_function.hpp
+++ b/inst/include/ensmallen_bits/problems/schwefel_function.hpp
@@ -47,9 +47,9 @@ class SchwefelFunction
   SchwefelFunction(const size_t n);
 
   /**
-  * Shuffle the order of function visitation. This may be called by the
-  * optimizer.
-  */
+   * Shuffle the order of function visitation. This may be called by the
+   * optimizer.
+   */
   void Shuffle();
 
   //! Return 1 (the number of functions).
@@ -58,7 +58,7 @@ class SchwefelFunction
   //! Get the starting point.
   arma::mat GetInitialPoint() const { return initialPoint; }
 
-  /*
+  /**
    * Evaluate a function for a particular batch-size.
    *
    * @param coordinates The function coordinates.
@@ -69,14 +69,14 @@ class SchwefelFunction
                   const size_t begin,
                   const size_t batchSize) const;
 
-  /*
+  /**
    * Evaluate a function with the given coordinates.
    *
    * @param coordinates The function coordinates.
    */
   double Evaluate(const arma::mat& coordinates) const;
 
-  /*
+  /**
    * Evaluate the gradient of a function for a particular batch-size.
    *
    * @param coordinates The function coordinates.
@@ -89,7 +89,7 @@ class SchwefelFunction
                 arma::mat& gradient,
                 const size_t batchSize) const;
 
-  /*
+  /**
    * Evaluate the gradient of a function with the given coordinates.
    *
    * @param coordinates The function coordinates.
diff --git a/inst/include/ensmallen_bits/problems/sgd_test_function.hpp b/inst/include/ensmallen_bits/problems/sgd_test_function.hpp
index 15bf5bd..ae528dc 100644
--- a/inst/include/ensmallen_bits/problems/sgd_test_function.hpp
+++ b/inst/include/ensmallen_bits/problems/sgd_test_function.hpp
@@ -29,8 +29,9 @@ class SGDTestFunction
   SGDTestFunction();
 
   /**
-  * Shuffle the order of function visitation.  This may be called by the optimizer.
-  */
+   * Shuffle the order of function visitation. This may be called by the
+   * optimizer.
+   */
   void Shuffle();
 
   //! Return 3 (the number of functions).
diff --git a/inst/include/ensmallen_bits/problems/sphere_function.hpp b/inst/include/ensmallen_bits/problems/sphere_function.hpp
index f7b6c33..03cb23f 100644
--- a/inst/include/ensmallen_bits/problems/sphere_function.hpp
+++ b/inst/include/ensmallen_bits/problems/sphere_function.hpp
@@ -48,9 +48,9 @@ class SphereFunction
   SphereFunction(const size_t n);
 
   /**
-  * Shuffle the order of function visitation. This may be called by the
-  * optimizer.
-  */
+   * Shuffle the order of function visitation. This may be called by the
+   * optimizer.
+   */
   void Shuffle();
 
   //! Return 1 (the number of functions).
@@ -59,7 +59,7 @@ class SphereFunction
   //! Get the starting point.
   arma::mat GetInitialPoint() const { return initialPoint; }
 
-  /*
+  /**
    * Evaluate a function for a particular batch-size.
    *
    * @param coordinates The function coordinates.
@@ -70,14 +70,14 @@ class SphereFunction
                   const size_t begin,
                   const size_t batchSize) const;
 
-  /*
+  /**
    * Evaluate a function with the given coordinates.
    *
    * @param coordinates The function coordinates.
    */
   double Evaluate(const arma::mat& coordinates) const;
 
-  /*
+  /**
    * Evaluate the gradient of a function for a particular batch-size.
    *
    * @param coordinates The function coordinates.
@@ -90,7 +90,7 @@ class SphereFunction
                 arma::mat& gradient,
                 const size_t batchSize) const;
 
-  /*
+  /**
    * Evaluate the gradient of a function with the given coordinates.
    *
    * @param coordinates The function coordinates.
diff --git a/inst/include/ensmallen_bits/problems/styblinski_tang_function.hpp b/inst/include/ensmallen_bits/problems/styblinski_tang_function.hpp
index d662d95..3e1fc33 100644
--- a/inst/include/ensmallen_bits/problems/styblinski_tang_function.hpp
+++ b/inst/include/ensmallen_bits/problems/styblinski_tang_function.hpp
@@ -49,9 +49,9 @@ class StyblinskiTangFunction
   StyblinskiTangFunction(const size_t n);
 
   /**
-  * Shuffle the order of function visitation. This may be called by the
-  * optimizer.
-  */
+   * Shuffle the order of function visitation. This may be called by the
+   * optimizer.
+   */
   void Shuffle();
 
   //! Return 1 (the number of functions).
@@ -60,7 +60,7 @@ class StyblinskiTangFunction
   //! Get the starting point.
   arma::mat GetInitialPoint() const { return initialPoint; }
 
-  /*
+  /**
    * Evaluate a function for a particular batch-size.
    *
    * @param coordinates The function coordinates.
@@ -71,14 +71,14 @@ class StyblinskiTangFunction
                   const size_t begin,
                   const size_t batchSize) const;
 
-  /*
+  /**
    * Evaluate a function with the given coordinates.
    *
    * @param coordinates The function coordinates.
    */
   double Evaluate(const arma::mat& coordinates) const;
 
-  /*
+  /**
    * Evaluate the gradient of a function for a particular batch-size.
    *
    * @param coordinates The function coordinates.
@@ -91,7 +91,7 @@ class StyblinskiTangFunction
                 arma::mat& gradient,
                 const size_t batchSize) const;
 
-  /*
+  /**
    * Evaluate the gradient of a function with the given coordinates.
    *
    * @param coordinates The function coordinates.
diff --git a/inst/include/ensmallen_bits/problems/three_hump_camel_function.hpp b/inst/include/ensmallen_bits/problems/three_hump_camel_function.hpp
new file mode 100644
index 0000000..73eb1e8
--- /dev/null
+++ b/inst/include/ensmallen_bits/problems/three_hump_camel_function.hpp
@@ -0,0 +1,110 @@
+/**
+ * @file three_hump_camel_function.hpp
+ * @author Suryoday Basak
+ *
+ * Definition of the Three-hump camel function.
+ *
+ * ensmallen 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 ensmallen.  If not, see
+ * http://www.opensource.org/licenses/BSD-3-Clause for more information.
+ */
+#ifndef ENSMALLEN_PROBLEMS_THREE_HUMP_CAMEL_FUNCTION_HPP
+#define ENSMALLEN_PROBLEMS_THREE_HUMP_CAMEL_FUNCTION_HPP
+
+namespace ens {
+namespace test {
+
+/**
+ * The Three-hump camel function, defined by
+ *
+ * \f[
+ * f(x_1,x_2) = 2 * x_1^2 - 1.05 * x_1^4 + (x_1^6)/6 + x_1 * x_2 + x_2^2
+ * \f]
+ *
+ * This should optimize to f(x) = 0, at x = [0, 0].
+ *
+ * For more information, please refer to:
+ *
+ * @code
+ * @article{Ali2005,
+ *   doi       = {10.1007/s10898-004-9972-2},
+ *   year      = {2005},
+ *   month     = apr,
+ *   publisher = {Springer Nature},
+ *   volume    = {31},
+ *   number    = {4},
+ *   pages     = {635--672},
+ *   author    = {M. Montaz Ali and Charoenchai Khompatraporn and
+ *                Zelda B. Zabinsky},
+ *   title     = {A Numerical Evaluation of Several Stochastic Algorithms
+ *                on Selected Continuous Global Optimization Test Problems},
+ *   journal   = {Journal of Global Optimization}
+ * }
+ * @endcode
+ */
+class ThreeHumpCamelFunction
+{
+ public:
+  //! Initialize the ThreeHumpCamelFunction.
+  ThreeHumpCamelFunction();
+
+  /**
+   * Shuffle the order of function visitation. This may be called by the
+   * optimizer.
+   */
+  void Shuffle();
+
+  //! Return 1 (the number of functions).
+  size_t NumFunctions() const { return 1; }
+
+  //! Get the starting point.
+  arma::mat GetInitialPoint() const { return arma::mat("-4.5; 4.5"); }
+
+  /**
+   * Evaluate a function for a particular batch-size.
+   *
+   * @param coordinates The function coordinates.
+   * @param begin The first function.
+   * @param batchSize Number of points to process.
+   */
+  double Evaluate(const arma::mat& coordinates,
+                  const size_t begin,
+                  const size_t batchSize) const;
+
+  /**
+   * Evaluate a function with the given coordinates.
+   *
+   * @param coordinates The function coordinates.
+   */
+  double Evaluate(const arma::mat& coordinates) const;
+
+  /**
+   * Evaluate the gradient of a function for a particular batch-size.
+   *
+   * @param coordinates The function coordinates.
+   * @param begin The first function.
+   * @param gradient The function gradient.
+   * @param batchSize Number of points to process.
+   */
+  void Gradient(const arma::mat& coordinates,
+                const size_t begin,
+                arma::mat& gradient,
+                const size_t batchSize) const;
+
+  /**
+   * Evaluate the gradient of a function with the given coordinates.
+   *
+   * @param coordinates The function coordinates.
+   * @param gradient The function gradient.
+   */
+  void Gradient(const arma::mat& coordinates, arma::mat& gradient);
+};
+
+} // namespace test
+} // namespace ens
+
+// Include implementation.
+#include "three_hump_camel_function_impl.hpp"
+
+#endif // ENSMALLEN_PROBLEMS_THREE_HUMP_CAMEL_FUNCTION_HPP
diff --git a/inst/include/ensmallen_bits/problems/three_hump_camel_function_impl.hpp b/inst/include/ensmallen_bits/problems/three_hump_camel_function_impl.hpp
new file mode 100644
index 0000000..91bb1d4
--- /dev/null
+++ b/inst/include/ensmallen_bits/problems/three_hump_camel_function_impl.hpp
@@ -0,0 +1,70 @@
+/**
+ * @file three_hump_camel_function_impl.hpp
+ * @author Suryoday Basak
+ *
+ * Implementation of the Three-hump camel function.
+ *
+ * ensmallen 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 ensmallen.  If not, see
+ * http://www.opensource.org/licenses/BSD-3-Clause for more information.
+ */
+#ifndef ENSMALLEN_PROBLEMS_THREE_HUMP_CAMEL_FUNCTION_IMPL_HPP
+#define ENSMALLEN_PROBLEMS_THREE_HUMP_CAMEL_FUNCTION_IMPL_HPP
+
+// In case it hasn't been included yet.
+#include "three_hump_camel_function.hpp"
+using namespace std;
+
+namespace ens {
+namespace test {
+
+inline ThreeHumpCamelFunction::ThreeHumpCamelFunction()
+{ /* Nothing to do here */ }
+
+inline void ThreeHumpCamelFunction::Shuffle() { /* Nothing to do here */ }
+
+inline double ThreeHumpCamelFunction::Evaluate(
+    const arma::mat& coordinates,
+    const size_t /* begin */,
+    const size_t /* batchSize */) const
+{
+  // For convenience; we assume these temporaries will be optimized out.
+  const double x1 = coordinates(0);
+  const double x2 = coordinates(1);
+
+  const double objective = (2 * pow(x1, 2)) - (1.05 * pow(x1, 4)) +
+      (pow(x1, 6) / 6) + (x1 * x2) + pow(x2, 2);
+  return objective;
+}
+
+inline double ThreeHumpCamelFunction::Evaluate(
+    const arma::mat& coordinates) const
+{
+  return Evaluate(coordinates, 0, NumFunctions());
+}
+
+inline void ThreeHumpCamelFunction::Gradient(const arma::mat& coordinates,
+                                             const size_t /* begin */,
+                                             arma::mat& gradient,
+                                             const size_t /* batchSize */) const
+{
+  // For convenience; we assume these temporaries will be optimized out.
+  const double x1 = coordinates(0);
+  const double x2 = coordinates(1);
+
+  gradient.set_size(2, 1);
+  gradient(0) = pow(x1, 5) - (4.2 * pow(x1, 3)) + (4 * x1) + x2;
+  gradient(1) = x1 + (2 * x2);
+}
+
+inline void ThreeHumpCamelFunction::Gradient(const arma::mat& coordinates,
+                                             arma::mat& gradient)
+{
+  Gradient(coordinates, 0, gradient, 1);
+}
+
+} // namespace test
+} // namespace ens
+
+#endif
diff --git a/inst/include/ensmallen_bits/problems/wood_function.hpp b/inst/include/ensmallen_bits/problems/wood_function.hpp
index b701061..c621395 100644
--- a/inst/include/ensmallen_bits/problems/wood_function.hpp
+++ b/inst/include/ensmallen_bits/problems/wood_function.hpp
@@ -51,9 +51,9 @@ class WoodFunction
   WoodFunction();
 
   /**
-  * Shuffle the order of function visitation. This may be called by the
-  * optimizer.
-  */
+   * Shuffle the order of function visitation. This may be called by the
+   * optimizer.
+   */
   void Shuffle();
 
   //! Return 1 (the number of functions).
@@ -62,7 +62,7 @@ class WoodFunction
   //! Get the starting point.
   arma::mat GetInitialPoint() const { return arma::mat("-3; -1; -3; -1"); }
 
-  /*
+  /**
    * Evaluate a function for a particular batch-size.
    *
    * @param coordinates The function coordinates.
@@ -73,14 +73,14 @@ class WoodFunction
                   const size_t begin,
                   const size_t batchSize) const;
 
-  /*
+  /**
    * Evaluate a function with the given coordinates.
    *
    * @param coordinates The function coordinates.
    */
   double Evaluate(const arma::mat& coordinates) const;
 
-  /*
+  /**
    * Evaluate the gradient of a function for a particular batch-size.
    *
    * @param coordinates The function coordinates.
@@ -93,7 +93,7 @@ class WoodFunction
                 arma::mat& gradient,
                 const size_t batchSize) const;
 
-  /*
+  /**
    * Evaluate the gradient of a function with the given coordinates.
    *
    * @param coordinates The function coordinates.
diff --git a/inst/include/ensmallen_bits/pso/init_policies/default_init.hpp b/inst/include/ensmallen_bits/pso/init_policies/default_init.hpp
new file mode 100644
index 0000000..48380a4
--- /dev/null
+++ b/inst/include/ensmallen_bits/pso/init_policies/default_init.hpp
@@ -0,0 +1,113 @@
+/**
+ * @file default_init.hpp
+ * @author Chintan Soni
+ * @author Suryoday Basak
+ *
+ * The default initialization policy used by the PSO optimizer.
+ *
+ * ensmallen 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 ENSMALLEN_PSO_INIT_POLICIES_DEFAULT_INIT_HPP
+#define ENSMALLEN_PSO_INIT_POLICIES_DEFAULT_INIT_HPP
+#include <assert.h>
+
+namespace ens {
+
+/**
+ * The default initialization policy used by the PSO optimizer. It initializes
+ * particle positions uniformly in [-1, 1], the velocities in [0, 1] personal
+ * bests of the particles to the initial positions, and all fitness values to
+ * std::numeric_limits<double>::max().
+ */
+class DefaultInit
+{
+ public:
+  /**
+   * Constructor for the DefaultInit policy. The policy initializes particle
+   * posiitons in the range [lowerBound, upperBound]. Defaults to [-1, 1].
+   */
+  DefaultInit()
+  {
+    /* Nothing to do.*/
+  }
+
+  /**
+   * The InitializeParticles method of the init policy. Any class that is used
+   * in place of this default must implement this method which is used by the
+   * optimizer.
+   *
+   * @param iterate Coordinates of the initial point for training.
+   * @param numParticles The number of particles in the swarm.
+   * @param lowerBound Lower bound of the position initialization range.
+   * @param upperBound Upper bound of the position initialization range.
+   * @param particlePositions Current positions of particles.
+   * @param particleVelocities Current velocities of particles.
+   * @param particleFitnesses Current fitness values of particles.
+   * @param particleBestPositions Best positions attained by each particle.
+   * @param particleBestFitnesses Best fitness values attained by each particle.
+   */
+  void Initialize(const arma::mat& iterate,
+                  const size_t numParticles,
+                  arma::vec& lowerBound,
+                  arma::vec& upperBound,
+                  arma::cube& particlePositions,
+                  arma::cube& particleVelocities,
+                  arma::mat& particleFitnesses,
+                  arma::cube& particleBestPositions,
+                  arma::mat& particleBestFitnesses)
+  {
+    // Randomly initialize the particle positions.
+    particlePositions.randu(iterate.n_rows, iterate.n_cols, numParticles);
+
+    // Check if lowerBound is equal to upperBound. If equal, reinitialize.
+    arma::umat lbEquality = (lowerBound == upperBound);
+    if (lbEquality.n_rows == 1 && lbEquality(0, 0) == 1)
+    {
+      lowerBound = -arma::ones<arma::vec>(iterate.n_rows);
+      upperBound = arma::ones<arma::vec>(iterate.n_rows);
+    }
+
+    // Check if lowerBound and upperBound are vectors of a single dimension.
+    else if (lbEquality.n_rows == 1 && lbEquality(0, 0) == 0)
+    {
+      const double lbScalar = lowerBound(0);
+      const double ubScalar = upperBound(0);
+      lowerBound = -lbScalar * arma::ones<arma::vec>(iterate.n_rows);
+      upperBound = ubScalar * arma::ones<arma::vec>(iterate.n_rows);
+    }
+
+    // Check the dimensions of lowerBound and upperBound.
+    assert(lowerBound.n_rows == iterate.n_rows && "The dimensions of "
+        "lowerBound are not the same as the dimensions of iterate.");
+    assert(upperBound.n_rows == iterate.n_rows && "The dimensions of "
+        "upperBound are not the same as the dimensions of iterate.");
+    
+    // Distribute particles in [lowerBound, upperBound].
+    for (size_t i = 0; i < numParticles; i++)
+    {
+      particlePositions.slice(i) = particlePositions.slice(i) % 
+          (upperBound - lowerBound) + lowerBound;
+    }
+
+    // Randomly initialize particle velocities.
+    particleVelocities.randu(iterate.n_rows, iterate.n_cols, numParticles);
+
+    // Initialize current fitness values to infinity.
+    particleFitnesses.set_size(numParticles);
+    particleFitnesses.fill(std::numeric_limits<double>::max());
+
+    // Copy to personal best values for first iteration.
+    particleBestPositions = particlePositions;
+    // Initialize personal best fitness values to infinity.
+    particleBestFitnesses.set_size(numParticles);
+    particleBestFitnesses.fill(std::numeric_limits<double>::max());
+  }
+
+};
+
+} // ens
+
+#endif
diff --git a/inst/include/ensmallen_bits/pso/pso.hpp b/inst/include/ensmallen_bits/pso/pso.hpp
new file mode 100644
index 0000000..681f97d
--- /dev/null
+++ b/inst/include/ensmallen_bits/pso/pso.hpp
@@ -0,0 +1,242 @@
+/**
+ * @file pso.hpp
+ * @author Chintan Soni
+ * @author Suryoday Basak
+ *
+ * Particle swarm optimization.
+ *
+ * ensmallen 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 ENSMALLEN_PSO_PSO_HPP
+#define ENSMALLEN_PSO_PSO_HPP
+
+#include "update_policies/lbest_update.hpp"
+#include "init_policies/default_init.hpp"
+
+namespace ens {
+
+/**
+ * Particle Swarm Optimization (PSO) is an evolutionary approach to optimization
+ * that is inspired by flocks or birds or fishes. The fundamental analogy is
+ * that every creature (particle in a swarm) is at a measurable position of
+ * `goodness' (in the context of PSO, called `fitness') in terms of being able
+ * to find food, and this information can be shared amongst the creatures in the
+ * flock, so that iteratively, the entire flock can get close to the optimal
+ * food source. In a more technical respect, the means by which the fitness
+ * information is shared determines the way in which the global optimum is
+ * approached.
+ *
+ * When this information is shared among particles whose fitness is close to
+ * each other (in a sense, the `nearest' neighbors in the fitness space), the
+ * variant of the approach is called the `local-best' or `lbest' PSO
+ * (consequently, it follows a ring-topology in an information-communication
+ * sense); and when this information is globally shared, the variant is called
+ * the `global-best' or `gbest' PSO (consequently, it follows a star-topology in
+ * an information-communication sense).
+ *
+ * For more information, refer to:
+ *
+ * @inproceedings{Kennedy,
+ *    doi = {10.1109/icnn.1995.488968},
+ *    url = {https://doi.org/10.1109/icnn.1995.488968},
+ *    publisher = {{IEEE}},
+ *    author = {J. Kennedy and R. Eberhart},
+ *    title = {Particle swarm optimization},
+ *    booktitle = {Proceedings of {ICNN}{\textquotesingle}95 -
+ *                 International Conference on Neural Networks}
+ * }
+ *
+ * PSO can optimize arbitrary functions. For more details, see the documentation
+ * on function types included with this distribution or on the ensmallen
+ * website.
+ *
+ * For PSO to work, the function being optimized must implement an
+ * ArbitraryFunctionType template parameter. The respective class must implement
+ * the following function:
+ *
+ *    double Evaluate(const arma::mat& x);
+ */
+template<typename VelocityUpdatePolicy = LBestUpdate,
+         typename InitPolicy = DefaultInit>
+class PSOType
+{
+ public:
+  /**
+   * Construct the particle swarm optimizer with the given function and
+   * parameters. 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 numParticles Number of particles in the swarm.
+   * @param lowerBound Lower bound of the coordinates of the initial population.
+   * @param upperBound Upper bound of the coordinates of the initial population.
+   * @param maxIterations Number of iterations allowed.
+   * @param horizonSize Size of the lookback-horizon for computing improvement.
+   * @param impTolerance Improvement threshold for termination.
+   * @param exploitationFactor Influence of the personal best of the particle.
+   * @param explorationFactor Influence of the neighbours of the particle.
+   */
+  PSOType(const size_t numParticles = 64,
+          const arma::vec& lowerBound = arma::ones<arma::vec>(1),
+          const arma::vec& upperBound = arma::ones<arma::vec>(1),
+          const size_t maxIterations = 3000,
+          const size_t horizonSize = 350,
+          const double impTolerance = 1e-10,
+          const double exploitationFactor = 2.05,
+          const double explorationFactor = 2.05,
+          const VelocityUpdatePolicy& velocityUpdatePolicy =
+              VelocityUpdatePolicy(),
+          const InitPolicy& initPolicy = InitPolicy()) :
+          numParticles(numParticles),
+          lowerBound(lowerBound),
+          upperBound(upperBound),
+          maxIterations(maxIterations),
+          horizonSize(horizonSize),
+          impTolerance(impTolerance),
+          exploitationFactor(exploitationFactor),
+          explorationFactor(explorationFactor),
+          velocityUpdatePolicy(velocityUpdatePolicy),
+          initPolicy(initPolicy) { /* Nothing to do. */ }
+
+  /**
+   * Construct the particle swarm optimizer with the given function and
+   * parameters. 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 numParticles Number of particles in the swarm.
+   * @param lowerBound Lower bound of the coordinates of the initial population.
+   * @param upperBound Upper bound of the coordinates of the initial population.
+   * @param maxIterations Number of iterations allowed.
+   * @param horizonSize Size of the lookback-horizon for computing improvement.
+   * @param impTolerance Improvement threshold for termination.
+   * @param exploitationFactor Influence of the personal best of the particle.
+   * @param explorationFactor Influence of the neighbours of the particle.
+   */
+  PSOType(const size_t numParticles,
+          const double lowerBound,
+          const double upperBound,
+          const size_t maxIterations = 3000,
+          const size_t horizonSize = 350,
+          const double impTolerance = 1e-10,
+          const double exploitationFactor = 2.05,
+          const double explorationFactor = 2.05,
+          const VelocityUpdatePolicy& velocityUpdatePolicy =
+              VelocityUpdatePolicy(),
+          const InitPolicy& initPolicy = InitPolicy()) :
+          numParticles(numParticles),
+          lowerBound(lowerBound * arma::ones<arma::vec>(1)),
+          upperBound(upperBound * arma::ones<arma::vec>(1)),
+          maxIterations(maxIterations),
+          horizonSize(horizonSize),
+          impTolerance(impTolerance),
+          exploitationFactor(exploitationFactor),
+          explorationFactor(explorationFactor),
+          velocityUpdatePolicy(velocityUpdatePolicy),
+          initPolicy(initPolicy) {/* Nothing to do. */ }
+
+  /**
+   * Optimize the input function using PSO. The given variable that holds the
+   * initial point will be modified to store the value of the optimum, or the
+   * point where the PSO method stops, and the final objective value is
+   * returned.
+   *
+   * @param FunctionType Type of the function to be optimized.
+   * @param function Function to be optimized.
+   * @param iterate Initial point (will be modified).
+   * @return Objective value of the final point.
+   */
+  template<typename FunctionType>
+  double Optimize(FunctionType& function, arma::mat& iterate);
+
+  //! Retrieve value of numParticles.
+  size_t NumParticles() const { return numParticles; }
+
+  //! Modify value of numParticles.
+  size_t& NumParticles() { return numParticles; }
+
+  //! Retrieve value of lowerBound.
+  const arma::vec& LowerBound() const { return lowerBound; }
+
+  //! Modify value of lowerBound.
+  arma::vec& LowerBound() { return lowerBound; }
+
+  //! Retrieve value of upperBound.
+  const arma::vec& UpperBound() const { return upperBound; }
+
+  //! Modify value of upperBound.
+  arma::vec& UpperBound() { return upperBound; }
+
+  //! Retrieve value of maxIterations.
+  size_t MaxIterations() const { return maxIterations; }
+
+  //! Modify value of maxIterations.
+  size_t& MaxIterations() { return maxIterations; }
+
+  //! Retrieve value of horizonSize.
+  size_t HorizonSize() const { return horizonSize; }
+
+  //! Modify value of horizonSize.
+  size_t& HorizonSize() { return horizonSize; }
+
+  //! Retrieve value of impTolerance.
+  double ImpTolerance() const { return impTolerance; }
+
+  //! Modify value of impTolerance.
+  double& ImpTolerance() { return impTolerance; }
+
+  //! Retrieve value of exploitationFactor.
+  double ExploitationFactor() const { return exploitationFactor; }
+
+  //! Modify value of exploitationFactor.
+  double& ExploitationFactor() { return exploitationFactor; }
+
+  //! Retrieve value of explorationFactor.
+  double ExplorationFactor() const { return explorationFactor; }
+
+  //! Modify value of explorationFactor.
+  double& ExplorationFactor() { return explorationFactor; }
+
+ private:
+
+  //! Number of particles in the swarm.
+  size_t numParticles;
+  //! Lower bound of the initial swarm.
+  arma::vec lowerBound;
+  //! Upper bound of the initial swarm.
+  arma::vec upperBound;
+  //! Maximum number of iterations for which the optimizer will run.
+  size_t maxIterations;
+  //! The number of iterations looked back at for improvement analysis.
+  size_t horizonSize;
+  //! The tolerance for improvement over the horizon.
+  double impTolerance;
+  //! Exploitation factor for lbest version.
+  double exploitationFactor;
+  //! Exploration factor for lbest version.
+  double explorationFactor;
+  //! Particle positions.
+  arma::cube particlePositions;
+  //! Particle velocities.
+  arma::cube particleVelocities;
+  //! Particle fitness values.
+  arma::vec particleFitnesses;
+  //! Best fitness attained by particle so far.
+  arma::vec particleBestFitnesses;
+  //! Position corresponding to the best fitness of particle.
+  arma::cube particleBestPositions;
+  //! Velocity update policy used.
+  VelocityUpdatePolicy velocityUpdatePolicy;
+  //! Particle initialization policy used.
+  InitPolicy initPolicy;
+};
+
+using LBestPSO = PSOType<LBestUpdate>;
+} // ens
+
+#include "pso_impl.hpp"
+
+#endif
diff --git a/inst/include/ensmallen_bits/pso/pso_impl.hpp b/inst/include/ensmallen_bits/pso/pso_impl.hpp
new file mode 100644
index 0000000..7204f43
--- /dev/null
+++ b/inst/include/ensmallen_bits/pso/pso_impl.hpp
@@ -0,0 +1,177 @@
+/**
+ * @file pso_impl.hpp
+ * @author Chintan Soni
+ * @author Suryoday Basak
+ *
+ * Implementation of the particle swarm optimization algorithm.
+ *
+ * ensmallen 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 ENSMALLEN_PSO_PSO_IMPL_HPP
+#define ENSMALLEN_PSO_PSO_IMPL_HPP
+
+#include "pso.hpp"
+#include <ensmallen_bits/function.hpp>
+#include <queue>
+
+namespace ens {
+/* After the velocity of each particle is updated at the end of each iteration
+ * in PSO, the position of particle i (in iteration j) is updated as:
+ * 
+ * \f[
+ *     p_{i, j + 1} = p_{i, j} + v_{i, j}
+ * \f]
+ *
+ * pso_impl.hpp implements the position-updating procedure. The velocity update
+ * may be done using either the lbest or gbest methods, by using the
+ * appropriate templates.
+ */
+
+//! Optimize the function (minimize).
+template<typename VelocityUpdatePolicy,
+         typename InitPolicy>
+template<typename FunctionType>
+double PSOType<VelocityUpdatePolicy, InitPolicy>::Optimize(
+    FunctionType& function, arma::mat& iterate)
+{
+  /* The following cast is made to make sure that PSO can run on arbitrary
+   * continuous functions and to ensure that we have all the necessary functions.
+   */
+  typedef Function<FunctionType> ArbitraryFunctionType;
+  ArbitraryFunctionType& f(static_cast<ArbitraryFunctionType&>(function));
+  
+  // Make sure we have the methods that we need.
+  traits::CheckEvaluate<ArbitraryFunctionType>();
+
+  // Initialize particles using the init policy.
+  initPolicy.Initialize(iterate,
+      numParticles,
+      lowerBound,
+      upperBound,
+      particlePositions,
+      particleVelocities,
+      particleFitnesses,
+      particleBestPositions,
+      particleBestFitnesses);
+
+  // Initialize the update policy.
+  velocityUpdatePolicy.Initialize(exploitationFactor,
+      explorationFactor,
+      numParticles,
+      iterate);
+
+  // Calculate initial fitness of population.
+  for (size_t i = 0; i < numParticles; i++)
+  {
+    // Calculate fitness value.
+    particleFitnesses(i) = f.Evaluate(particlePositions.slice(i));
+    particleBestFitnesses(i) = particleFitnesses(i);
+  }
+  
+  // Declare queue to keep track of improvements over a number of iterations.
+  queue <double> performanceHorizon;
+  // Variable to store the position of the best particle.
+  size_t bestParticle = 0;
+  // Find the best fitness.
+  double bestFitness = std::numeric_limits<double>::max();
+  
+  // Run PSO for horizonSize number of iterations.
+  // This will allow the performanceHorizon to be updated.
+  // With some initial values in this, we may proceed with the remaining steps
+  // in the PSO method.
+  // The performanceHorizon will be updated with the best particle
+  // in a FIFO manner.
+  for (size_t i = 0; i < horizonSize; i++)
+  {
+    // Calculate fitness and evaluate personal best.
+    for (size_t j = 0; j < numParticles; j++)
+    {
+      particleFitnesses(j) = f.Evaluate(particlePositions.slice(j));
+      // Compare and copy fitness and position to particle best.
+      if (particleFitnesses(j) < particleBestFitnesses(j))
+      {
+        particleBestFitnesses(j) = particleFitnesses(j);
+        particleBestPositions.slice(j) = particlePositions.slice(j);
+      }
+    }
+
+    // Evaluate local best and update velocity.
+    velocityUpdatePolicy.Update(particlePositions,
+        particleVelocities,
+        particleBestPositions,
+        particleBestFitnesses);
+
+    // In-place update of particle positions.
+    particlePositions += particleVelocities;
+
+    // Find the best particle.
+    for (size_t j = 0; j < numParticles; j++)
+    {
+      if (particleBestFitnesses(j) < bestFitness)
+      {
+        bestParticle = j;
+        bestFitness = particleBestFitnesses(bestParticle);
+      }
+    }
+
+    // Append bestFitness to performanceHorizon.
+    performanceHorizon.push(bestFitness);
+  }
+
+  // Run the remaining iterations of PSO.
+  for (size_t i = 0; i < maxIterations - horizonSize; i++)
+  {
+    // Check if there is any improvement over the horizon.
+    // If there is no significant improvement, terminate.
+    if (performanceHorizon.front() - performanceHorizon.back() < impTolerance)
+      break;
+    
+    // Calculate fitness and evaluate personal best.
+    for (size_t j = 0; j < numParticles; j++)
+    {
+      particleFitnesses(j) = f.Evaluate(particlePositions.slice(j));
+      // Compare and copy fitness and position to particle best.
+      if (particleFitnesses(j) < particleBestFitnesses(j))
+      {
+        particleBestFitnesses(j) = particleFitnesses(j);
+        particleBestPositions.slice(j) = particlePositions.slice(j);
+      }
+    }
+
+    // Evaluate local best and update velocity.
+    velocityUpdatePolicy.Update(particlePositions,
+        particleVelocities,
+        particleBestPositions,
+        particleBestFitnesses);
+
+    // In-place update of particle positions.
+    particlePositions += particleVelocities;
+
+    //Find the best particle.
+    for (size_t j = 0; j < numParticles; j++)
+    {
+      if (particleBestFitnesses(j) < bestFitness)
+      {
+        bestParticle = j;
+        bestFitness = particleBestFitnesses(bestParticle);
+      }
+    }
+
+    // Pop the oldest value from performanceHorizon.
+    performanceHorizon.pop();
+    // Push most recent bestFitness to performanceHorizon.
+    performanceHorizon.push(bestFitness);
+  }
+
+  // Copy results back.
+  iterate = particleBestPositions.slice(bestParticle);
+
+  return bestFitness;
+}
+
+} // ens
+
+#endif
diff --git a/inst/include/ensmallen_bits/pso/update_policies/lbest_update.hpp b/inst/include/ensmallen_bits/pso/update_policies/lbest_update.hpp
new file mode 100644
index 0000000..fd0988e
--- /dev/null
+++ b/inst/include/ensmallen_bits/pso/update_policies/lbest_update.hpp
@@ -0,0 +1,160 @@
+/**
+ * @file lbest_update.hpp
+ * @author Chintan Soni
+ * @author Suryoday Basak
+ *
+ * Implementation of the lbest update policy for particle swarm optimization.
+ *
+ * ensmallen 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 ENSMALLEN_PSO_UPDATE_POLICIES_LBEST_UPDATE_HPP
+#define ENSMALLEN_PSO_UPDATE_POLICIES_LBEST_UPDATE_HPP
+#include <assert.h>
+
+namespace ens {
+
+/**
+ * The local best version (lbest) of PSO in which particles communicate with
+ * only two neighbours each, thus forming a ring topology amongst them. This
+ * approach allows PSO to converge at the global minimum, but takes
+ * significantly more iterations to do so.
+ *
+ * The lbest update scheme is described as follows:
+ *
+ * \f{eqation}{
+ * v_{i+1} = \phi (v_i + c_1 * r_1 * (p_{best} - p_{current}) +
+ *           c_1 * r_1 * (l_{best} - p_{current}))
+ * \f}
+ *
+ * where \f$ v_i \f$ is the velocity of a particle in iteration \f$ i \f$,
+ *       \f$ p_{best} \f$ is the best position of an individual particle,
+ *       \f$ p_{current} \f$ is the current position of the particle,
+ *       \f$ l_{best} \f$ is the local best position,
+ *       \f$ r_1 \f$ and \f$ r_2 \f$  are standard uniform random variables,
+ *       \f$ c_1 \f$ is the exploitation factor,
+ *       \f$ c_2 \f$ is the exploration factor, and
+ *       \f$ \phi \f$ is the constriction factor.
+ *
+ * For more information, refer the following:
+ *
+ * @code
+ * @article{Poli2007,
+ *          doi = {10.1007/s11721-007-0002-0},
+ *          url = {https://doi.org/10.1007/s11721-007-0002-0},
+ *          year = {2007},
+ *          month = aug,
+ *          publisher = {Springer},
+ *          volume = {1},
+ *          number = {1},
+ *          pages = {33--57},
+ *          author = {Riccardo Poli and James Kennedy and Tim Blackwell},
+ *          title = {Particle swarm optimization},
+ *          journal = {Swarm Intelligence}
+ * }
+ * @endcod
+ */
+class LBestUpdate
+{
+ public:
+  /**
+   * The Initialize method is called by PSO Optimizer method before the start of
+   * the iteration process. It calculates the value of the constriction
+   * coefficent, initializes the local best indices of each particle to itself,
+   * and sets the shape of the r1 and r2 vectors.
+   *
+   * @param exploitationFactor Influence of personal best achieved.
+   * @param explorationFactor Influence of neighbouring particles.
+   * @param numParticles The number of particles in the swarm.
+   * @param iterate The user input, used for shaping intermediate vectors.
+   */
+  void Initialize(const double exploitationFactor,
+                  const double explorationFactor,
+                  const size_t numParticles,
+                  arma::mat& iterate)
+  {
+    // Copy values to aliases.
+    n = numParticles;
+    c1 = exploitationFactor;
+    c2 = explorationFactor;
+
+    // Calculate the constriction factor
+    static double phi = c1 + c2;
+    assert(phi > 4.0 && "The sum of the exploitation and exploration factors "
+        "must be greater than 4.");
+
+    chi = 2.0 / std::abs(2.0 - phi - std::sqrt((phi - 4.0) * phi));
+
+    // Initialize local best indices to self indices of particles.
+    localBestIndices = arma::linspace<arma::vec>(0, n-1, n);
+
+    // Set sizes r1 and r2.
+    r1.set_size(iterate.n_rows, iterate.n_cols);
+    r2.set_size(iterate.n_rows, iterate.n_cols);
+  }
+
+  /**
+   * Update step for LBestPSO. Compares personal best of each particle with that
+   * of its neighbours, and sets the best of the 3 as the lobal best. This
+   * particle is then used for calculating the velocity for the update step.
+   *
+   * @param particlePositions The current coordinates of particles.
+   * @param particleVelocities The current velocities (will be modified).
+   * @param particleFitnesses The current fitness values or particles.
+   * @param particleBestPositions The personal best coordinates of particles.
+   * @param particleBestFitnesses The personal best fitness values of particles.
+   */
+  void Update(arma::cube& particlePositions,
+              arma::cube& particleVelocities,
+              arma::cube& particleBestPositions,
+              arma::vec& particleBestFitnesses)
+  {
+    // Velocity update logic.
+    for (size_t i = 0; i < n; i++)
+    {
+      localBestIndices(i) =
+          particleBestFitnesses(left(i)) < particleBestFitnesses(i) ?
+          left(i) : i;
+      localBestIndices(i) =
+          particleBestFitnesses(right(i)) < particleBestFitnesses(i) ?
+          right(i) : i;
+    }
+
+    for (size_t i = 0; i < n; i++)
+    {
+      // Generate random numbers for current particle.
+      r1.randu();
+      r2.randu();
+      particleVelocities.slice(i) = chi * (particleVelocities.slice(i) +
+          c1 * r1 % (particleBestPositions.slice(i) -
+          particlePositions.slice(i)) + c2 * r2 %
+          (particleBestPositions.slice(localBestIndices(i)) -
+          particlePositions.slice(i)));
+    }
+  }
+
+ private:
+  //! Number of particles.
+  size_t n;
+  //! Exploitation factor.
+  double c1;
+  //! Exploration factor.
+  double c2;
+  //! Constriction factor chi.
+  double chi;
+  //! Vectors of random numbers.
+  arma::mat r1;
+  arma::mat r2;
+  //! Indices of each particle's best neighbour.
+  arma::vec localBestIndices;
+
+  // Helper functions for calculating neighbours.
+  inline size_t left(size_t index) { return (index + n - 1) % n; }
+  inline size_t right(size_t index) { return (index + 1) % n; }
+};
+
+} // ens
+
+#endif
diff --git a/inst/include/ensmallen_bits/qhadam/qhadam.hpp b/inst/include/ensmallen_bits/qhadam/qhadam.hpp
index fc10aac..16510af 100644
--- a/inst/include/ensmallen_bits/qhadam/qhadam.hpp
+++ b/inst/include/ensmallen_bits/qhadam/qhadam.hpp
@@ -66,6 +66,8 @@ class QHAdam
    *     function is visited in linear order.
    * @param resetPolicy If true, parameters are reset before every Optimize
    *     call; otherwise, their values are retained.
+   * @param exactObjective Calculate the exact objective (Default: estimate the
+   *        final objective obtained on the last pass over the data).
    */
   QHAdam(const double stepSize = 0.001,
          const size_t batchSize = 32,
@@ -77,7 +79,8 @@ class QHAdam
          const size_t maxIterations = 100000,
          const double tolerance = 1e-5,
          const bool shuffle = true,
-         const bool resetPolicy = true);
+         const bool resetPolicy = true,
+         const bool exactObjective = false);
 
   /**
    * Optimize the given function using QHAdam. The given starting point will be
@@ -135,6 +138,11 @@ class QHAdam
   //! Modify whether or not the individual functions are shuffled.
   bool& Shuffle() { return optimizer.Shuffle(); }
 
+  //! Get whether or not the actual objective is calculated.
+  bool ExactObjective() const { return optimizer.ExactObjective(); }
+  //! Modify whether or not the actual objective is calculated.
+  bool& ExactObjective() { return optimizer.ExactObjective(); }
+
   //! Get whether or not the update policy parameters are reset before
   //! Optimize call.
   bool ResetPolicy() const { return optimizer.ResetPolicy(); }
diff --git a/inst/include/ensmallen_bits/qhadam/qhadam_impl.hpp b/inst/include/ensmallen_bits/qhadam/qhadam_impl.hpp
index a559aee..f511e16 100644
--- a/inst/include/ensmallen_bits/qhadam/qhadam_impl.hpp
+++ b/inst/include/ensmallen_bits/qhadam/qhadam_impl.hpp
@@ -28,7 +28,8 @@ inline QHAdam::QHAdam(
     const size_t maxIterations,
     const double tolerance,
     const bool shuffle,
-    const bool resetPolicy) :
+    const bool resetPolicy,
+    const bool exactObjective) :
     optimizer(stepSize,
               batchSize,
               maxIterations,
@@ -36,7 +37,8 @@ inline QHAdam::QHAdam(
               shuffle,
               QHAdamUpdate(epsilon, beta1, beta2, v1, v2),
               NoDecay(),
-              resetPolicy)
+              resetPolicy,
+              exactObjective)
 { /* Nothing to do. */ }
 
 } // namespace ens
diff --git a/inst/include/ensmallen_bits/qhadam/qhadam_update.hpp b/inst/include/ensmallen_bits/qhadam/qhadam_update.hpp
index 1189214..1fb4bd5 100644
--- a/inst/include/ensmallen_bits/qhadam/qhadam_update.hpp
+++ b/inst/include/ensmallen_bits/qhadam/qhadam_update.hpp
@@ -34,7 +34,7 @@ namespace ens {
  */
 class QHAdamUpdate
 {
-public:
+ public:
   /**
    * Construct the QHAdam update policy with the given parameters.
    *
@@ -50,12 +50,12 @@ class QHAdamUpdate
                const double beta2 = 0.999,
                const double v1 = 0.7,
                const double v2 = 1) :
-  epsilon(epsilon),
-  beta1(beta1),
-  beta2(beta2),
-  v1(v1),
-  v2(v2),
-  iteration(0)
+    epsilon(epsilon),
+    beta1(beta1),
+    beta2(beta2),
+    v1(v1),
+    v2(v2),
+    iteration(0)
   {
     // Nothing to do.
   }
@@ -102,8 +102,8 @@ class QHAdamUpdate
 
     // QHAdam recovers Adam when v2 = v1 = 1.
     iterate -= stepSize * ((((1 - v1) * gradient) + v1 * mDash) /
-      (arma::sqrt(((1 - v2) * (gradient % gradient)) +
-        v2 * vDash) + epsilon));
+               (arma::sqrt(((1 - v2) * (gradient % gradient)) +
+               v2 * vDash) + epsilon));
   }
 
   //! Get the value used to initialise the squared gradient parameter.
@@ -131,7 +131,7 @@ class QHAdamUpdate
   //! Modify the second quasi-hyperbolic term.
   double& V2() { return v2; }
 
-private:
+ private:
   // The epsilon value used to initialise the squared gradient parameter.
   double epsilon;
 
@@ -152,7 +152,7 @@ class QHAdamUpdate
 
   // The second quasi-hyperbolic term.
   double v2;
-
+  
   // The number of iterations.
   double iteration;
 };
diff --git a/inst/include/ensmallen_bits/rmsprop/rmsprop.hpp b/inst/include/ensmallen_bits/rmsprop/rmsprop.hpp
index 965c3a7..9de0752 100644
--- a/inst/include/ensmallen_bits/rmsprop/rmsprop.hpp
+++ b/inst/include/ensmallen_bits/rmsprop/rmsprop.hpp
@@ -68,6 +68,8 @@ class RMSProp
    *        function is visited in linear order.
    * @param resetPolicy If true, parameters are reset before every Optimize
    *        call; otherwise, their values are retained.
+   * @param exactObjective Calculate the exact objective (Default: estimate the
+   *        final objective obtained on the last pass over the data).
    */
   RMSProp(const double stepSize = 0.01,
           const size_t batchSize = 32,
@@ -76,7 +78,8 @@ class RMSProp
           const size_t maxIterations = 100000,
           const double tolerance = 1e-5,
           const bool shuffle = true,
-          const bool resetPolicy = true) :
+          const bool resetPolicy = true,
+          const bool exactObjective = false) :
       optimizer(stepSize,
                 batchSize,
                 maxIterations,
@@ -84,7 +87,8 @@ class RMSProp
                 shuffle,
                 RMSPropUpdate(epsilon, alpha),
                 NoDecay(),
-                resetPolicy)
+                resetPolicy,
+                exactObjective)
   { /* Nothing to do. */ }
 
   /**
@@ -138,6 +142,11 @@ class RMSProp
   //! Modify whether or not the individual functions are shuffled.
   bool& Shuffle() { return optimizer.Shuffle(); }
 
+  //! Get whether or not the actual objective is calculated.
+  bool ExactObjective() const { return optimizer.ExactObjective(); }
+  //! Modify whether or not the actual objective is calculated.
+  bool& ExactObjective() { return optimizer.ExactObjective(); }
+
   //! Get whether or not the update policy parameters
   //! are reset before Optimize call.
   bool ResetPolicy() const { return optimizer.ResetPolicy(); }
diff --git a/inst/include/ensmallen_bits/sarah/sarah.hpp b/inst/include/ensmallen_bits/sarah/sarah.hpp
index 91194d3..17b9c20 100644
--- a/inst/include/ensmallen_bits/sarah/sarah.hpp
+++ b/inst/include/ensmallen_bits/sarah/sarah.hpp
@@ -67,6 +67,8 @@ class SARAHType
    *     function is visited in linear order.
    * @param updatePolicy Instantiated update policy used to adjust the given
    *     parameters.
+   * @param exactObjective Calculate the exact objective (Default: estimate the
+   *        final objective obtained on the last pass over the data).
    */
   SARAHType(const double stepSize = 0.01,
             const size_t batchSize = 32,
@@ -74,7 +76,8 @@ class SARAHType
             const size_t innerIterations = 0,
             const double tolerance = 1e-5,
             const bool shuffle = true,
-            const UpdatePolicyType& updatePolicy = UpdatePolicyType());
+            const UpdatePolicyType& updatePolicy = UpdatePolicyType(),
+            const bool exactObjective = false);
 
   /**
    * Optimize the given function using SARAH. The given starting point will be
@@ -119,6 +122,11 @@ class SARAHType
   //! Modify whether or not the individual functions are shuffled.
   bool& Shuffle() { return shuffle; }
 
+  //! Get whether or not the actual objective is calculated.
+  bool ExactObjective() const { return exactObjective; }
+  //! Modify whether or not the actual objective is calculated.
+  bool& ExactObjective() { return exactObjective; }
+
   //! Get the update policy.
   const UpdatePolicyType& UpdatePolicy() const { return updatePolicy; }
   //! Modify the update policy.
@@ -144,6 +152,9 @@ class SARAHType
   //! iterating.
   bool shuffle;
 
+  //! Controls whether or not the actual Objective value is calculated.
+  bool exactObjective;
+
   //! The update policy used to update the parameters in each iteration.
   UpdatePolicyType updatePolicy;
 };
diff --git a/inst/include/ensmallen_bits/sarah/sarah_impl.hpp b/inst/include/ensmallen_bits/sarah/sarah_impl.hpp
index e6e6414..c8c60b6 100644
--- a/inst/include/ensmallen_bits/sarah/sarah_impl.hpp
+++ b/inst/include/ensmallen_bits/sarah/sarah_impl.hpp
@@ -27,13 +27,15 @@ SARAHType<UpdatePolicyType>::SARAHType(
     const size_t innerIterations,
     const double tolerance,
     const bool shuffle,
-    const UpdatePolicyType& updatePolicy) :
+    const UpdatePolicyType& updatePolicy,
+    const bool exactObjective) :
     stepSize(stepSize),
     batchSize(batchSize),
     maxIterations(maxIterations),
     innerIterations(innerIterations),
     tolerance(tolerance),
     shuffle(shuffle),
+    exactObjective(exactObjective),
     updatePolicy(updatePolicy)
 { /* Nothing to do. */ }
 
@@ -177,12 +179,15 @@ double SARAHType<UpdatePolicyType>::Optimize(
       << "terminating optimization." << std::endl;
 
   // Calculate final objective.
-  overallObjective = 0;
-  for (size_t i = 0; i < numFunctions; i += batchSize)
+  if (exactObjective)
   {
-    const size_t effectiveBatchSize = std::min(batchSize, numFunctions - i);
-    overallObjective += function.Evaluate(iterate, i, effectiveBatchSize);
-  }
+    overallObjective = 0;
+    for (size_t i = 0; i < numFunctions; i += batchSize)
+    {
+      const size_t effectiveBatchSize = std::min(batchSize, numFunctions - i);
+      overallObjective += function.Evaluate(iterate, i, effectiveBatchSize);
+    }
+  }  
   return overallObjective;
 }
 
diff --git a/inst/include/ensmallen_bits/sgd/sgd.hpp b/inst/include/ensmallen_bits/sgd/sgd.hpp
index 97a7f85..720b66c 100644
--- a/inst/include/ensmallen_bits/sgd/sgd.hpp
+++ b/inst/include/ensmallen_bits/sgd/sgd.hpp
@@ -88,6 +88,8 @@ class SGD
    * @param decayPolicy Instantiated decay policy used to adjust the step size.
    * @param resetPolicy Flag that determines whether update policy parameters
    *                    are reset before every Optimize call.
+   * @param exactObjective Calculate the exact objective (Default: estimate the
+   *        final objective obtained on the last pass over the data).
    */
   SGD(const double stepSize = 0.01,
       const size_t batchSize = 32,
@@ -96,7 +98,8 @@ class SGD
       const bool shuffle = true,
       const UpdatePolicyType& updatePolicy = UpdatePolicyType(),
       const DecayPolicyType& decayPolicy = DecayPolicyType(),
-      const bool resetPolicy = true);
+      const bool resetPolicy = true,
+      const bool exactObjective = false);
 
   /**
    * Optimize the given function using stochastic gradient descent.  The given
@@ -137,6 +140,11 @@ class SGD
   //! Modify whether or not the individual functions are shuffled.
   bool& Shuffle() { return shuffle; }
 
+  //! Get whether or not the actual objective is calculated.
+  bool ExactObjective() const { return exactObjective; }
+  //! Modify whether or not the actual objective is calculated.
+  bool& ExactObjective() { return exactObjective; }
+
   //! Get whether or not the update policy parameters
   //! are reset before Optimize call.
   bool ResetPolicy() const { return resetPolicy; }
@@ -171,6 +179,9 @@ class SGD
   //! iterating.
   bool shuffle;
 
+  //! Controls whether or not the actual Objective value is calculated.
+  bool exactObjective;
+
   //! The update policy used to update the parameters in each iteration.
   UpdatePolicyType updatePolicy;
 
diff --git a/inst/include/ensmallen_bits/sgd/sgd_impl.hpp b/inst/include/ensmallen_bits/sgd/sgd_impl.hpp
index 6d5ed41..ed8c602 100644
--- a/inst/include/ensmallen_bits/sgd/sgd_impl.hpp
+++ b/inst/include/ensmallen_bits/sgd/sgd_impl.hpp
@@ -30,12 +30,14 @@ SGD<UpdatePolicyType, DecayPolicyType>::SGD(
     const bool shuffle,
     const UpdatePolicyType& updatePolicy,
     const DecayPolicyType& decayPolicy,
-    const bool resetPolicy) :
+    const bool resetPolicy,
+    const bool exactObjective) :
     stepSize(stepSize),
     batchSize(batchSize),
     maxIterations(maxIterations),
     tolerance(tolerance),
     shuffle(shuffle),
+    exactObjective(exactObjective),
     updatePolicy(updatePolicy),
     decayPolicy(decayPolicy),
     resetPolicy(resetPolicy),
@@ -134,12 +136,15 @@ double SGD<UpdatePolicyType, DecayPolicyType>::Optimize(
   Info << "SGD: maximum iterations (" << maxIterations << ") reached; "
       << "terminating optimization." << std::endl;
 
-  // Calculate final objective.
-  overallObjective = 0;
-  for (size_t i = 0; i < numFunctions; i += batchSize)
+  // Calculate final objective if exactObjective is set to true.
+  if (exactObjective)
   {
-    const size_t effectiveBatchSize = std::min(batchSize, numFunctions - i);
-    overallObjective += f.Evaluate(iterate, i, effectiveBatchSize);
+    overallObjective = 0;
+    for (size_t i = 0; i < numFunctions; i += batchSize)
+    {
+      const size_t effectiveBatchSize = std::min(batchSize, numFunctions - i);
+      overallObjective += f.Evaluate(iterate, i, effectiveBatchSize);
+    }
   }
   return overallObjective;
 }
diff --git a/inst/include/ensmallen_bits/sgdr/sgdr.hpp b/inst/include/ensmallen_bits/sgdr/sgdr.hpp
index f181783..3b1a1ec 100644
--- a/inst/include/ensmallen_bits/sgdr/sgdr.hpp
+++ b/inst/include/ensmallen_bits/sgdr/sgdr.hpp
@@ -71,6 +71,8 @@ class SGDR
    *        parameters.
    * @param resetPolicy If true, parameters are reset before every Optimize
    *        call; otherwise, their values are retained.
+   * @param exactObjective Calculate the exact objective (Default: estimate the
+   *        final objective obtained on the last pass over the data).
    */
   SGDR(const size_t epochRestart = 50,
        const double multFactor = 2.0,
@@ -80,7 +82,8 @@ class SGDR
        const double tolerance = 1e-5,
        const bool shuffle = true,
        const UpdatePolicyType& updatePolicy = UpdatePolicyType(),
-       const bool resetPolicy = true);
+       const bool resetPolicy = true,
+       const bool exactObjective = false);
 
   /**
    * Optimize the given function using SGDR.  The given starting point
@@ -120,6 +123,11 @@ class SGDR
   //! Modify whether or not the individual functions are shuffled.
   bool& Shuffle() { return optimizer.Shuffle(); }
 
+  //! Get whether or not the actual objective is calculated.
+  bool ExactObjective() const { return optimizer.ExactObjective(); }
+  //! Modify whether or not the actual objective is calculated.
+  bool& ExactObjective() { return optimizer.ExactObjective(); }
+
   //! Get the update policy.
   const UpdatePolicyType& UpdatePolicy() const
   {
diff --git a/inst/include/ensmallen_bits/sgdr/sgdr_impl.hpp b/inst/include/ensmallen_bits/sgdr/sgdr_impl.hpp
index 0392d00..f8fefe0 100644
--- a/inst/include/ensmallen_bits/sgdr/sgdr_impl.hpp
+++ b/inst/include/ensmallen_bits/sgdr/sgdr_impl.hpp
@@ -27,7 +27,8 @@ SGDR<UpdatePolicyType>::SGDR(
     const double tolerance,
     const bool shuffle,
     const UpdatePolicyType& updatePolicy,
-    const bool resetPolicy) :
+    const bool resetPolicy,
+    const bool exactObjective) :
     batchSize(batchSize),
     optimizer(OptimizerType(stepSize,
                             batchSize,
@@ -39,7 +40,8 @@ SGDR<UpdatePolicyType>::SGDR(
                                 epochRestart,
                                 multFactor,
                                 stepSize),
-                            resetPolicy))
+                            resetPolicy,
+                            exactObjective))
 {
   /* Nothing to do here */
 }
diff --git a/inst/include/ensmallen_bits/sgdr/snapshot_sgdr.hpp b/inst/include/ensmallen_bits/sgdr/snapshot_sgdr.hpp
index e5050dd..37f03c2 100644
--- a/inst/include/ensmallen_bits/sgdr/snapshot_sgdr.hpp
+++ b/inst/include/ensmallen_bits/sgdr/snapshot_sgdr.hpp
@@ -86,6 +86,8 @@ class SnapshotSGDR
    *        parameters.
    * @param resetPolicy If true, parameters are reset before every Optimize
    *        call; otherwise, their values are retained.
+   * @param exactObjective Calculate the exact objective (Default: estimate the
+   *        final objective obtained on the last pass over the data).
    */
   SnapshotSGDR(const size_t epochRestart = 50,
                const double multFactor = 2.0,
@@ -97,7 +99,8 @@ class SnapshotSGDR
                const size_t snapshots = 5,
                const bool accumulate = true,
                const UpdatePolicyType& updatePolicy = UpdatePolicyType(),
-               const bool resetPolicy = true);
+               const bool resetPolicy = true,
+               const bool exactObjective = false);
 
   /**
    * Optimize the given function using SGDR.  The given starting point
@@ -136,6 +139,11 @@ class SnapshotSGDR
   //! Modify whether or not the individual functions are shuffled.
   bool& Shuffle() { return optimizer.Shuffle(); }
 
+  //! Get whether or not the actual objective is calculated.
+  bool ExactObjective() const { return optimizer.ExactObjective(); }
+  //! Modify whether or not the actual objective is calculated.
+  bool& ExactObjective() { return optimizer.ExactObjective(); }
+
   //! Get the snapshots.
   std::vector<arma::mat> Snapshots() const
   {
@@ -177,6 +185,9 @@ class SnapshotSGDR
   //! Whether or not to accumulate the snapshots.
   bool accumulate;
 
+  //! Controls whether or not the actual Objective value is calculated.
+  bool exactObjective;
+
   //! Locally-stored optimizer instance.
   OptimizerType optimizer;
 };
diff --git a/inst/include/ensmallen_bits/sgdr/snapshot_sgdr_impl.hpp b/inst/include/ensmallen_bits/sgdr/snapshot_sgdr_impl.hpp
index fee3035..51128ef 100644
--- a/inst/include/ensmallen_bits/sgdr/snapshot_sgdr_impl.hpp
+++ b/inst/include/ensmallen_bits/sgdr/snapshot_sgdr_impl.hpp
@@ -29,9 +29,11 @@ SnapshotSGDR<UpdatePolicyType>::SnapshotSGDR(
     const size_t snapshots,
     const bool accumulate,
     const UpdatePolicyType& updatePolicy,
-    const bool resetPolicy) :
+    const bool resetPolicy,
+    const bool exactObjective) :
     batchSize(batchSize),
     accumulate(accumulate),
+    exactObjective(exactObjective),
     optimizer(OptimizerType(stepSize,
                             batchSize,
                             maxIterations,
@@ -44,7 +46,8 @@ SnapshotSGDR<UpdatePolicyType>::SnapshotSGDR(
                                 stepSize,
                                 maxIterations,
                                 snapshots),
-                            resetPolicy))
+                            resetPolicy,
+                            exactObjective))
 {
   /* Nothing to do here */
 }
@@ -83,10 +86,13 @@ double SnapshotSGDR<UpdatePolicyType>::Optimize(
     }
     iterate /= (optimizer.DecayPolicy().Snapshots().size() + 1);
 
-    // Calculate final objective.
-    overallObjective = 0;
-    for (size_t i = 0; i < function.NumFunctions(); ++i)
-      overallObjective += function.Evaluate(iterate, i, 1);
+    // Calculate final objective if exactObjective is true.
+    if (exactObjective)
+    {
+      overallObjective = 0;
+      for (size_t i = 0; i < function.NumFunctions(); ++i)
+        overallObjective += function.Evaluate(iterate, i, 1);
+    }
   }
 
   return overallObjective;
diff --git a/inst/include/ensmallen_bits/smorms3/smorms3.hpp b/inst/include/ensmallen_bits/smorms3/smorms3.hpp
index a6da0a1..f723306 100644
--- a/inst/include/ensmallen_bits/smorms3/smorms3.hpp
+++ b/inst/include/ensmallen_bits/smorms3/smorms3.hpp
@@ -62,6 +62,8 @@ class SMORMS3
    *        function is visited in linear order.
    * @param resetPolicy If true, parameters are reset before every Optimize
    *        call; otherwise, their values are retained.
+   * @param exactObjective Calculate the exact objective (Default: estimate the
+   *        final objective obtained on the last pass over the data).
    */
   SMORMS3(const double stepSize = 0.001,
           const size_t batchSize = 32,
@@ -69,7 +71,8 @@ class SMORMS3
           const size_t maxIterations = 100000,
           const double tolerance = 1e-5,
           const bool shuffle = true,
-          const bool resetPolicy = true);
+          const bool resetPolicy = true,
+          const bool exactObjective = false);
 
   /**
    * Optimize the given function using SMORMS3. The given starting point will
@@ -117,6 +120,11 @@ class SMORMS3
   //! Modify whether or not the individual functions are shuffled.
   bool& Shuffle() { return optimizer.Shuffle(); }
 
+  //! Get whether or not the actual objective is calculated.
+  bool ExactObjective() const { return optimizer.ExactObjective(); }
+  //! Modify whether or not the actual objective is calculated.
+  bool& ExactObjective() { return optimizer.ExactObjective(); }
+
   //! Get whether or not the update policy parameters
   //! are reset before Optimize call.
   bool ResetPolicy() const { return optimizer.ResetPolicy(); }
diff --git a/inst/include/ensmallen_bits/smorms3/smorms3_impl.hpp b/inst/include/ensmallen_bits/smorms3/smorms3_impl.hpp
index 66d0e5e..9a47343 100644
--- a/inst/include/ensmallen_bits/smorms3/smorms3_impl.hpp
+++ b/inst/include/ensmallen_bits/smorms3/smorms3_impl.hpp
@@ -23,7 +23,8 @@ inline SMORMS3::SMORMS3(const double stepSize,
                         const size_t maxIterations,
                         const double tolerance,
                         const bool shuffle,
-                        const bool resetPolicy) :
+                        const bool resetPolicy,
+                        const bool exactObjective) :
     optimizer(stepSize,
               batchSize,
               maxIterations,
@@ -31,7 +32,8 @@ inline SMORMS3::SMORMS3(const double stepSize,
               shuffle,
               SMORMS3Update(epsilon),
               NoDecay(),
-              resetPolicy)
+              resetPolicy,
+              exactObjective)
 { /* Nothing to do. */ }
 
 } // namespace ens
diff --git a/inst/include/ensmallen_bits/spalera_sgd/spalera_sgd.hpp b/inst/include/ensmallen_bits/spalera_sgd/spalera_sgd.hpp
index b1fb148..b40bb45 100644
--- a/inst/include/ensmallen_bits/spalera_sgd/spalera_sgd.hpp
+++ b/inst/include/ensmallen_bits/spalera_sgd/spalera_sgd.hpp
@@ -96,6 +96,8 @@ class SPALeRASGD
    * @param decayPolicy Instantiated decay policy used to adjust the step size.
    * @param resetPolicy Flag that determines whether update policy parameters
    *    are reset before every Optimize call.
+   * @param exactObjective Calculate the exact objective (Default: estimate the
+   *        final objective obtained on the last pass over the data).
    */
   SPALeRASGD(const double stepSize = 0.01,
              const size_t batchSize = 32,
@@ -107,7 +109,8 @@ class SPALeRASGD
              const double adaptRate = 3.10e-8,
              const bool shuffle = true,
              const DecayPolicyType& decayPolicy = DecayPolicyType(),
-             const bool resetPolicy = true);
+             const bool resetPolicy = true,
+             const bool exactObjective = false);
 
   /**
    * Optimize the given function using SPALeRA SGD.  The given starting point
@@ -157,6 +160,11 @@ class SPALeRASGD
   //! Modify whether or not the individual functions are shuffled.
   bool& Shuffle() { return shuffle; }
 
+  //! Get whether or not the actual objective is calculated.
+  bool ExactObjective() const { return exactObjective; }
+  //! Modify whether or not the actual objective is calculated.
+  bool& ExactObjective() { return exactObjective; }
+
   //! Get whether or not the update policy parameters
   //! are reset before Optimize call.
   bool ResetPolicy() const { return resetPolicy; }
@@ -194,6 +202,9 @@ class SPALeRASGD
   //! iterating.
   bool shuffle;
 
+  //! Controls whether or not the actual Objective value is calculated.
+  bool exactObjective;
+
   //! The update policy used to update the parameters in each iteration.
   SPALeRAStepsize updatePolicy;
 
diff --git a/inst/include/ensmallen_bits/spalera_sgd/spalera_sgd_impl.hpp b/inst/include/ensmallen_bits/spalera_sgd/spalera_sgd_impl.hpp
index 3109a4f..d2087e9 100644
--- a/inst/include/ensmallen_bits/spalera_sgd/spalera_sgd_impl.hpp
+++ b/inst/include/ensmallen_bits/spalera_sgd/spalera_sgd_impl.hpp
@@ -30,13 +30,15 @@ SPALeRASGD<DecayPolicyType>::SPALeRASGD(const double stepSize,
                                         const double adaptRate,
                                         const bool shuffle,
                                         const DecayPolicyType& decayPolicy,
-                                        const bool resetPolicy) :
+                                        const bool resetPolicy,
+                                        const bool exactObjective) :
     stepSize(stepSize),
     batchSize(batchSize),
     maxIterations(maxIterations),
     tolerance(tolerance),
     lambda(lambda),
     shuffle(shuffle),
+    exactObjective(exactObjective),
     updatePolicy(SPALeRAStepsize(alpha, epsilon, adaptRate)),
     decayPolicy(decayPolicy),
     resetPolicy(resetPolicy)
@@ -150,12 +152,15 @@ double SPALeRASGD<DecayPolicyType>::Optimize(DecomposableFunctionType& function,
   Info << "SPALeRA SGD: maximum iterations (" << maxIterations
       << ") reached; terminating optimization." << std::endl;
 
-  // Calculate final objective.
-  overallObjective = 0;
-  for (size_t i = 0; i < numFunctions; i += batchSize)
+  // Calculate final objective if exactObjective is set to true.
+  if (exactObjective)
   {
-    const size_t effectiveBatchSize = std::min(batchSize, numFunctions - i);
-    overallObjective += f.Evaluate(iterate, i, effectiveBatchSize);
+    overallObjective = 0;
+    for (size_t i = 0; i < numFunctions; i += batchSize)
+    {
+      const size_t effectiveBatchSize = std::min(batchSize, numFunctions - i);
+      overallObjective += f.Evaluate(iterate, i, effectiveBatchSize);
+    }
   }
 
   return overallObjective;
diff --git a/inst/include/ensmallen_bits/svrg/svrg.hpp b/inst/include/ensmallen_bits/svrg/svrg.hpp
index c9383f9..7170915 100644
--- a/inst/include/ensmallen_bits/svrg/svrg.hpp
+++ b/inst/include/ensmallen_bits/svrg/svrg.hpp
@@ -106,6 +106,8 @@ class SVRGType
    * @param decayPolicy Instantiated decay policy used to adjust the step size.
    * @param resetPolicy Flag that determines whether update policy parameters
    *     are reset before every Optimize call.
+   * @param exactObjective Calculate the exact objective (Default: estimate the
+   *        final objective obtained on the last pass over the data).
    */
   SVRGType(const double stepSize = 0.01,
            const size_t batchSize = 32,
@@ -115,7 +117,8 @@ class SVRGType
            const bool shuffle = true,
            const UpdatePolicyType& updatePolicy = UpdatePolicyType(),
            const DecayPolicyType& decayPolicy = DecayPolicyType(),
-           const bool resetPolicy = true);
+           const bool resetPolicy = true,
+           const bool exactObjective = false);
 
   /**
    * Optimize the given function using SVRG. The given starting point will be
@@ -160,6 +163,11 @@ class SVRGType
   //! Modify whether or not the individual functions are shuffled.
   bool& Shuffle() { return shuffle; }
 
+  //! Get whether or not the actual objective is calculated.
+  bool ExactObjective() const { return exactObjective; }
+  //! Modify whether or not the actual objective is calculated.
+  bool& ExactObjective() { return exactObjective; }
+
   //! Get whether or not the update policy parameters
   //! are reset before Optimize call.
   bool ResetPolicy() const { return resetPolicy; }
@@ -197,6 +205,9 @@ class SVRGType
   //! iterating.
   bool shuffle;
 
+  //! Controls whether or not the actual Objective value is calculated.
+  bool exactObjective;
+
   //! The update policy used to update the parameters in each iteration.
   UpdatePolicyType updatePolicy;
 
diff --git a/inst/include/ensmallen_bits/svrg/svrg_impl.hpp b/inst/include/ensmallen_bits/svrg/svrg_impl.hpp
index 75c5ad6..a5c3b73 100644
--- a/inst/include/ensmallen_bits/svrg/svrg_impl.hpp
+++ b/inst/include/ensmallen_bits/svrg/svrg_impl.hpp
@@ -27,13 +27,15 @@ SVRGType<UpdatePolicyType, DecayPolicyType>::SVRGType(
     const bool shuffle,
     const UpdatePolicyType& updatePolicy,
     const DecayPolicyType& decayPolicy,
-    const bool resetPolicy) :
+    const bool resetPolicy,
+    const bool exactObjective) :
     stepSize(stepSize),
     batchSize(batchSize),
     maxIterations(maxIterations),
     innerIterations(innerIterations),
     tolerance(tolerance),
     shuffle(shuffle),
+    exactObjective(exactObjective),
     updatePolicy(updatePolicy),
     decayPolicy(decayPolicy),
     resetPolicy(resetPolicy)
@@ -159,12 +161,15 @@ double SVRGType<UpdatePolicyType, DecayPolicyType>::Optimize(
   Info << "SVRG: maximum iterations (" << maxIterations << ") reached; "
       << "terminating optimization." << std::endl;
 
-  // Calculate final objective.
-  overallObjective = 0;
-  for (size_t i = 0; i < numFunctions; i += batchSize)
+  // Calculate final objective if exactObjective is set to true.
+  if (exactObjective)
   {
-    const size_t effectiveBatchSize = std::min(batchSize, numFunctions - i);
-    overallObjective += function.Evaluate(iterate, i, effectiveBatchSize);
+    overallObjective = 0;
+    for (size_t i = 0; i < numFunctions; i += batchSize)
+    {
+      const size_t effectiveBatchSize = std::min(batchSize, numFunctions - i);
+      overallObjective += function.Evaluate(iterate, i, effectiveBatchSize);
+    }
   }
   return overallObjective;
 }
diff --git a/inst/include/ensmallen_bits/swats/swats.hpp b/inst/include/ensmallen_bits/swats/swats.hpp
index 26f9545..09da031 100644
--- a/inst/include/ensmallen_bits/swats/swats.hpp
+++ b/inst/include/ensmallen_bits/swats/swats.hpp
@@ -66,6 +66,8 @@ class SWATS
    *        function is visited in linear order.
    * @param resetPolicy If true, parameters are reset before every Optimize
    *        call; otherwise, their values are retained.
+   * @param exactObjective Calculate the exact objective (Default: estimate the
+   *        final objective obtained on the last pass over the data).
    */
   SWATS(const double stepSize = 0.001,
         const size_t batchSize = 32,
@@ -75,7 +77,8 @@ class SWATS
         const size_t maxIterations = 100000,
         const double tolerance = 1e-5,
         const bool shuffle = true,
-        const bool resetPolicy = true);
+        const bool resetPolicy = true,
+        const bool exactObjective = false);
 
   /**
    * Optimize the given function using SWATS. The given starting point will
@@ -133,6 +136,11 @@ class SWATS
   //! Modify whether or not the individual functions are shuffled.
   bool& Shuffle() { return optimizer.Shuffle(); }
 
+  //! Get whether or not the actual objective is calculated.
+  bool ExactObjective() const { return optimizer.ExactObjective(); }
+  //! Modify whether or not the actual objective is calculated.
+  bool& ExactObjective() { return optimizer.ExactObjective(); }
+
   //! Get whether or not the update policy parameters
   //! are reset before Optimize call.
   bool ResetPolicy() const { return optimizer.ResetPolicy(); }
diff --git a/inst/include/ensmallen_bits/swats/swats_impl.hpp b/inst/include/ensmallen_bits/swats/swats_impl.hpp
index 56d5645..c69dfd5 100644
--- a/inst/include/ensmallen_bits/swats/swats_impl.hpp
+++ b/inst/include/ensmallen_bits/swats/swats_impl.hpp
@@ -26,7 +26,8 @@ inline SWATS::SWATS(
     const size_t maxIterations,
     const double tolerance,
     const bool shuffle,
-    const bool resetPolicy) :
+    const bool resetPolicy,
+    const bool exactObjective) :
     optimizer(stepSize,
               batchSize,
               maxIterations,
@@ -34,7 +35,8 @@ inline SWATS::SWATS(
               shuffle,
               SWATSUpdate(epsilon, beta1, beta2),
               NoDecay(),
-              resetPolicy)
+              resetPolicy,
+              exactObjective)
 { /* Nothing to do. */ }
 
 } // namespace ens
diff --git a/inst/include/ensmallen_bits/wn_grad/wn_grad.hpp b/inst/include/ensmallen_bits/wn_grad/wn_grad.hpp
index ba661ad..6ce2c49 100644
--- a/inst/include/ensmallen_bits/wn_grad/wn_grad.hpp
+++ b/inst/include/ensmallen_bits/wn_grad/wn_grad.hpp
@@ -58,13 +58,16 @@ class WNGrad
    *        function is visited in linear order.
    * @param resetPolicy If true, parameters are reset before every Optimize
    *        call; otherwise, their values are retained.
+   * @param exactObjective Calculate the exact objective (Default: estimate the
+   *        final objective obtained on the last pass over the data).
    */
   WNGrad(const double stepSize = 0.562,
          const size_t batchSize = 32,
          const size_t maxIterations = 100000,
          const double tolerance = 1e-5,
          const bool shuffle = true,
-         const bool resetPolicy = true);
+         const bool resetPolicy = true,
+         const bool exactObjective = false);
 
   /**
    * Optimize the given function using WNGrad. The given starting point will
@@ -107,6 +110,11 @@ class WNGrad
   //! Modify whether or not the individual functions are shuffled.
   bool& Shuffle() { return optimizer.Shuffle(); }
 
+  //! Get whether or not the actual objective is calculated.
+  bool ExactObjective() const { return optimizer.ExactObjective(); }
+  //! Modify whether or not the actual objective is calculated.
+  bool& ExactObjective() { return optimizer.ExactObjective(); }
+
   //! Get whether or not the update policy parameters
   //! are reset before Optimize call.
   bool ResetPolicy() const { return optimizer.ResetPolicy(); }
diff --git a/inst/include/ensmallen_bits/wn_grad/wn_grad_impl.hpp b/inst/include/ensmallen_bits/wn_grad/wn_grad_impl.hpp
index 20de70c..01bc2d8 100644
--- a/inst/include/ensmallen_bits/wn_grad/wn_grad_impl.hpp
+++ b/inst/include/ensmallen_bits/wn_grad/wn_grad_impl.hpp
@@ -23,7 +23,8 @@ inline WNGrad::WNGrad(
     const size_t maxIterations,
     const double tolerance,
     const bool shuffle,
-    const bool resetPolicy) :
+    const bool resetPolicy,
+    const bool exactObjective) :
     optimizer(stepSize,
               batchSize,
               maxIterations,
@@ -31,7 +32,8 @@ inline WNGrad::WNGrad(
               shuffle,
               WNGradUpdate(),
               NoDecay(),
-              resetPolicy)
+              resetPolicy,
+              exactObjective)
 { /* Nothing to do. */ }
 
 } // namespace ens