internal/ceres/inner_product_computer.h

   1 // Ceres Solver - A fast non-linear least squares minimizer
   2 // Copyright 2017 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_INNER_PRODUCT_COMPUTER_H_
  32 #define CERES_INTERNAL_INNER_PRODUCT_COMPUTER_H_
  33
  34 #include <vector>
  35
  36 #include "ceres/block_sparse_matrix.h"
  37 #include "ceres/compressed_row_sparse_matrix.h"
  38 #include "ceres/internal/scoped_ptr.h"
  39
  40 namespace ceres {
  41 namespace internal {
  42
  43 // This class is used to repeatedly compute the inner product
  44 //
  45 //   result = m' * m
  46 //
  47 // where the sparsity structure of m remains constant across calls.
  48 //
  49 // Upon creation, the class computes and caches information needed to
  50 // compute v, and then uses it to efficiently compute the product
  51 // every time InnerProductComputer::Compute is called.
  52 //
  53 // See sparse_normal_cholesky_solver.cc for example usage.
  54 //
  55 // Note that the result matrix is a block upper or lower-triangular
  56 // matrix, i.e., it will contain entries in the upper or lower
  57 // triangular part of the matrix corresponding to the block that occur
  58 // along its diagonal.
  59 //
  60 // This is not a problem as sparse linear algebra libraries can ignore
  61 // these entries with ease and the space used is minimal/linear in the
  62 // size of the matrices.
  63 class InnerProductComputer {
  64  public:
  65   // Factory
  66   //
  67   // m is the input matrix
  68   //
  69   // Since m' * m is a symmetric matrix, we only compute half of the
  70   // matrix and the value of storage_type which must be
  71   // UPPER_TRIANGULAR or LOWER_TRIANGULAR determines which half is
  72   // computed.
  73   //
  74   // The user must ensure that the matrix m is valid for the life time
  75   // of this object.
  76   static InnerProductComputer* Create(
  77       const BlockSparseMatrix& m,
  78       CompressedRowSparseMatrix::StorageType storage_type);
  79
  80   // This factory method allows the user control over range of row
  81   // blocks of m that should be used to compute the inner product.
  82   //
  83   // a = m(start_row_block : end_row_block, :);
  84   // result = a' * a;
  85   static InnerProductComputer* Create(
  86       const BlockSparseMatrix& m,
  87       int start_row_block,
  88       int end_row_block,
  89       CompressedRowSparseMatrix::StorageType storage_type);
  90
  91   // Update result_ to be numerically equal to m' * m.
  92   void Compute();
  93
  94   // Accessors for the result containing the inner product.
  95   //
  96   // Compute must be called before accessing this result for
  97   // the first time.
  98   const CompressedRowSparseMatrix& result() const { return *result_; }
  99   CompressedRowSparseMatrix* mutable_result() const { return result_.get(); }
 100
 101  private:
 102   // A ProductTerm is a term in the block inner product of a matrix
 103   // with itself.
 104   struct ProductTerm {
 105     ProductTerm(const int row, const int col, const int index)
 106         : row(row), col(col), index(index) {}
 107
 108     bool operator<(const ProductTerm& right) const {
 109       if (row == right.row) {
 110         if (col == right.col) {
 111           return index < right.index;
 112         }
 113         return col < right.col;
 114       }
 115       return row < right.row;
 116     }
 117
 118     int row;
 119     int col;
 120     int index;
 121   };
 122
 123   InnerProductComputer(const BlockSparseMatrix& m,
 124                        int start_row_block,
 125                        int end_row_block);
 126
 127   void Init(CompressedRowSparseMatrix::StorageType storage_type);
 128
 129   CompressedRowSparseMatrix* CreateResultMatrix(
 130       const CompressedRowSparseMatrix::StorageType storage_type,
 131       int num_nonzeros);
 132
 133   int ComputeNonzeros(const std::vector<ProductTerm>& product_terms,
 134                       std::vector<int>* row_block_nnz);
 135
 136   void ComputeOffsetsAndCreateResultMatrix(
 137       const CompressedRowSparseMatrix::StorageType storage_type,
 138       const std::vector<ProductTerm>& product_terms);
 139
 140   const BlockSparseMatrix& m_;
 141   const int start_row_block_;
 142   const int end_row_block_;
 143   scoped_ptr<CompressedRowSparseMatrix> result_;
 144
 145   // For each term in the inner product, result_offsets_ contains the
 146   // location in the values array of the result_ matrix where it
 147   // should be stored.
 148   //
 149   // This is the principal look up table that allows this class to
 150   // compute the inner product fast.
 151   std::vector<int> result_offsets_;
 152 };
 153
 154 }  // namespace internal
 155 }  // namespace ceres
 156
 157 #endif  // CERES_INTERNAL_INNER_PRODUCT_COMPUTER_H_