blob: 1c6c9925f1c45b12799fbe15b03cf4feafe75123 [file] [log] [blame]
Sameer Agarwal3e8d1922012-11-28 17:20:22 -08001// Ceres Solver - A fast non-linear least squares minimizer
Keir Mierle7492b0d2015-03-17 22:30:16 -07002// Copyright 2015 Google Inc. All rights reserved.
3// http://ceres-solver.org/
Sameer Agarwal3e8d1922012-11-28 17:20:22 -08004//
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
Sameer Agarwal80a53ee2014-01-09 12:40:54 -080031#include <list>
32
Sameer Agarwal3e8d1922012-11-28 17:20:22 -080033#include "ceres/internal/eigen.h"
Sameer Agarwal9883fc32012-11-30 12:32:43 -080034#include "ceres/low_rank_inverse_hessian.h"
35#include "glog/logging.h"
Sameer Agarwal3e8d1922012-11-28 17:20:22 -080036
37namespace ceres {
38namespace internal {
39
Sameer Agarwalbcc865f2014-12-17 07:35:09 -080040using std::list;
41
Alex Stewart3fca2c42013-11-18 10:26:49 +000042// 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//
Sameer Agarwal89a592f2013-11-26 11:35:49 -080073// TODO(alexs.mac): Consider using Damped BFGS update instead of
74// skipping update.
Alex Stewart3fca2c42013-11-18 10:26:49 +000075const double kLBFGSSecantConditionHessianUpdateTolerance = 1e-14;
76
Alex Stewart9aa0e3c2013-07-05 20:22:37 +010077LowRankInverseHessian::LowRankInverseHessian(
78 int num_parameters,
79 int max_num_corrections,
80 bool use_approximate_eigenvalue_scaling)
Sameer Agarwal3e8d1922012-11-28 17:20:22 -080081 : num_parameters_(num_parameters),
82 max_num_corrections_(max_num_corrections),
Alex Stewart9aa0e3c2013-07-05 20:22:37 +010083 use_approximate_eigenvalue_scaling_(use_approximate_eigenvalue_scaling),
Alex Stewart9aa0e3c2013-07-05 20:22:37 +010084 approximate_eigenvalue_scale_(1.0),
Sameer Agarwal3e8d1922012-11-28 17:20:22 -080085 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
Sameer Agarwal9883fc32012-11-30 12:32:43 -080090bool LowRankInverseHessian::Update(const Vector& delta_x,
91 const Vector& delta_gradient) {
Sameer Agarwal3e8d1922012-11-28 17:20:22 -080092 const double delta_x_dot_delta_gradient = delta_x.dot(delta_gradient);
Alex Stewart3fca2c42013-11-18 10:26:49 +000093 if (delta_x_dot_delta_gradient <=
94 kLBFGSSecantConditionHessianUpdateTolerance) {
Alex Stewart331ff092013-11-25 13:44:53 +000095 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).";
Sameer Agarwal3e8d1922012-11-28 17:20:22 -080099 return false;
100 }
101
Sameer Agarwal3e8d1922012-11-28 17:20:22 -0800102
Sameer Agarwal80a53ee2014-01-09 12:40:54 -0800103 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();
Sameer Agarwal3e8d1922012-11-28 17:20:22 -0800110 }
111
Sameer Agarwal80a53ee2014-01-09 12:40:54 -0800112 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;
Alex Stewart9aa0e3c2013-07-05 20:22:37 +0100116 approximate_eigenvalue_scale_ =
117 delta_x_dot_delta_gradient / delta_gradient.squaredNorm();
Sameer Agarwal3e8d1922012-11-28 17:20:22 -0800118 return true;
119}
120
Sameer Agarwal9883fc32012-11-30 12:32:43 -0800121void LowRankInverseHessian::RightMultiply(const double* x_ptr,
122 double* y_ptr) const {
Sameer Agarwal3e8d1922012-11-28 17:20:22 -0800123 ConstVectorRef gradient(x_ptr, num_parameters_);
124 VectorRef search_direction(y_ptr, num_parameters_);
125
126 search_direction = gradient;
127
Sameer Agarwal80a53ee2014-01-09 12:40:54 -0800128 const int num_corrections = indices_.size();
129 Vector alpha(num_corrections);
Sameer Agarwal3e8d1922012-11-28 17:20:22 -0800130
Sameer Agarwalbcc865f2014-12-17 07:35:09 -0800131 for (list<int>::const_reverse_iterator it = indices_.rbegin();
Sameer Agarwal80a53ee2014-01-09 12:40:54 -0800132 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;
Sameer Agarwal3e8d1922012-11-28 17:20:22 -0800138 }
139
Alex Stewart9aa0e3c2013-07-05 20:22:37 +0100140 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_;
Alex Stewart40ef9032013-11-25 16:36:40 +0000172
173 VLOG(4) << "Applying approximate_eigenvalue_scale: "
174 << approximate_eigenvalue_scale_ << " to initial inverse Hessian "
175 << "approximation.";
Alex Stewart9aa0e3c2013-07-05 20:22:37 +0100176 }
Sameer Agarwal3e8d1922012-11-28 17:20:22 -0800177
Sameer Agarwalbcc865f2014-12-17 07:35:09 -0800178 for (list<int>::const_iterator it = indices_.begin();
Sameer Agarwal80a53ee2014-01-09 12:40:54 -0800179 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);
Sameer Agarwal3e8d1922012-11-28 17:20:22 -0800184 }
185}
186
187} // namespace internal
188} // namespace ceres