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