Home | History | Annotate | Download | only in ceres 
      1 // Ceres Solver - A fast non-linear least squares minimizer
      2 // Copyright 2010, 2011, 2012 Google Inc. All rights reserved.
      3 // http://code.google.com/p/ceres-solver/
      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 (at) google.com (Sameer Agarwal)
     30 //
     31 // An iterative solver for solving the Schur complement/reduced camera
     32 // linear system that arise in SfM problems.
     33 
     34 #ifndef CERES_INTERNAL_IMPLICIT_SCHUR_COMPLEMENT_H_
     35 #define CERES_INTERNAL_IMPLICIT_SCHUR_COMPLEMENT_H_
     36 
     37 #include "ceres/linear_operator.h"
     38 #include "ceres/partitioned_matrix_view.h"
     39 #include "ceres/internal/eigen.h"
     40 #include "ceres/internal/scoped_ptr.h"
     41 #include "ceres/types.h"
     42 
     43 namespace ceres {
     44 namespace internal {
     45 
     46 class BlockSparseMatrix;
     47 
     48 // This class implements various linear algebraic operations related
     49 // to the Schur complement without explicitly forming it.
     50 //
     51 //
     52 // Given a reactangular linear system Ax = b, where
     53 //
     54 //   A = [E F]
     55 //
     56 // The normal equations are given by
     57 //
     58 //   A'Ax = A'b
     59 //
     60 //  |E'E E'F||y| = |E'b|
     61 //  |F'E F'F||z|   |F'b|
     62 //
     63 // and the Schur complement system is given by
     64 //
     65 //  [F'F - F'E (E'E)^-1 E'F] z = F'b - F'E (E'E)^-1 E'b
     66 //
     67 // Now if we wish to solve Ax = b in the least squares sense, one way
     68 // is to form this Schur complement system and solve it using
     69 // Preconditioned Conjugate Gradients.
     70 //
     71 // The key operation in a conjugate gradient solver is the evaluation of the
     72 // matrix vector product with the Schur complement
     73 //
     74 //   S = F'F - F'E (E'E)^-1 E'F
     75 //
     76 // It is straightforward to see that matrix vector products with S can
     77 // be evaluated without storing S in memory. Instead, given (E'E)^-1
     78 // (which for our purposes is an easily inverted block diagonal
     79 // matrix), it can be done in terms of matrix vector products with E,
     80 // F and (E'E)^-1. This class implements this functionality and other
     81 // auxilliary bits needed to implement a CG solver on the Schur
     82 // complement using the PartitionedMatrixView object.
     83 //
     84 // THREAD SAFETY: This class is nqot thread safe. In particular, the
     85 // RightMultiply (and the LeftMultiply) methods are not thread safe as
     86 // they depend on mutable arrays used for the temporaries needed to
     87 // compute the product y += Sx;
     88 class ImplicitSchurComplement : public LinearOperator {
     89  public:
     90   // num_eliminate_blocks is the number of E blocks in the matrix
     91   // A.
     92   //
     93   // preconditioner indicates whether the inverse of the matrix F'F
     94   // should be computed or not as a preconditioner for the Schur
     95   // Complement.
     96   //
     97   // TODO(sameeragarwal): Get rid of the two bools below and replace
     98   // them with enums.
     99   ImplicitSchurComplement(int num_eliminate_blocks, bool preconditioner);
    100   virtual ~ImplicitSchurComplement();
    101 
    102   // Initialize the Schur complement for a linear least squares
    103   // problem of the form
    104   //
    105   //   |A      | x = |b|
    106   //   |diag(D)|     |0|
    107   //
    108   // If D is null, then it is treated as a zero dimensional matrix. It
    109   // is important that the matrix A have a BlockStructure object
    110   // associated with it and has a block structure that is compatible
    111   // with the SchurComplement solver.
    112   void Init(const BlockSparseMatrix& A, const double* D, const double* b);
    113 
    114   // y += Sx, where S is the Schur complement.
    115   virtual void RightMultiply(const double* x, double* y) const;
    116 
    117   // The Schur complement is a symmetric positive definite matrix,
    118   // thus the left and right multiply operators are the same.
    119   virtual void LeftMultiply(const double* x, double* y) const {
    120     RightMultiply(x, y);
    121   }
    122 
    123   // y = (E'E)^-1 (E'b - E'F x). Given an estimate of the solution to
    124   // the Schur complement system, this method computes the value of
    125   // the e_block variables that were eliminated to form the Schur
    126   // complement.
    127   void BackSubstitute(const double* x, double* y);
    128 
    129   virtual int num_rows() const { return A_->num_cols_f(); }
    130   virtual int num_cols() const { return A_->num_cols_f(); }
    131   const Vector& rhs()    const { return rhs_;             }
    132 
    133   const BlockSparseMatrix* block_diagonal_EtE_inverse() const {
    134     return block_diagonal_EtE_inverse_.get();
    135   }
    136 
    137   const BlockSparseMatrix* block_diagonal_FtF_inverse() const {
    138     return block_diagonal_FtF_inverse_.get();
    139   }
    140 
    141  private:
    142   void AddDiagonalAndInvert(const double* D, BlockSparseMatrix* matrix);
    143   void UpdateRhs();
    144 
    145   int num_eliminate_blocks_;
    146   bool preconditioner_;
    147 
    148   scoped_ptr<PartitionedMatrixView> A_;
    149   const double* D_;
    150   const double* b_;
    151 
    152   scoped_ptr<BlockSparseMatrix> block_diagonal_EtE_inverse_;
    153   scoped_ptr<BlockSparseMatrix> block_diagonal_FtF_inverse_;
    154 
    155   Vector rhs_;
    156 
    157   // Temporary storage vectors used to implement RightMultiply.
    158   mutable Vector tmp_rows_;
    159   mutable Vector tmp_e_cols_;
    160   mutable Vector tmp_e_cols_2_;
    161   mutable Vector tmp_f_cols_;
    162 };
    163 
    164 }  // namespace internal
    165 }  // namespace ceres
    166 
    167 #endif  // CERES_INTERNAL_IMPLICIT_SCHUR_COMPLEMENT_H_
    168