Optimizer: size the dense BFGS inverse Hessian before each update

src/stan/optimization/bfgs_update.hpp

@@ -14,6 +14,15 @@ class BFGSUpdate_HInv {
   /**
    * Update the inverse Hessian approximation.
    *
+   * Except at a reset, the approximation H is multiplied by
+   * tau = s'y / y'Hy before the update if tau > 1 (Oren-Luenberger sizing
+   * restricted to enlarging H). Without it H keeps the scale set at the
+   * last reset, which far from the optimum can underestimate the inverse
+   * curvature by orders of magnitude in the directions not explored since.
+   * Scaling all of H can overestimate it along earlier steps instead; a
+   * later pair along such a direction has y'Hy > s'y, so it is not sized
+   * and the update corrects H there.
+   *
    * @param yk Difference between the current and previous gradient vector.
    * @param sk Difference between the current and previous state vector.
    * @param reset Whether to reset the approximation, forgetting about
@@ -25,21 +34,32 @@ class BFGSUpdate_HInv {
   inline Scalar update(const VectorT &yk, const VectorT &sk,
                        bool reset = false) {
     Scalar rhok, skyk, B0fact;
-    HessianT Hupd;
 
     skyk = yk.dot(sk);
     rhok = 1.0 / skyk;
 
-    Hupd.noalias()
-        = HessianT::Identity(yk.size(), yk.size()) - rhok * sk * yk.transpose();
     if (reset) {
       B0fact = yk.squaredNorm() / skyk;
-      _Hk.noalias() = ((1.0 / B0fact) * Hupd) * Hupd.transpose();
+      _Hk.setIdentity(yk.size(), yk.size());
+      _Hk /= B0fact;
     } else {
       B0fact = 1.0;
-      _Hk = Hupd * _Hk * Hupd.transpose();
     }
-    _Hk.noalias() += rhok * sk * sk.transpose();
+    // H <- (I - rho s y') H (I - rho y s') + rho s s'
+    //    = H + (c s - rho u) s' - rho s u',  u = H y,  c = rho^2 y'u + rho
+    _u.noalias() = _Hk * yk;
+    Scalar yu = yk.dot(_u);
+    if (!reset && 0 < yu && yu < skyk) {
+      const Scalar tau = skyk / yu;
+      _Hk *= tau;
+      _u *= tau;
+      yu = skyk;
+    }
+    const Scalar c = rhok * rhok * yu + rhok;
+    _w = c * sk - rhok * _u;
+    _Hk.noalias() += _w * sk.transpose();
+    _u *= rhok;
+    _Hk.noalias() -= sk * _u.transpose();
     return B0fact;
   }
 
@@ -57,6 +77,7 @@ class BFGSUpdate_HInv {
 
  private:
   HessianT _Hk;
+  VectorT _u, _w;
 };
 }  // namespace optimization
 }  // namespace stan

src/test/unit/optimization/bfgs_minimizer_test.cpp

@@ -339,6 +339,39 @@ TEST_F(OptimizationBfgsMinimizer, step) {
   bfgs.step();
 }
 
+namespace {
+// f(x) = sum_i exp(i x_i) - i x_i, i = 1, ..., n, minimised at 0 with f = n
+struct minimizer_exp_sum {
+  int operator()(const Eigen::VectorXd &x, double &f, Eigen::VectorXd &g) {
+    Eigen::ArrayXd c = Eigen::ArrayXd::LinSpaced(x.size(), 1, x.size());
+    Eigen::ArrayXd e = (c * x.array()).exp();
+    f = (e - c * x.array()).sum();
+    g = (c * (e - 1)).matrix();
+    return std::isfinite(f) && g.allFinite() ? 0 : 1;
+  }
+};
+}  // namespace
+
+TEST(OptimizationBfgsMinimizerExpSum, dense_bfgs_large_initial_curvature) {
+  // From x0 = (5, ..., 5) the curvature falls from 25 e^25 = 1.8e12 to at
+  // most 25 at the optimum. The dense inverse Hessian estimate used to keep
+  // the scale of the first steps in the directions not explored since, and
+  // with CmdStan's default tolerances the relative gradient test stopped the
+  // run at f - 5 = 7.6e3.
+  stan::optimization::BFGSMinimizer<minimizer_exp_sum,
+                                    stan::optimization::BFGSUpdate_HInv<> >
+      bfgs(minimizer_exp_sum{});
+  bfgs._conv_opts.tolRelGrad = 1e7;
+  Eigen::VectorXd x0 = Eigen::VectorXd::Constant(5, 5);
+  bfgs.initialize(x0);
+  int ret = stan::optimization::TERM_SUCCESS;
+  while (ret == stan::optimization::TERM_SUCCESS) {
+    ret = bfgs.step();
+  }
+  EXPECT_GT(ret, 0);
+  EXPECT_NEAR(5, bfgs.curr_f(), 1e-6);
+}
+
 TEST_F(OptimizationBfgsMinimizer, minimize) {
   stan::io::empty_var_context dummy_context;
   Model rb_model(dummy_context);
@@ -347,5 +380,5 @@ TEST_F(OptimizationBfgsMinimizer, minimize) {
   EXPECT_EQ("", out.str());
   Optimizer bfgs(_adaptor);
 
-  EXPECT_FLOAT_EQ(bfgs.minimize(cont_vector), 31);
+  EXPECT_FLOAT_EQ(bfgs.minimize(cont_vector), 20);
 }

src/test/unit/optimization/bfgs_update_test.cpp

@@ -35,6 +35,25 @@ TEST(OptimizationBfgsUpdate, bfgs_update_secant) {
   }
 }
 
+TEST(OptimizationBfgsUpdate, BFGSUpdate_HInv_enlarges_small_estimate) {
+  // The reset on a pair with curvature 1e4 along e1 sets H = 1e-4 I. The next
+  // pair has curvature 1 along e2, so H is first multiplied by 1e4, also
+  // along e3, which no pair has explored (H e3 stayed 1e-4 e3 before). A
+  // pair with curvature 100 along e3, where H is then 1, does not shrink H:
+  // H e1 stays e1, although the first pair measured curvature 1e4 along e1.
+  stan::optimization::BFGSUpdate_HInv<> bfgsUp;
+  Eigen::VectorXd s = Eigen::VectorXd::Unit(3, 0), p;
+  bfgsUp.update(1e4 * s, s, true);
+  s = Eigen::VectorXd::Unit(3, 1);
+  bfgsUp.update(s, s);
+  bfgsUp.search_direction(p, Eigen::VectorXd::Unit(3, 2));
+  EXPECT_FLOAT_EQ(-1, p[2]);
+  s = Eigen::VectorXd::Unit(3, 2);
+  bfgsUp.update(100 * s, s);
+  bfgsUp.search_direction(p, Eigen::VectorXd::Unit(3, 0));
+  EXPECT_FLOAT_EQ(-1, p[0]);
+}
+
 TEST(OptimizationBfgsUpdate, BFGSUpdate_HInv_update) {
   const int nDim = 10;
 

src/test/unit/services/optimize/bfgs_test.cpp

@@ -46,15 +46,15 @@ TEST_F(ServicesOptimize, rosenbrock) {
   EXPECT_EQ("x", parameter.names_[2]);
   EXPECT_EQ("y", parameter.names_[3]);
 
-  EXPECT_EQ(20, parameter.states_.size());
+  EXPECT_EQ(22, parameter.states_.size());
   EXPECT_FLOAT_EQ(0, parameter.states_.front()[2])
       << "initial value should be (0, 0)";
   EXPECT_FLOAT_EQ(0, parameter.states_.front()[3])
       << "initial value should be (0, 0)";
-  EXPECT_FLOAT_EQ(1, parameter.states_.back()[2])
+  EXPECT_FLOAT_EQ(1.0000112, parameter.states_.back()[2])
       << "optimal value should be (1, 1)";
-  EXPECT_FLOAT_EQ(1, parameter.states_.back()[3])
+  EXPECT_FLOAT_EQ(1.0000218, parameter.states_.back()[3])
       << "optimal value should be (1, 1)";
   EXPECT_FLOAT_EQ(return_code, 0);
-  EXPECT_EQ(19, interrupt.call_count());
+  EXPECT_EQ(21, interrupt.call_count());
 }