diff --git a/include/boost/math/differentiation/detail/reverse_mode_autodiff_memory_management.hpp b/include/boost/math/differentiation/detail/reverse_mode_autodiff_memory_management.hpp
index a626ce5d6..84078c4ac 100644
--- a/include/boost/math/differentiation/detail/reverse_mode_autodiff_memory_management.hpp
+++ b/include/boost/math/differentiation/detail/reverse_mode_autodiff_memory_management.hpp
@@ -282,7 +282,7 @@ public:
   void add_checkpoint()
   {
     if (total_size_ > 0) {
-      checkpoints_.push_back(total_size_ - 1);
+      checkpoints_.push_back(total_size_); // - 1);
     } else {
       checkpoints_.push_back(0);
     }
diff --git a/include/boost/math/optimization/minimizer.hpp b/include/boost/math/optimization/minimizer.hpp
index ea3edfed3..1b3e5b9ae 100644
--- a/include/boost/math/optimization/minimizer.hpp
+++ b/include/boost/math/optimization/minimizer.hpp
@@ -275,8 +275,7 @@ struct unit_sphere_constraint
 {
   void operator()(ArgumentContainer& x) const
   {
-    RealType norm2v = norm_2(x);
-    RealType norm = sqrt(norm2v);
+    RealType norm = norm_2(x);
     if (norm > RealType{ 0 }) {
       for (auto& xi : x)
         xi /= norm;
diff --git a/include/boost/math/optimization/nesterov.hpp b/include/boost/math/optimization/nesterov.hpp
index 961dfa691..efd13a502 100644
--- a/include/boost/math/optimization/nesterov.hpp
+++ b/include/boost/math/optimization/nesterov.hpp
@@ -32,8 +32,9 @@ struct nesterov_update_policy
                ArgumentType>::value>::type>
   void operator()(ArgumentType& x, RealType& g, RealType& v)
   {
+    RealType v_prev = v;
     v = mu_ * v - lr_ * g;
-    x.get_value() += v;
+    x.get_value() += -mu_ * v_prev + (static_cast<RealType>(1) + mu_) * v;
   }
   template<typename ArgumentType,
            typename std::enable_if<!boost::math::differentiation::reverse_mode::
@@ -41,8 +42,9 @@ struct nesterov_update_policy
                                    int>::type = 0>
   void operator()(ArgumentType& x, RealType& g, RealType& v) const
   {
+    const RealType v_prev = v;
     v = mu_ * v - lr_ * g;
-    x += v;
+    x += -mu_ * v_prev + (static_cast<RealType>(1) + mu_) * v;
   }
   RealType lr() const noexcept { return lr_; }
   RealType mu() const noexcept { return mu_; }
diff --git a/test/test_functions_for_optimization.hpp b/test/test_functions_for_optimization.hpp
index a7fbd7cdc..da989fc09 100644
--- a/test/test_functions_for_optimization.hpp
+++ b/test/test_functions_for_optimization.hpp
@@ -14,92 +14,101 @@
 
 /* simple n-d quadratic function */
 template<typename RealType>
-RealType quadratic(std::vector<RealType> &x)
+RealType
+quadratic(std::vector<RealType>& x)
 {
-    RealType res{0.0};
-    for (auto &item : x) {
-        res += item * item;
-    }
-    return res;
+  RealType res{ 0.0 };
+  for (auto& item : x) {
+    res += item * item;
+  }
+  return res;
 }
 
 template<typename RealType>
-RealType quadratic_high_cond_2D(std::vector<RealType> &x)
+RealType
+quadratic_high_cond_2D(std::vector<RealType>& x)
 {
-    return 1000 * x[0] * x[0] + x[1] * x[1];
+  return 1000 * x[0] * x[0] + x[1] * x[1];
 }
 
 // Taken from: https://en.wikipedia.org/wiki/Test_functions_for_optimization
 template<typename Real>
-Real ackley(std::array<Real, 2> const &v)
+Real
+ackley(std::array<Real, 2> const& v)
 {
-    using boost::math::constants::e;
-    using boost::math::constants::two_pi;
-    using std::cos;
-    using std::exp;
-    using std::sqrt;
-    Real x    = v[0];
-    Real y    = v[1];
-    Real arg1 = -sqrt((x * x + y * y) / 2) / 5;
-    Real arg2 = cos(two_pi<Real>() * x) + cos(two_pi<Real>() * y);
-    return -20 * exp(arg1) - exp(arg2 / 2) + 20 + e<Real>();
+  using boost::math::constants::e;
+  using boost::math::constants::two_pi;
+  using std::cos;
+  using std::exp;
+  using std::sqrt;
+  Real x = v[0];
+  Real y = v[1];
+  Real arg1 = -sqrt((x * x + y * y) / 2) / 5;
+  Real arg2 = cos(two_pi<Real>() * x) + cos(two_pi<Real>() * y);
+  return -20 * exp(arg1) - exp(arg2 / 2) + 20 + e<Real>();
 }
 
 template<typename Real>
-auto rosenbrock_saddle(std::array<Real, 2> const &v) -> Real
+auto
+rosenbrock_saddle(std::array<Real, 2> const& v) -> Real
 {
-    Real x{v[0]};
-    Real y{v[1]};
-    return static_cast<Real>(100 * (x * x - y) * (x * x - y) + (1 - x) * (1 - x));
+  Real x{ v[0] };
+  Real y{ v[1] };
+  return static_cast<Real>(100 * (x * x - y) * (x * x - y) + (1 - x) * (1 - x));
 }
 
 template<class Real>
-Real rastrigin(std::vector<Real> const &v)
+Real
+rastrigin(std::vector<Real> const& v)
 {
-    using boost::math::constants::two_pi;
-    using std::cos;
-    auto A = static_cast<Real>(10);
-    auto y = static_cast<Real>(10 * v.size());
-    for (auto x : v) {
-        y += x * x - A * cos(two_pi<Real>() * x);
-    }
-    return y;
+  using boost::math::constants::two_pi;
+  using std::cos;
+  auto A = static_cast<Real>(10);
+  auto y = static_cast<Real>(10 * v.size());
+  for (auto x : v) {
+    y += x * x - A * cos(two_pi<Real>() * x);
+  }
+  return y;
 }
 
 // Useful for testing return-type != scalar argument type,
 // and robustness to NaNs:
-double sphere(std::vector<float> const &v)
+double
+sphere(std::vector<float> const& v)
 {
-    double r = 0.0;
-    for (auto x : v) {
-        double x_  = static_cast<double>(x);
-        r         += x_ * x_;
-    }
-    if (r >= 1) {
-        return std::numeric_limits<double>::quiet_NaN();
-    }
-    return r;
+  double r = 0.0;
+  for (auto x : v) {
+    double x_ = static_cast<double>(x);
+    r += x_ * x_;
+  }
+  if (r >= 1) {
+    return std::numeric_limits<double>::quiet_NaN();
+  }
+  return r;
 }
 
 template<typename Real>
-Real three_hump_camel(std::array<Real, 2> const &v)
+Real
+three_hump_camel(std::array<Real, 2> const& v)
 {
-    Real x   = v[0];
-    Real y   = v[1];
-    auto xsq = x * x;
-    return 2 * xsq - (1 + Real(1) / Real(20)) * xsq * xsq + xsq * xsq * xsq / 6 + x * y + y * y;
+  Real x = v[0];
+  Real y = v[1];
+  auto xsq = x * x;
+  return 2 * xsq - (1 + Real(1) / Real(20)) * xsq * xsq + xsq * xsq * xsq / 6 +
+         x * y + y * y;
 }
 
 // Minima occurs at (3, 1/2) with value 0:
 template<typename Real>
-Real beale(std::array<Real, 2> const &v)
+Real
+beale(std::array<Real, 2> const& v)
 {
-    Real x  = v[0];
-    Real y  = v[1];
-    Real t1 = Real(3) / Real(2) - x + x * y;
-    Real t2 = Real(9) / Real(4) - x + x * y * y;
-    Real t3 = Real(21) / Real(8) - x + x * y * y * y;
-    return t1 * t1 + t2 * t2 + t3 * t3;
+  Real x = v[0];
+  Real y = v[1];
+  Real t1 = Real(3) / Real(2) - x + x * y;
+  Real t2 = Real(9) / Real(4) - x + x * y * y;
+  Real t3 = Real(21) / Real(8) - x + x * y * y * y;
+  return t1 * t1 + t2 * t2 + t3 * t3;
 }
 
 #endif
diff --git a/test/test_nesterov_optimizer.cpp b/test/test_nesterov_optimizer.cpp
index 8bf9642f5..83e43db25 100644
--- a/test/test_nesterov_optimizer.cpp
+++ b/test/test_nesterov_optimizer.cpp
@@ -9,12 +9,11 @@
 #include <boost/math/optimization/nesterov.hpp>
 namespace rdiff = boost::math::differentiation::reverse_mode;
 namespace bopt = boost::math::optimization;
-BOOST_AUTO_TEST_SUITE(basic_gradient_descent)
+BOOST_AUTO_TEST_SUITE(nesterov_descent)
 
 BOOST_AUTO_TEST_CASE_TEMPLATE(default_nesterov_test, T, all_float_types)
 {
-  size_t NITER = 5;
-  T lr = T{ 1e-3 };
+  T lr = T{ 1e-4 };
   T mu = T{ 0.95 };
   RandomSample<T> rng{ T(-10), (10) };
   std::vector<rdiff::rvar<T, 1>> x;