preconditioner.h 8.7 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223
  1. // Ceres Solver - A fast non-linear least squares minimizer
  2. // Copyright 2023 Google Inc. All rights reserved.
  3. // http://ceres-solver.org/
  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. #ifndef CERES_INTERNAL_PRECONDITIONER_H_
  31. #define CERES_INTERNAL_PRECONDITIONER_H_
  32. #include <vector>
  33. #include "ceres/casts.h"
  34. #include "ceres/compressed_row_sparse_matrix.h"
  35. #include "ceres/context_impl.h"
  36. #include "ceres/internal/disable_warnings.h"
  37. #include "ceres/internal/export.h"
  38. #include "ceres/linear_operator.h"
  39. #include "ceres/linear_solver.h"
  40. #include "ceres/sparse_matrix.h"
  41. #include "ceres/types.h"
  42. namespace ceres::internal {
  43. class BlockSparseMatrix;
  44. class SparseMatrix;
  45. class CERES_NO_EXPORT Preconditioner : public LinearOperator {
  46. public:
  47. struct Options {
  48. Options() = default;
  49. Options(const LinearSolver::Options& linear_solver_options)
  50. : type(linear_solver_options.preconditioner_type),
  51. visibility_clustering_type(
  52. linear_solver_options.visibility_clustering_type),
  53. sparse_linear_algebra_library_type(
  54. linear_solver_options.sparse_linear_algebra_library_type),
  55. num_threads(linear_solver_options.num_threads),
  56. row_block_size(linear_solver_options.row_block_size),
  57. e_block_size(linear_solver_options.e_block_size),
  58. f_block_size(linear_solver_options.f_block_size),
  59. elimination_groups(linear_solver_options.elimination_groups),
  60. context(linear_solver_options.context) {}
  61. PreconditionerType type = JACOBI;
  62. VisibilityClusteringType visibility_clustering_type = CANONICAL_VIEWS;
  63. SparseLinearAlgebraLibraryType sparse_linear_algebra_library_type =
  64. SUITE_SPARSE;
  65. OrderingType ordering_type = OrderingType::NATURAL;
  66. // When using the subset preconditioner, all row blocks starting
  67. // from this row block are used to construct the preconditioner.
  68. //
  69. // i.e., the Jacobian matrix A is horizontally partitioned as
  70. //
  71. // A = [P]
  72. // [Q]
  73. //
  74. // where P has subset_preconditioner_start_row_block row blocks,
  75. // and the preconditioner is the inverse of the matrix Q'Q.
  76. int subset_preconditioner_start_row_block = -1;
  77. // If possible, how many threads the preconditioner can use.
  78. int num_threads = 1;
  79. // Hints about the order in which the parameter blocks should be
  80. // eliminated by the linear solver.
  81. //
  82. // For example if elimination_groups is a vector of size k, then
  83. // the linear solver is informed that it should eliminate the
  84. // parameter blocks 0 ... elimination_groups[0] - 1 first, and
  85. // then elimination_groups[0] ... elimination_groups[1] - 1 and so
  86. // on. Within each elimination group, the linear solver is free to
  87. // choose how the parameter blocks are ordered. Different linear
  88. // solvers have differing requirements on elimination_groups.
  89. //
  90. // The most common use is for Schur type solvers, where there
  91. // should be at least two elimination groups and the first
  92. // elimination group must form an independent set in the normal
  93. // equations. The first elimination group corresponds to the
  94. // num_eliminate_blocks in the Schur type solvers.
  95. std::vector<int> elimination_groups;
  96. // If the block sizes in a BlockSparseMatrix are fixed, then in
  97. // some cases the Schur complement based solvers can detect and
  98. // specialize on them.
  99. //
  100. // It is expected that these parameters are set programmatically
  101. // rather than manually.
  102. //
  103. // Please see schur_complement_solver.h and schur_eliminator.h for
  104. // more details.
  105. int row_block_size = Eigen::Dynamic;
  106. int e_block_size = Eigen::Dynamic;
  107. int f_block_size = Eigen::Dynamic;
  108. ContextImpl* context = nullptr;
  109. };
  110. // If the optimization problem is such that there are no remaining
  111. // e-blocks, ITERATIVE_SCHUR with a Schur type preconditioner cannot
  112. // be used. This function returns JACOBI if a preconditioner for
  113. // ITERATIVE_SCHUR is used. The input preconditioner_type is
  114. // returned otherwise.
  115. static PreconditionerType PreconditionerForZeroEBlocks(
  116. PreconditionerType preconditioner_type);
  117. ~Preconditioner() override;
  118. // Update the numerical value of the preconditioner for the linear
  119. // system:
  120. //
  121. // | A | x = |b|
  122. // |diag(D)| |0|
  123. //
  124. // for some vector b. It is important that the matrix A have the
  125. // same block structure as the one used to construct this object.
  126. //
  127. // D can be nullptr, in which case its interpreted as a diagonal matrix
  128. // of size zero.
  129. virtual bool Update(const LinearOperator& A, const double* D) = 0;
  130. // LinearOperator interface. Since the operator is symmetric,
  131. // LeftMultiplyAndAccumulate and num_cols are just calls to
  132. // RightMultiplyAndAccumulate and num_rows respectively. Update() must be
  133. // called before RightMultiplyAndAccumulate can be called.
  134. void RightMultiplyAndAccumulate(const double* x,
  135. double* y) const override = 0;
  136. void LeftMultiplyAndAccumulate(const double* x, double* y) const override {
  137. return RightMultiplyAndAccumulate(x, y);
  138. }
  139. int num_rows() const override = 0;
  140. int num_cols() const override { return num_rows(); }
  141. };
  142. class CERES_NO_EXPORT IdentityPreconditioner : public Preconditioner {
  143. public:
  144. IdentityPreconditioner(int num_rows) : num_rows_(num_rows) {}
  145. bool Update(const LinearOperator& /*A*/, const double* /*D*/) final {
  146. return true;
  147. }
  148. void RightMultiplyAndAccumulate(const double* x, double* y) const final {
  149. VectorRef(y, num_rows_) += ConstVectorRef(x, num_rows_);
  150. }
  151. int num_rows() const final { return num_rows_; }
  152. private:
  153. int num_rows_ = -1;
  154. };
  155. // This templated subclass of Preconditioner serves as a base class for
  156. // other preconditioners that depend on the particular matrix layout of
  157. // the underlying linear operator.
  158. template <typename MatrixType>
  159. class CERES_NO_EXPORT TypedPreconditioner : public Preconditioner {
  160. public:
  161. bool Update(const LinearOperator& A, const double* D) final {
  162. return UpdateImpl(*down_cast<const MatrixType*>(&A), D);
  163. }
  164. private:
  165. virtual bool UpdateImpl(const MatrixType& A, const double* D) = 0;
  166. };
  167. // Preconditioners that depend on access to the low level structure
  168. // of a SparseMatrix.
  169. // clang-format off
  170. using SparseMatrixPreconditioner = TypedPreconditioner<SparseMatrix>;
  171. using BlockSparseMatrixPreconditioner = TypedPreconditioner<BlockSparseMatrix>;
  172. using CompressedRowSparseMatrixPreconditioner = TypedPreconditioner<CompressedRowSparseMatrix>;
  173. // clang-format on
  174. // Wrap a SparseMatrix object as a preconditioner.
  175. class CERES_NO_EXPORT SparseMatrixPreconditionerWrapper final
  176. : public SparseMatrixPreconditioner {
  177. public:
  178. // Wrapper does NOT take ownership of the matrix pointer.
  179. explicit SparseMatrixPreconditionerWrapper(
  180. const SparseMatrix* matrix, const Preconditioner::Options& options);
  181. ~SparseMatrixPreconditionerWrapper() override;
  182. // Preconditioner interface
  183. void RightMultiplyAndAccumulate(const double* x, double* y) const override;
  184. int num_rows() const override;
  185. private:
  186. bool UpdateImpl(const SparseMatrix& A, const double* D) override;
  187. const SparseMatrix* matrix_;
  188. const Preconditioner::Options options_;
  189. };
  190. } // namespace ceres::internal
  191. #include "ceres/internal/reenable_warnings.h"
  192. #endif // CERES_INTERNAL_PRECONDITIONER_H_