diff options
Diffstat (limited to 'extern/ceres/internal/ceres/low_rank_inverse_hessian.cc')
| -rw-r--r-- | extern/ceres/internal/ceres/low_rank_inverse_hessian.cc | 188 |
1 files changed, 188 insertions, 0 deletions
diff --git a/extern/ceres/internal/ceres/low_rank_inverse_hessian.cc b/extern/ceres/internal/ceres/low_rank_inverse_hessian.cc new file mode 100644 index 00000000000..1c6c9925f1c --- /dev/null +++ b/extern/ceres/internal/ceres/low_rank_inverse_hessian.cc @@ -0,0 +1,188 @@ +// Ceres Solver - A fast non-linear least squares minimizer +// Copyright 2015 Google Inc. All rights reserved. +// http://ceres-solver.org/ +// +// Redistribution and use in source and binary forms, with or without +// modification, are permitted provided that the following conditions are met: +// +// * Redistributions of source code must retain the above copyright notice, +// this list of conditions and the following disclaimer. +// * Redistributions in binary form must reproduce the above copyright notice, +// this list of conditions and the following disclaimer in the documentation +// and/or other materials provided with the distribution. +// * Neither the name of Google Inc. nor the names of its contributors may be +// used to endorse or promote products derived from this software without +// specific prior written permission. +// +// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" +// AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE +// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE +// ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE +// LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +// CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF +// SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS +// INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +// CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +// ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +// POSSIBILITY OF SUCH DAMAGE. +// +// Author: sameeragarwal@google.com (Sameer Agarwal) + +#include <list> + +#include "ceres/internal/eigen.h" +#include "ceres/low_rank_inverse_hessian.h" +#include "glog/logging.h" + +namespace ceres { +namespace internal { + +using std::list; + +// The (L)BFGS algorithm explicitly requires that the secant equation: +// +// B_{k+1} * s_k = y_k +// +// Is satisfied at each iteration, where B_{k+1} is the approximated +// Hessian at the k+1-th iteration, s_k = (x_{k+1} - x_{k}) and +// y_k = (grad_{k+1} - grad_{k}). As the approximated Hessian must be +// positive definite, this is equivalent to the condition: +// +// s_k^T * y_k > 0 [s_k^T * B_{k+1} * s_k = s_k^T * y_k > 0] +// +// This condition would always be satisfied if the function was strictly +// convex, alternatively, it is always satisfied provided that a Wolfe line +// search is used (even if the function is not strictly convex). See [1] +// (p138) for a proof. +// +// Although Ceres will always use a Wolfe line search when using (L)BFGS, +// practical implementation considerations mean that the line search +// may return a point that satisfies only the Armijo condition, and thus +// could violate the Secant equation. As such, we will only use a step +// to update the Hessian approximation if: +// +// s_k^T * y_k > tolerance +// +// It is important that tolerance is very small (and >=0), as otherwise we +// might skip the update too often and fail to capture important curvature +// information in the Hessian. For example going from 1e-10 -> 1e-14 improves +// the NIST benchmark score from 43/54 to 53/54. +// +// [1] Nocedal J., Wright S., Numerical Optimization, 2nd Ed. Springer, 1999. +// +// TODO(alexs.mac): Consider using Damped BFGS update instead of +// skipping update. +const double kLBFGSSecantConditionHessianUpdateTolerance = 1e-14; + +LowRankInverseHessian::LowRankInverseHessian( + int num_parameters, + int max_num_corrections, + bool use_approximate_eigenvalue_scaling) + : num_parameters_(num_parameters), + max_num_corrections_(max_num_corrections), + use_approximate_eigenvalue_scaling_(use_approximate_eigenvalue_scaling), + approximate_eigenvalue_scale_(1.0), + delta_x_history_(num_parameters, max_num_corrections), + delta_gradient_history_(num_parameters, max_num_corrections), + delta_x_dot_delta_gradient_(max_num_corrections) { +} + +bool LowRankInverseHessian::Update(const Vector& delta_x, + const Vector& delta_gradient) { + const double delta_x_dot_delta_gradient = delta_x.dot(delta_gradient); + if (delta_x_dot_delta_gradient <= + kLBFGSSecantConditionHessianUpdateTolerance) { + VLOG(2) << "Skipping L-BFGS Update, delta_x_dot_delta_gradient too " + << "small: " << delta_x_dot_delta_gradient << ", tolerance: " + << kLBFGSSecantConditionHessianUpdateTolerance + << " (Secant condition)."; + return false; + } + + + int next = indices_.size(); + // Once the size of the list reaches max_num_corrections_, simulate + // a circular buffer by removing the first element of the list and + // making it the next position where the LBFGS history is stored. + if (next == max_num_corrections_) { + next = indices_.front(); + indices_.pop_front(); + } + + indices_.push_back(next); + delta_x_history_.col(next) = delta_x; + delta_gradient_history_.col(next) = delta_gradient; + delta_x_dot_delta_gradient_(next) = delta_x_dot_delta_gradient; + approximate_eigenvalue_scale_ = + delta_x_dot_delta_gradient / delta_gradient.squaredNorm(); + return true; +} + +void LowRankInverseHessian::RightMultiply(const double* x_ptr, + double* y_ptr) const { + ConstVectorRef gradient(x_ptr, num_parameters_); + VectorRef search_direction(y_ptr, num_parameters_); + + search_direction = gradient; + + const int num_corrections = indices_.size(); + Vector alpha(num_corrections); + + for (list<int>::const_reverse_iterator it = indices_.rbegin(); + it != indices_.rend(); + ++it) { + const double alpha_i = delta_x_history_.col(*it).dot(search_direction) / + delta_x_dot_delta_gradient_(*it); + search_direction -= alpha_i * delta_gradient_history_.col(*it); + alpha(*it) = alpha_i; + } + + if (use_approximate_eigenvalue_scaling_) { + // Rescale the initial inverse Hessian approximation (H_0) to be iteratively + // updated so that it is of similar 'size' to the true inverse Hessian along + // the most recent search direction. As shown in [1]: + // + // \gamma_k = (delta_gradient_{k-1}' * delta_x_{k-1}) / + // (delta_gradient_{k-1}' * delta_gradient_{k-1}) + // + // Satisfies: + // + // (1 / \lambda_m) <= \gamma_k <= (1 / \lambda_1) + // + // Where \lambda_1 & \lambda_m are the smallest and largest eigenvalues of + // the true Hessian (not the inverse) along the most recent search direction + // respectively. Thus \gamma is an approximate eigenvalue of the true + // inverse Hessian, and choosing: H_0 = I * \gamma will yield a starting + // point that has a similar scale to the true inverse Hessian. This + // technique is widely reported to often improve convergence, however this + // is not universally true, particularly if there are errors in the initial + // jacobians, or if there are significant differences in the sensitivity + // of the problem to the parameters (i.e. the range of the magnitudes of + // the components of the gradient is large). + // + // The original origin of this rescaling trick is somewhat unclear, the + // earliest reference appears to be Oren [1], however it is widely discussed + // without specific attributation in various texts including [2] (p143/178). + // + // [1] Oren S.S., Self-scaling variable metric (SSVM) algorithms Part II: + // Implementation and experiments, Management Science, + // 20(5), 863-874, 1974. + // [2] Nocedal J., Wright S., Numerical Optimization, Springer, 1999. + search_direction *= approximate_eigenvalue_scale_; + + VLOG(4) << "Applying approximate_eigenvalue_scale: " + << approximate_eigenvalue_scale_ << " to initial inverse Hessian " + << "approximation."; + } + + for (list<int>::const_iterator it = indices_.begin(); + it != indices_.end(); + ++it) { + const double beta = delta_gradient_history_.col(*it).dot(search_direction) / + delta_x_dot_delta_gradient_(*it); + search_direction += delta_x_history_.col(*it) * (alpha(*it) - beta); + } +} + +} // namespace internal +} // namespace ceres |