10ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong// Ceres Solver - A fast non-linear least squares minimizer 20ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong// Copyright 2010, 2011, 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// Abstract interface for objects solving linear systems of various 320ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong// kinds. 330ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong 340ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong#ifndef CERES_INTERNAL_LINEAR_SOLVER_H_ 350ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong#define CERES_INTERNAL_LINEAR_SOLVER_H_ 360ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong 370ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong#include <cstddef> 381d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling#include <map> 391d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling#include <string> 400ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong#include <vector> 410ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong#include "ceres/block_sparse_matrix.h" 420ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong#include "ceres/casts.h" 430ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong#include "ceres/compressed_row_sparse_matrix.h" 440ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong#include "ceres/dense_sparse_matrix.h" 451d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling#include "ceres/execution_summary.h" 460ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong#include "ceres/triplet_sparse_matrix.h" 470ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong#include "ceres/types.h" 481d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling#include "glog/logging.h" 490ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong 500ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kongnamespace ceres { 510ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kongnamespace internal { 520ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong 5379397c21138f54fcff6ec067b44b847f1f7e0e98Carlos Hernandezenum LinearSolverTerminationType { 5479397c21138f54fcff6ec067b44b847f1f7e0e98Carlos Hernandez // Termination criterion was met. 5579397c21138f54fcff6ec067b44b847f1f7e0e98Carlos Hernandez LINEAR_SOLVER_SUCCESS, 5679397c21138f54fcff6ec067b44b847f1f7e0e98Carlos Hernandez 5779397c21138f54fcff6ec067b44b847f1f7e0e98Carlos Hernandez // Solver ran for max_num_iterations and terminated before the 5879397c21138f54fcff6ec067b44b847f1f7e0e98Carlos Hernandez // termination tolerance could be satisfied. 5979397c21138f54fcff6ec067b44b847f1f7e0e98Carlos Hernandez LINEAR_SOLVER_NO_CONVERGENCE, 6079397c21138f54fcff6ec067b44b847f1f7e0e98Carlos Hernandez 6179397c21138f54fcff6ec067b44b847f1f7e0e98Carlos Hernandez // Solver was terminated due to numerical problems, generally due to 6279397c21138f54fcff6ec067b44b847f1f7e0e98Carlos Hernandez // the linear system being poorly conditioned. 6379397c21138f54fcff6ec067b44b847f1f7e0e98Carlos Hernandez LINEAR_SOLVER_FAILURE, 6479397c21138f54fcff6ec067b44b847f1f7e0e98Carlos Hernandez 6579397c21138f54fcff6ec067b44b847f1f7e0e98Carlos Hernandez // Solver failed with a fatal error that cannot be recovered from, 6679397c21138f54fcff6ec067b44b847f1f7e0e98Carlos Hernandez // e.g. CHOLMOD ran out of memory when computing the symbolic or 6779397c21138f54fcff6ec067b44b847f1f7e0e98Carlos Hernandez // numeric factorization or an underlying library was called with 6879397c21138f54fcff6ec067b44b847f1f7e0e98Carlos Hernandez // the wrong arguments. 6979397c21138f54fcff6ec067b44b847f1f7e0e98Carlos Hernandez LINEAR_SOLVER_FATAL_ERROR 7079397c21138f54fcff6ec067b44b847f1f7e0e98Carlos Hernandez}; 7179397c21138f54fcff6ec067b44b847f1f7e0e98Carlos Hernandez 7279397c21138f54fcff6ec067b44b847f1f7e0e98Carlos Hernandez 730ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kongclass LinearOperator; 740ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong 750ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong// Abstract base class for objects that implement algorithms for 760ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong// solving linear systems 770ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong// 780ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong// Ax = b 790ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong// 800ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong// It is expected that a single instance of a LinearSolver object 810ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong// maybe used multiple times for solving multiple linear systems with 820ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong// the same sparsity structure. This allows them to cache and reuse 830ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong// information across solves. This means that calling Solve on the 840ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong// same LinearSolver instance with two different linear systems will 850ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong// result in undefined behaviour. 860ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong// 870ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong// Subclasses of LinearSolver use two structs to configure themselves. 880ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong// The Options struct configures the LinearSolver object for its 890ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong// lifetime. The PerSolveOptions struct is used to specify options for 900ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong// a particular Solve call. 910ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kongclass LinearSolver { 920ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong public: 930ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong struct Options { 940ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong Options() 950ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong : type(SPARSE_NORMAL_CHOLESKY), 960ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong preconditioner_type(JACOBI), 9779397c21138f54fcff6ec067b44b847f1f7e0e98Carlos Hernandez visibility_clustering_type(CANONICAL_VIEWS), 98399f7d09e0c45af54b77b4ab9508d6f23759b927Scott Ettinger dense_linear_algebra_library_type(EIGEN), 99399f7d09e0c45af54b77b4ab9508d6f23759b927Scott Ettinger sparse_linear_algebra_library_type(SUITE_SPARSE), 1001d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling use_postordering(false), 10179397c21138f54fcff6ec067b44b847f1f7e0e98Carlos Hernandez dynamic_sparsity(false), 1020ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong min_num_iterations(1), 1030ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong max_num_iterations(1), 1040ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong num_threads(1), 1050ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong residual_reset_period(10), 1061d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling row_block_size(Eigen::Dynamic), 1071d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling e_block_size(Eigen::Dynamic), 1081d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling f_block_size(Eigen::Dynamic) { 1090ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong } 1100ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong 1110ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong LinearSolverType type; 1120ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong PreconditionerType preconditioner_type; 11379397c21138f54fcff6ec067b44b847f1f7e0e98Carlos Hernandez VisibilityClusteringType visibility_clustering_type; 114399f7d09e0c45af54b77b4ab9508d6f23759b927Scott Ettinger DenseLinearAlgebraLibraryType dense_linear_algebra_library_type; 115399f7d09e0c45af54b77b4ab9508d6f23759b927Scott Ettinger SparseLinearAlgebraLibraryType sparse_linear_algebra_library_type; 1160ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong 1171d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling // See solver.h for information about this flag. 1181d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling bool use_postordering; 11979397c21138f54fcff6ec067b44b847f1f7e0e98Carlos Hernandez bool dynamic_sparsity; 1200ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong 1210ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong // Number of internal iterations that the solver uses. This 1220ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong // parameter only makes sense for iterative solvers like CG. 1230ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong int min_num_iterations; 1240ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong int max_num_iterations; 1250ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong 1260ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong // If possible, how many threads can the solver use. 1270ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong int num_threads; 1280ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong 1290ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong // Hints about the order in which the parameter blocks should be 1300ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong // eliminated by the linear solver. 1310ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong // 1320ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong // For example if elimination_groups is a vector of size k, then 1330ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong // the linear solver is informed that it should eliminate the 1341d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling // parameter blocks 0 ... elimination_groups[0] - 1 first, and 1351d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling // then elimination_groups[0] ... elimination_groups[1] - 1 and so 1361d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling // on. Within each elimination group, the linear solver is free to 1371d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling // choose how the parameter blocks are ordered. Different linear 1381d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling // solvers have differing requirements on elimination_groups. 1390ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong // 1400ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong // The most common use is for Schur type solvers, where there 1410ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong // should be at least two elimination groups and the first 1420ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong // elimination group must form an independent set in the normal 1430ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong // equations. The first elimination group corresponds to the 1440ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong // num_eliminate_blocks in the Schur type solvers. 1450ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong vector<int> elimination_groups; 1460ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong 1470ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong // Iterative solvers, e.g. Preconditioned Conjugate Gradients 1480ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong // maintain a cheap estimate of the residual which may become 1490ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong // inaccurate over time. Thus for non-zero values of this 1500ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong // parameter, the solver can be told to recalculate the value of 1510ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong // the residual using a |b - Ax| evaluation. 1520ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong int residual_reset_period; 1530ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong 1540ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong // If the block sizes in a BlockSparseMatrix are fixed, then in 1550ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong // some cases the Schur complement based solvers can detect and 1560ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong // specialize on them. 1570ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong // 1580ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong // It is expected that these parameters are set programmatically 1590ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong // rather than manually. 1600ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong // 1610ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong // Please see schur_complement_solver.h and schur_eliminator.h for 1620ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong // more details. 1630ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong int row_block_size; 1640ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong int e_block_size; 1650ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong int f_block_size; 1660ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong }; 1670ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong 1680ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong // Options for the Solve method. 1690ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong struct PerSolveOptions { 1700ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong PerSolveOptions() 1710ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong : D(NULL), 1720ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong preconditioner(NULL), 1730ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong r_tolerance(0.0), 1740ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong q_tolerance(0.0) { 1750ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong } 1760ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong 1770ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong // This option only makes sense for unsymmetric linear solvers 1780ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong // that can solve rectangular linear systems. 1790ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong // 1800ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong // Given a matrix A, an optional diagonal matrix D as a vector, 1810ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong // and a vector b, the linear solver will solve for 1820ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong // 1830ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong // | A | x = | b | 1840ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong // | D | | 0 | 1850ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong // 1860ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong // If D is null, then it is treated as zero, and the solver returns 1870ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong // the solution to 1880ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong // 1890ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong // A x = b 1900ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong // 1910ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong // In either case, x is the vector that solves the following 1920ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong // optimization problem. 1930ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong // 1940ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong // arg min_x ||Ax - b||^2 + ||Dx||^2 1950ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong // 1960ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong // Here A is a matrix of size m x n, with full column rank. If A 1970ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong // does not have full column rank, the results returned by the 1980ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong // solver cannot be relied on. D, if it is not null is an array of 1990ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong // size n. b is an array of size m and x is an array of size n. 2000ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong double * D; 2010ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong 2020ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong // This option only makes sense for iterative solvers. 2030ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong // 2040ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong // In general the performance of an iterative linear solver 2050ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong // depends on the condition number of the matrix A. For example 2060ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong // the convergence rate of the conjugate gradients algorithm 2070ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong // is proportional to the square root of the condition number. 2080ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong // 2090ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong // One particularly useful technique for improving the 2100ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong // conditioning of a linear system is to precondition it. In its 2110ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong // simplest form a preconditioner is a matrix M such that instead 2120ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong // of solving Ax = b, we solve the linear system AM^{-1} y = b 2130ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong // instead, where M is such that the condition number k(AM^{-1}) 2140ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong // is smaller than the conditioner k(A). Given the solution to 2150ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong // this system, x = M^{-1} y. The iterative solver takes care of 2160ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong // the mechanics of solving the preconditioned system and 2170ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong // returning the corrected solution x. The user only needs to 2180ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong // supply a linear operator. 2190ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong // 2200ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong // A null preconditioner is equivalent to an identity matrix being 2210ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong // used a preconditioner. 2220ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong LinearOperator* preconditioner; 2230ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong 2240ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong 2250ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong // The following tolerance related options only makes sense for 2260ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong // iterative solvers. Direct solvers ignore them. 2270ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong 2280ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong // Solver terminates when 2290ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong // 2300ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong // |Ax - b| <= r_tolerance * |b|. 2310ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong // 2320ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong // This is the most commonly used termination criterion for 2330ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong // iterative solvers. 2340ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong double r_tolerance; 2350ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong 2360ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong // For PSD matrices A, let 2370ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong // 2380ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong // Q(x) = x'Ax - 2b'x 2390ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong // 2400ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong // be the cost of the quadratic function defined by A and b. Then, 2410ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong // the solver terminates at iteration i if 2420ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong // 2430ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong // i * (Q(x_i) - Q(x_i-1)) / Q(x_i) < q_tolerance. 2440ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong // 2450ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong // This termination criterion is more useful when using CG to 2460ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong // solve the Newton step. This particular convergence test comes 2470ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong // from Stephen Nash's work on truncated Newton 2480ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong // methods. References: 2490ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong // 2500ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong // 1. Stephen G. Nash & Ariela Sofer, Assessing A Search 2510ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong // Direction Within A Truncated Newton Method, Operation 2520ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong // Research Letters 9(1990) 219-221. 2530ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong // 2540ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong // 2. Stephen G. Nash, A Survey of Truncated Newton Methods, 2550ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong // Journal of Computational and Applied Mathematics, 2560ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong // 124(1-2), 45-59, 2000. 2570ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong // 2580ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong double q_tolerance; 2590ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong }; 2600ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong 2610ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong // Summary of a call to the Solve method. We should move away from 2620ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong // the true/false method for determining solver success. We should 2630ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong // let the summary object do the talking. 2640ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong struct Summary { 2650ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong Summary() 2660ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong : residual_norm(0.0), 2670ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong num_iterations(-1), 26879397c21138f54fcff6ec067b44b847f1f7e0e98Carlos Hernandez termination_type(LINEAR_SOLVER_FAILURE) { 2690ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong } 2700ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong 2710ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong double residual_norm; 2720ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong int num_iterations; 2730ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong LinearSolverTerminationType termination_type; 27479397c21138f54fcff6ec067b44b847f1f7e0e98Carlos Hernandez string message; 2750ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong }; 2760ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong 27779397c21138f54fcff6ec067b44b847f1f7e0e98Carlos Hernandez // If the optimization problem is such that there are no remaining 27879397c21138f54fcff6ec067b44b847f1f7e0e98Carlos Hernandez // e-blocks, a Schur type linear solver cannot be used. If the 27979397c21138f54fcff6ec067b44b847f1f7e0e98Carlos Hernandez // linear solver is of Schur type, this function implements a policy 28079397c21138f54fcff6ec067b44b847f1f7e0e98Carlos Hernandez // to select an alternate nearest linear solver to the one selected 28179397c21138f54fcff6ec067b44b847f1f7e0e98Carlos Hernandez // by the user. The input linear_solver_type is returned otherwise. 28279397c21138f54fcff6ec067b44b847f1f7e0e98Carlos Hernandez static LinearSolverType LinearSolverForZeroEBlocks( 28379397c21138f54fcff6ec067b44b847f1f7e0e98Carlos Hernandez LinearSolverType linear_solver_type); 28479397c21138f54fcff6ec067b44b847f1f7e0e98Carlos Hernandez 2850ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong virtual ~LinearSolver(); 2860ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong 2870ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong // Solve Ax = b. 2880ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong virtual Summary Solve(LinearOperator* A, 2890ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong const double* b, 2900ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong const PerSolveOptions& per_solve_options, 2910ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong double* x) = 0; 2920ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong 2931d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling // The following two methods return copies instead of references so 2941d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling // that the base class implementation does not have to worry about 2951d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling // life time issues. Further, these calls are not expected to be 2961d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling // frequent or performance sensitive. 2971d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling virtual map<string, int> CallStatistics() const { 2981d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling return map<string, int>(); 2991d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling } 3001d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling 3011d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling virtual map<string, double> TimeStatistics() const { 3021d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling return map<string, double>(); 3031d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling } 3041d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling 3050ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong // Factory 3060ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong static LinearSolver* Create(const Options& options); 3070ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong}; 3080ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong 3090ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong// This templated subclass of LinearSolver serves as a base class for 3100ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong// other linear solvers that depend on the particular matrix layout of 3110ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong// the underlying linear operator. For example some linear solvers 3120ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong// need low level access to the TripletSparseMatrix implementing the 3130ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong// LinearOperator interface. This class hides those implementation 3140ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong// details behind a private virtual method, and has the Solve method 3150ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong// perform the necessary upcasting. 3160ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kongtemplate <typename MatrixType> 3170ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kongclass TypedLinearSolver : public LinearSolver { 3180ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong public: 3190ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong virtual ~TypedLinearSolver() {} 3200ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong virtual LinearSolver::Summary Solve( 3210ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong LinearOperator* A, 3220ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong const double* b, 3230ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong const LinearSolver::PerSolveOptions& per_solve_options, 3240ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong double* x) { 3251d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling ScopedExecutionTimer total_time("LinearSolver::Solve", &execution_summary_); 3260ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong CHECK_NOTNULL(A); 3270ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong CHECK_NOTNULL(b); 3280ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong CHECK_NOTNULL(x); 3290ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong return SolveImpl(down_cast<MatrixType*>(A), b, per_solve_options, x); 3300ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong } 3310ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong 3321d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling virtual map<string, int> CallStatistics() const { 3331d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling return execution_summary_.calls(); 3341d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling } 3351d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling 3361d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling virtual map<string, double> TimeStatistics() const { 3371d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling return execution_summary_.times(); 3381d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling } 3391d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling 3400ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong private: 3410ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong virtual LinearSolver::Summary SolveImpl( 3420ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong MatrixType* A, 3430ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong const double* b, 3440ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong const LinearSolver::PerSolveOptions& per_solve_options, 3450ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong double* x) = 0; 3461d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling 3471d2624a10e2c559f8ba9ef89eaa30832c0a83a96Sascha Haeberling ExecutionSummary execution_summary_; 3480ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong}; 3490ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong 3500ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong// Linear solvers that depend on acccess to the low level structure of 3510ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong// a SparseMatrix. 3520ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kongtypedef TypedLinearSolver<BlockSparseMatrix> BlockSparseMatrixSolver; // NOLINT 3530ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kongtypedef TypedLinearSolver<CompressedRowSparseMatrix> CompressedRowSparseMatrixSolver; // NOLINT 3540ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kongtypedef TypedLinearSolver<DenseSparseMatrix> DenseSparseMatrixSolver; // NOLINT 3550ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kongtypedef TypedLinearSolver<TripletSparseMatrix> TripletSparseMatrixSolver; // NOLINT 3560ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong 3570ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong} // namespace internal 3580ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong} // namespace ceres 3590ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong 3600ae28bd5885b5daa526898fcf7c323dc2c3e1963Angus Kong#endif // CERES_INTERNAL_LINEAR_SOLVER_H_ 361