internal/ceres/preconditioner.h

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