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 "ceres/internal/eigen.h"
32#include "ceres/low_rank_inverse_hessian.h"
33#include "glog/logging.h"
34
35namespace ceres {
36namespace internal {
37
38LowRankInverseHessian::LowRankInverseHessian(
39    int num_parameters,
40    int max_num_corrections,
41    bool use_approximate_eigenvalue_scaling)
42    : num_parameters_(num_parameters),
43      max_num_corrections_(max_num_corrections),
44      use_approximate_eigenvalue_scaling_(use_approximate_eigenvalue_scaling),
45      num_corrections_(0),
46      approximate_eigenvalue_scale_(1.0),
47      delta_x_history_(num_parameters, max_num_corrections),
48      delta_gradient_history_(num_parameters, max_num_corrections),
49      delta_x_dot_delta_gradient_(max_num_corrections) {
50}
51
52bool LowRankInverseHessian::Update(const Vector& delta_x,
53                                   const Vector& delta_gradient) {
54  const double delta_x_dot_delta_gradient = delta_x.dot(delta_gradient);
55  if (delta_x_dot_delta_gradient <= 1e-10) {
56    VLOG(2) << "Skipping LBFGS Update, delta_x_dot_delta_gradient too small: "
57            << delta_x_dot_delta_gradient;
58    return false;
59  }
60
61  if (num_corrections_ == max_num_corrections_) {
62    // TODO(sameeragarwal): This can be done more efficiently using
63    // a circular buffer/indexing scheme, but for simplicity we will
64    // do the expensive copy for now.
65    delta_x_history_.block(0, 0, num_parameters_, max_num_corrections_ - 1) =
66        delta_x_history_
67        .block(0, 1, num_parameters_, max_num_corrections_ - 1);
68
69    delta_gradient_history_
70        .block(0, 0, num_parameters_, max_num_corrections_ - 1) =
71        delta_gradient_history_
72        .block(0, 1, num_parameters_, max_num_corrections_ - 1);
73
74    delta_x_dot_delta_gradient_.head(num_corrections_ - 1) =
75        delta_x_dot_delta_gradient_.tail(num_corrections_ - 1);
76  } else {
77    ++num_corrections_;
78  }
79
80  delta_x_history_.col(num_corrections_ - 1) = delta_x;
81  delta_gradient_history_.col(num_corrections_ - 1) = delta_gradient;
82  delta_x_dot_delta_gradient_(num_corrections_ - 1) =
83      delta_x_dot_delta_gradient;
84  approximate_eigenvalue_scale_ =
85      delta_x_dot_delta_gradient / delta_gradient.squaredNorm();
86  return true;
87}
88
89void LowRankInverseHessian::RightMultiply(const double* x_ptr,
90                                          double* y_ptr) const {
91  ConstVectorRef gradient(x_ptr, num_parameters_);
92  VectorRef search_direction(y_ptr, num_parameters_);
93
94  search_direction = gradient;
95
96  Vector alpha(num_corrections_);
97
98  for (int i = num_corrections_ - 1; i >= 0; --i) {
99    alpha(i) = delta_x_history_.col(i).dot(search_direction) /
100        delta_x_dot_delta_gradient_(i);
101    search_direction -= alpha(i) * delta_gradient_history_.col(i);
102  }
103
104  if (use_approximate_eigenvalue_scaling_) {
105    // Rescale the initial inverse Hessian approximation (H_0) to be iteratively
106    // updated so that it is of similar 'size' to the true inverse Hessian along
107    // the most recent search direction.  As shown in [1]:
108    //
109    //   \gamma_k = (delta_gradient_{k-1}' * delta_x_{k-1}) /
110    //              (delta_gradient_{k-1}' * delta_gradient_{k-1})
111    //
112    // Satisfies:
113    //
114    //   (1 / \lambda_m) <= \gamma_k <= (1 / \lambda_1)
115    //
116    // Where \lambda_1 & \lambda_m are the smallest and largest eigenvalues of
117    // the true Hessian (not the inverse) along the most recent search direction
118    // respectively.  Thus \gamma is an approximate eigenvalue of the true
119    // inverse Hessian, and choosing: H_0 = I * \gamma will yield a starting
120    // point that has a similar scale to the true inverse Hessian.  This
121    // technique is widely reported to often improve convergence, however this
122    // is not universally true, particularly if there are errors in the initial
123    // jacobians, or if there are significant differences in the sensitivity
124    // of the problem to the parameters (i.e. the range of the magnitudes of
125    // the components of the gradient is large).
126    //
127    // The original origin of this rescaling trick is somewhat unclear, the
128    // earliest reference appears to be Oren [1], however it is widely discussed
129    // without specific attributation in various texts including [2] (p143/178).
130    //
131    // [1] Oren S.S., Self-scaling variable metric (SSVM) algorithms Part II:
132    //     Implementation and experiments, Management Science,
133    //     20(5), 863-874, 1974.
134    // [2] Nocedal J., Wright S., Numerical Optimization, Springer, 1999.
135    search_direction *= approximate_eigenvalue_scale_;
136  }
137
138  for (int i = 0; i < num_corrections_; ++i) {
139    const double beta = delta_gradient_history_.col(i).dot(search_direction) /
140        delta_x_dot_delta_gradient_(i);
141    search_direction += delta_x_history_.col(i) * (alpha(i) - beta);
142  }
143}
144
145}  // namespace internal
146}  // namespace ceres
147