10ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong// Ceres Solver - A fast non-linear least squares minimizer
20ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong// Copyright 2012 Google Inc. All rights reserved.
30ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong// http://code.google.com/p/ceres-solver/
40ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong//
50ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong// Redistribution and use in source and binary forms, with or without
60ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong// modification, are permitted provided that the following conditions are met:
70ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong//
80ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong// * Redistributions of source code must retain the above copyright notice,
90ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong//   this list of conditions and the following disclaimer.
100ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong// * Redistributions in binary form must reproduce the above copyright notice,
110ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong//   this list of conditions and the following disclaimer in the documentation
120ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong//   and/or other materials provided with the distribution.
130ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong// * Neither the name of Google Inc. nor the names of its contributors may be
140ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong//   used to endorse or promote products derived from this software without
150ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong//   specific prior written permission.
160ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong//
170ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
180ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong// AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
190ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
200ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong// ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
210ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong// LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
220ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong// CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
230ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong// SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
240ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong// INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
250ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong// CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
260ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong// ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
270ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong// POSSIBILITY OF SUCH DAMAGE.
280ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong//
290ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong// Author: sameeragarwal@google.com (Sameer Agarwal)
300ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong
310ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong#ifndef CERES_INTERNAL_TRUST_REGION_STRATEGY_H_
320ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong#define CERES_INTERNAL_TRUST_REGION_STRATEGY_H_
330ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong
341d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling#include <string>
351d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling#include "ceres/internal/port.h"
3679397c21138f54fcff6ec067b44b847f1f7e0e98Carlos Hernandez#include "ceres/linear_solver.h"
370ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong
380ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kongnamespace ceres {
390ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kongnamespace internal {
400ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong
410ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kongclass LinearSolver;
420ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kongclass SparseMatrix;
430ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong
440ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong// Interface for classes implementing various trust region strategies
450ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong// for nonlinear least squares problems.
460ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong//
470ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong// The object is expected to maintain and update a trust region
480ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong// radius, which it then uses to solve for the trust region step using
490ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong// the jacobian matrix and residual vector.
500ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong//
510ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong// Here the term trust region radius is used loosely, as the strategy
520ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong// is free to treat it as guidance and violate it as need be. e.g.,
530ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong// the LevenbergMarquardtStrategy uses the inverse of the trust region
540ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong// radius to scale the damping term, which controls the step size, but
550ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong// does not set a hard limit on its size.
560ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kongclass TrustRegionStrategy {
571d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling public:
580ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong  struct Options {
590ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong    Options()
600ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong        : trust_region_strategy_type(LEVENBERG_MARQUARDT),
610ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong          initial_radius(1e4),
620ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong          max_radius(1e32),
631d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling          min_lm_diagonal(1e-6),
641d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling          max_lm_diagonal(1e32),
650ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong          dogleg_type(TRADITIONAL_DOGLEG) {
660ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong    }
670ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong
680ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong    TrustRegionStrategyType trust_region_strategy_type;
690ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong    // Linear solver used for actually solving the trust region step.
700ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong    LinearSolver* linear_solver;
710ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong    double initial_radius;
720ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong    double max_radius;
730ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong
740ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong    // Minimum and maximum values of the diagonal damping matrix used
750ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong    // by LevenbergMarquardtStrategy. The DoglegStrategy also uses
760ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong    // these bounds to construct a regularizing diagonal to ensure
770ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong    // that the Gauss-Newton step computation is of full rank.
781d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling    double min_lm_diagonal;
791d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling    double max_lm_diagonal;
800ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong
810ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong    // Further specify which dogleg method to use
820ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong    DoglegType dogleg_type;
830ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong  };
840ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong
850ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong  // Per solve options.
860ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong  struct PerSolveOptions {
871d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling    PerSolveOptions()
881d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling        : eta(0),
891d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling          dump_filename_base(""),
901d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling          dump_format_type(TEXTFILE) {
911d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling    }
921d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling
930ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong    // Forcing sequence for inexact solves.
940ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong    double eta;
951d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling
961d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling    // If non-empty and dump_format_type is not CONSOLE, the trust
971d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling    // regions strategy will write the linear system to file(s) with
981d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling    // name starting with dump_filename_base.  If dump_format_type is
991d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling    // CONSOLE then dump_filename_base will be ignored and the linear
1001d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling    // system will be written to the standard error.
1011d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling    string dump_filename_base;
1021d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling    DumpFormatType dump_format_type;
1030ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong  };
1040ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong
1050ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong  struct Summary {
1060ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong    Summary()
1070ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong        : residual_norm(0.0),
1080ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong          num_iterations(-1),
10979397c21138f54fcff6ec067b44b847f1f7e0e98Carlos Hernandez          termination_type(LINEAR_SOLVER_FAILURE) {
1100ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong    }
1110ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong
1120ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong    // If the trust region problem is,
1130ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong    //
1140ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong    //   1/2 x'Ax + b'x + c,
1150ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong    //
1160ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong    // then
1170ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong    //
1180ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong    //   residual_norm = |Ax -b|
1190ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong    double residual_norm;
1200ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong
1210ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong    // Number of iterations used by the linear solver. If a linear
1220ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong    // solver was not called (e.g., DogLegStrategy after an
1230ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong    // unsuccessful step), then this would be zero.
1240ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong    int num_iterations;
1250ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong
1260ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong    // Status of the linear solver used to solve the Newton system.
1270ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong    LinearSolverTerminationType termination_type;
1280ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong  };
1290ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong
1300ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong  virtual ~TrustRegionStrategy();
1310ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong
1320ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong  // Use the current radius to solve for the trust region step.
1330ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong  virtual Summary ComputeStep(const PerSolveOptions& per_solve_options,
1340ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong                              SparseMatrix* jacobian,
1350ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong                              const double* residuals,
1360ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong                              double* step) = 0;
1370ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong
1380ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong  // Inform the strategy that the current step has been accepted, and
1390ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong  // that the ratio of the decrease in the non-linear objective to the
1400ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong  // decrease in the trust region model is step_quality.
1410ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong  virtual void StepAccepted(double step_quality) = 0;
1420ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong
1430ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong  // Inform the strategy that the current step has been rejected, and
1440ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong  // that the ratio of the decrease in the non-linear objective to the
1450ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong  // decrease in the trust region model is step_quality.
1460ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong  virtual void StepRejected(double step_quality) = 0;
1470ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong
1480ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong  // Inform the strategy that the current step has been rejected
1490ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong  // because it was found to be numerically invalid.
1500ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong  // StepRejected/StepAccepted will not be called for this step, and
1510ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong  // the strategy is free to do what it wants with this information.
1520ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong  virtual void StepIsInvalid() = 0;
1530ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong
1540ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong  // Current trust region radius.
1550ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong  virtual double Radius() const = 0;
1560ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong
1570ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong  // Factory.
1580ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong  static TrustRegionStrategy* Create(const Options& options);
1590ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong};
1600ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong
1610ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong}  // namespace internal
1620ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong}  // namespace ceres
1630ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong
1640ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong#endif  // CERES_INTERNAL_TRUST_REGION_STRATEGY_H_
165