1// Ceres Solver - A fast non-linear least squares minimizer
2// Copyright 2012 Google Inc. All rights reserved.
3// http://code.google.com/p/ceres-solver/
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
37namespace ceres {
38namespace internal {
39
40// The (L)BFGS algorithm explicitly requires that the secant equation:
41//
42//   B_{k+1} * s_k = y_k
43//
44// Is satisfied at each iteration, where B_{k+1} is the approximated
45// Hessian at the k+1-th iteration, s_k = (x_{k+1} - x_{k}) and
46// y_k = (grad_{k+1} - grad_{k}). As the approximated Hessian must be
47// positive definite, this is equivalent to the condition:
48//
49//   s_k^T * y_k > 0     [s_k^T * B_{k+1} * s_k = s_k^T * y_k > 0]
50//
51// This condition would always be satisfied if the function was strictly
52// convex, alternatively, it is always satisfied provided that a Wolfe line
53// search is used (even if the function is not strictly convex).  See [1]
54// (p138) for a proof.
55//
56// Although Ceres will always use a Wolfe line search when using (L)BFGS,
57// practical implementation considerations mean that the line search
58// may return a point that satisfies only the Armijo condition, and thus
59// could violate the Secant equation.  As such, we will only use a step
60// to update the Hessian approximation if:
61//
62//   s_k^T * y_k > tolerance
63//
64// It is important that tolerance is very small (and >=0), as otherwise we
65// might skip the update too often and fail to capture important curvature
66// information in the Hessian.  For example going from 1e-10 -> 1e-14 improves
67// the NIST benchmark score from 43/54 to 53/54.
68//
69// [1] Nocedal J., Wright S., Numerical Optimization, 2nd Ed. Springer, 1999.
70//
71// TODO(alexs.mac): Consider using Damped BFGS update instead of
72// skipping update.
73const double kLBFGSSecantConditionHessianUpdateTolerance = 1e-14;
74
75LowRankInverseHessian::LowRankInverseHessian(
76    int num_parameters,
77    int max_num_corrections,
78    bool use_approximate_eigenvalue_scaling)
79    : num_parameters_(num_parameters),
80      max_num_corrections_(max_num_corrections),
81      use_approximate_eigenvalue_scaling_(use_approximate_eigenvalue_scaling),
82      approximate_eigenvalue_scale_(1.0),
83      delta_x_history_(num_parameters, max_num_corrections),
84      delta_gradient_history_(num_parameters, max_num_corrections),
85      delta_x_dot_delta_gradient_(max_num_corrections) {
86}
87
88bool LowRankInverseHessian::Update(const Vector& delta_x,
89                                   const Vector& delta_gradient) {
90  const double delta_x_dot_delta_gradient = delta_x.dot(delta_gradient);
91  if (delta_x_dot_delta_gradient <=
92      kLBFGSSecantConditionHessianUpdateTolerance) {
93    VLOG(2) << "Skipping L-BFGS Update, delta_x_dot_delta_gradient too "
94            << "small: " << delta_x_dot_delta_gradient << ", tolerance: "
95            << kLBFGSSecantConditionHessianUpdateTolerance
96            << " (Secant condition).";
97    return false;
98  }
99
100
101  int next = indices_.size();
102  // Once the size of the list reaches max_num_corrections_, simulate
103  // a circular buffer by removing the first element of the list and
104  // making it the next position where the LBFGS history is stored.
105  if (next == max_num_corrections_) {
106    next = indices_.front();
107    indices_.pop_front();
108  }
109
110  indices_.push_back(next);
111  delta_x_history_.col(next) = delta_x;
112  delta_gradient_history_.col(next) = delta_gradient;
113  delta_x_dot_delta_gradient_(next) = delta_x_dot_delta_gradient;
114  approximate_eigenvalue_scale_ =
115      delta_x_dot_delta_gradient / delta_gradient.squaredNorm();
116  return true;
117}
118
119void LowRankInverseHessian::RightMultiply(const double* x_ptr,
120                                          double* y_ptr) const {
121  ConstVectorRef gradient(x_ptr, num_parameters_);
122  VectorRef search_direction(y_ptr, num_parameters_);
123
124  search_direction = gradient;
125
126  const int num_corrections = indices_.size();
127  Vector alpha(num_corrections);
128
129  for (std::list<int>::const_reverse_iterator it = indices_.rbegin();
130       it != indices_.rend();
131       ++it) {
132    const double alpha_i = delta_x_history_.col(*it).dot(search_direction) /
133        delta_x_dot_delta_gradient_(*it);
134    search_direction -= alpha_i * delta_gradient_history_.col(*it);
135    alpha(*it) = alpha_i;
136  }
137
138  if (use_approximate_eigenvalue_scaling_) {
139    // Rescale the initial inverse Hessian approximation (H_0) to be iteratively
140    // updated so that it is of similar 'size' to the true inverse Hessian along
141    // the most recent search direction.  As shown in [1]:
142    //
143    //   \gamma_k = (delta_gradient_{k-1}' * delta_x_{k-1}) /
144    //              (delta_gradient_{k-1}' * delta_gradient_{k-1})
145    //
146    // Satisfies:
147    //
148    //   (1 / \lambda_m) <= \gamma_k <= (1 / \lambda_1)
149    //
150    // Where \lambda_1 & \lambda_m are the smallest and largest eigenvalues of
151    // the true Hessian (not the inverse) along the most recent search direction
152    // respectively.  Thus \gamma is an approximate eigenvalue of the true
153    // inverse Hessian, and choosing: H_0 = I * \gamma will yield a starting
154    // point that has a similar scale to the true inverse Hessian.  This
155    // technique is widely reported to often improve convergence, however this
156    // is not universally true, particularly if there are errors in the initial
157    // jacobians, or if there are significant differences in the sensitivity
158    // of the problem to the parameters (i.e. the range of the magnitudes of
159    // the components of the gradient is large).
160    //
161    // The original origin of this rescaling trick is somewhat unclear, the
162    // earliest reference appears to be Oren [1], however it is widely discussed
163    // without specific attributation in various texts including [2] (p143/178).
164    //
165    // [1] Oren S.S., Self-scaling variable metric (SSVM) algorithms Part II:
166    //     Implementation and experiments, Management Science,
167    //     20(5), 863-874, 1974.
168    // [2] Nocedal J., Wright S., Numerical Optimization, Springer, 1999.
169    search_direction *= approximate_eigenvalue_scale_;
170
171    VLOG(4) << "Applying approximate_eigenvalue_scale: "
172            << approximate_eigenvalue_scale_ << " to initial inverse Hessian "
173            << "approximation.";
174  }
175
176  for (std::list<int>::const_iterator it = indices_.begin();
177       it != indices_.end();
178       ++it) {
179    const double beta = delta_gradient_history_.col(*it).dot(search_direction) /
180        delta_x_dot_delta_gradient_(*it);
181    search_direction += delta_x_history_.col(*it) * (alpha(*it) - beta);
182  }
183}
184
185}  // namespace internal
186}  // namespace ceres
187