internal/ceres/low_rank_inverse_hessian.cc

   1 // Ceres Solver - A fast non-linear least squares minimizer
   2 // Copyright 2015 Google Inc. All rights reserved.
   3 // http://ceres-solver.org/
   4 //
   5 // Redistribution and use in source and binary forms, with or without
   6 // modification, are permitted provided that the following conditions are met:
   7 //
   8 // * Redistributions of source code must retain the above copyright notice,
   9 //   this list of conditions and the following disclaimer.
  10 // * Redistributions in binary form must reproduce the above copyright notice,
  11 //   this list of conditions and the following disclaimer in the documentation
  12 //   and/or other materials provided with the distribution.
  13 // * Neither the name of Google Inc. nor the names of its contributors may be
  14 //   used to endorse or promote products derived from this software without
  15 //   specific prior written permission.
  16 //
  17 // THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
  18 // AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
  19 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
  20 // ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
  21 // LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
  22 // CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
  23 // SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
  24 // INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
  25 // CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
  26 // ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
  27 // POSSIBILITY OF SUCH DAMAGE.
  28 //
  29 // Author: sameeragarwal@google.com (Sameer Agarwal)
  30
  31 #include <list>
  32
  33 #include "ceres/internal/eigen.h"
  34 #include "ceres/low_rank_inverse_hessian.h"
  35 #include "glog/logging.h"
  36
  37 namespace ceres {
  38 namespace internal {
  39
  40 using std::list;
  41
  42 // The (L)BFGS algorithm explicitly requires that the secant equation:
  43 //
  44 //   B_{k+1} * s_k = y_k
  45 //
  46 // Is satisfied at each iteration, where B_{k+1} is the approximated
  47 // Hessian at the k+1-th iteration, s_k = (x_{k+1} - x_{k}) and
  48 // y_k = (grad_{k+1} - grad_{k}). As the approximated Hessian must be
  49 // positive definite, this is equivalent to the condition:
  50 //
  51 //   s_k^T * y_k > 0     [s_k^T * B_{k+1} * s_k = s_k^T * y_k > 0]
  52 //
  53 // This condition would always be satisfied if the function was strictly
  54 // convex, alternatively, it is always satisfied provided that a Wolfe line
  55 // search is used (even if the function is not strictly convex).  See [1]
  56 // (p138) for a proof.
  57 //
  58 // Although Ceres will always use a Wolfe line search when using (L)BFGS,
  59 // practical implementation considerations mean that the line search
  60 // may return a point that satisfies only the Armijo condition, and thus
  61 // could violate the Secant equation.  As such, we will only use a step
  62 // to update the Hessian approximation if:
  63 //
  64 //   s_k^T * y_k > tolerance
  65 //
  66 // It is important that tolerance is very small (and >=0), as otherwise we
  67 // might skip the update too often and fail to capture important curvature
  68 // information in the Hessian.  For example going from 1e-10 -> 1e-14 improves
  69 // the NIST benchmark score from 43/54 to 53/54.
  70 //
  71 // [1] Nocedal J., Wright S., Numerical Optimization, 2nd Ed. Springer, 1999.
  72 //
  73 // TODO(alexs.mac): Consider using Damped BFGS update instead of
  74 // skipping update.
  75 const double kLBFGSSecantConditionHessianUpdateTolerance = 1e-14;
  76
  77 LowRankInverseHessian::LowRankInverseHessian(
  78     int num_parameters,
  79     int max_num_corrections,
  80     bool use_approximate_eigenvalue_scaling)
  81     : num_parameters_(num_parameters),
  82       max_num_corrections_(max_num_corrections),
  83       use_approximate_eigenvalue_scaling_(use_approximate_eigenvalue_scaling),
  84       approximate_eigenvalue_scale_(1.0),
  85       delta_x_history_(num_parameters, max_num_corrections),
  86       delta_gradient_history_(num_parameters, max_num_corrections),
  87       delta_x_dot_delta_gradient_(max_num_corrections) {
  88 }
  89
  90 bool LowRankInverseHessian::Update(const Vector& delta_x,
  91                                    const Vector& delta_gradient) {
  92   const double delta_x_dot_delta_gradient = delta_x.dot(delta_gradient);
  93   if (delta_x_dot_delta_gradient <=
  94       kLBFGSSecantConditionHessianUpdateTolerance) {
  95     VLOG(2) << "Skipping L-BFGS Update, delta_x_dot_delta_gradient too "
  96             << "small: " << delta_x_dot_delta_gradient << ", tolerance: "
  97             << kLBFGSSecantConditionHessianUpdateTolerance
  98             << " (Secant condition).";
  99     return false;
 100   }
 101
 102
 103   int next = indices_.size();
 104   // Once the size of the list reaches max_num_corrections_, simulate
 105   // a circular buffer by removing the first element of the list and
 106   // making it the next position where the LBFGS history is stored.
 107   if (next == max_num_corrections_) {
 108     next = indices_.front();
 109     indices_.pop_front();
 110   }
 111
 112   indices_.push_back(next);
 113   delta_x_history_.col(next) = delta_x;
 114   delta_gradient_history_.col(next) = delta_gradient;
 115   delta_x_dot_delta_gradient_(next) = delta_x_dot_delta_gradient;
 116   approximate_eigenvalue_scale_ =
 117       delta_x_dot_delta_gradient / delta_gradient.squaredNorm();
 118   return true;
 119 }
 120
 121 void LowRankInverseHessian::RightMultiply(const double* x_ptr,
 122                                           double* y_ptr) const {
 123   ConstVectorRef gradient(x_ptr, num_parameters_);
 124   VectorRef search_direction(y_ptr, num_parameters_);
 125
 126   search_direction = gradient;
 127
 128   const int num_corrections = indices_.size();
 129   Vector alpha(num_corrections);
 130
 131   for (list<int>::const_reverse_iterator it = indices_.rbegin();
 132        it != indices_.rend();
 133        ++it) {
 134     const double alpha_i = delta_x_history_.col(*it).dot(search_direction) /
 135         delta_x_dot_delta_gradient_(*it);
 136     search_direction -= alpha_i * delta_gradient_history_.col(*it);
 137     alpha(*it) = alpha_i;
 138   }
 139
 140   if (use_approximate_eigenvalue_scaling_) {
 141     // Rescale the initial inverse Hessian approximation (H_0) to be iteratively
 142     // updated so that it is of similar 'size' to the true inverse Hessian along
 143     // the most recent search direction.  As shown in [1]:
 144     //
 145     //   \gamma_k = (delta_gradient_{k-1}' * delta_x_{k-1}) /
 146     //              (delta_gradient_{k-1}' * delta_gradient_{k-1})
 147     //
 148     // Satisfies:
 149     //
 150     //   (1 / \lambda_m) <= \gamma_k <= (1 / \lambda_1)
 151     //
 152     // Where \lambda_1 & \lambda_m are the smallest and largest eigenvalues of
 153     // the true Hessian (not the inverse) along the most recent search direction
 154     // respectively.  Thus \gamma is an approximate eigenvalue of the true
 155     // inverse Hessian, and choosing: H_0 = I * \gamma will yield a starting
 156     // point that has a similar scale to the true inverse Hessian.  This
 157     // technique is widely reported to often improve convergence, however this
 158     // is not universally true, particularly if there are errors in the initial
 159     // jacobians, or if there are significant differences in the sensitivity
 160     // of the problem to the parameters (i.e. the range of the magnitudes of
 161     // the components of the gradient is large).
 162     //
 163     // The original origin of this rescaling trick is somewhat unclear, the
 164     // earliest reference appears to be Oren [1], however it is widely discussed
 165     // without specific attributation in various texts including [2] (p143/178).
 166     //
 167     // [1] Oren S.S., Self-scaling variable metric (SSVM) algorithms Part II:
 168     //     Implementation and experiments, Management Science,
 169     //     20(5), 863-874, 1974.
 170     // [2] Nocedal J., Wright S., Numerical Optimization, Springer, 1999.
 171     search_direction *= approximate_eigenvalue_scale_;
 172
 173     VLOG(4) << "Applying approximate_eigenvalue_scale: "
 174             << approximate_eigenvalue_scale_ << " to initial inverse Hessian "
 175             << "approximation.";
 176   }
 177
 178   for (list<int>::const_iterator it = indices_.begin();
 179        it != indices_.end();
 180        ++it) {
 181     const double beta = delta_gradient_history_.col(*it).dot(search_direction) /
 182         delta_x_dot_delta_gradient_(*it);
 183     search_direction += delta_x_history_.col(*it) * (alpha(*it) - beta);
 184   }
 185 }
 186
 187 }  // namespace internal
 188 }  // namespace ceres