Add InnerProductComputer
Add a class that given a block sparse matrix m will compute
the product m'*m efficiently.
This code is refactoring and cleanup of the code in
CompressedRowSparseMatrix devoted to computing the inner product.
In that class, the code is mistakenly said to be computing
the outer product. It is also devoted to computing the inner
product of a CompressedRowSparseMatrix with itself.
This code works with BlockSparseMatrix objects instead, which
are simpler to deal with as they are better structured to handle
block sparse matrices.
Change-Id: I920fee1a396bb0fcae9e6f7e46a308c7391d21aa
diff --git a/internal/ceres/CMakeLists.txt b/internal/ceres/CMakeLists.txt
index 8a56ca0..052e7a0 100644
--- a/internal/ceres/CMakeLists.txt
+++ b/internal/ceres/CMakeLists.txt
@@ -70,6 +70,7 @@
gradient_problem.cc
gradient_problem_solver.cc
implicit_schur_complement.cc
+ inner_product_computer.cc
is_close.cc
iterative_schur_complement_solver.cc
levenberg_marquardt_strategy.cc
@@ -316,6 +317,7 @@
ceres_test(graph_algorithms)
ceres_test(householder_vector)
ceres_test(implicit_schur_complement)
+ ceres_test(inner_product_computer)
ceres_test(invert_psd_matrix)
ceres_test(is_close)
ceres_test(iterative_schur_complement_solver)
diff --git a/internal/ceres/inner_product_computer.cc b/internal/ceres/inner_product_computer.cc
new file mode 100644
index 0000000..052bd5a
--- /dev/null
+++ b/internal/ceres/inner_product_computer.cc
@@ -0,0 +1,373 @@
+// Ceres Solver - A fast non-linear least squares minimizer
+// Copyright 2017 Google Inc. All rights reserved.
+// http://ceres-solver.org/
+//
+// Redistribution and use in source and binary forms, with or without
+// modification, are permitted provided that the following conditions are met:
+//
+// * Redistributions of source code must retain the above copyright notice,
+// this list of conditions and the following disclaimer.
+// * Redistributions in binary form must reproduce the above copyright notice,
+// this list of conditions and the following disclaimer in the documentation
+// and/or other materials provided with the distribution.
+// * Neither the name of Google Inc. nor the names of its contributors may be
+// used to endorse or promote products derived from this software without
+// specific prior written permission.
+//
+// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
+// AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
+// ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
+// LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
+// CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
+// SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
+// INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
+// CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
+// ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
+// POSSIBILITY OF SUCH DAMAGE.
+//
+// Author: sameeragarwal@google.com (Sameer Agarwal)
+
+#include "ceres/inner_product_computer.h"
+
+#include <algorithm>
+
+namespace ceres {
+namespace internal {
+namespace {
+
+// Compute the product (in MATLAB notation)
+//
+// c(0:a_cols, 0:b_cols) = a' * b
+//
+// Where:
+//
+// a is ab_rows x a_cols
+// b is ab_rows x b_cols
+// c is a_cos x c_col_stride
+//
+// a, b and c are row-major matrices.
+//
+// Performance note:
+// ----------------
+//
+// Technically this function is a repeat of a similarly named function
+// in small_blas.h but its performance is considerably better than
+// that of the version there due to the way it accesses memory.
+//
+// TODO(sameeragarwal): Measure and tune the performance of
+// small_blas.h based on the insights gained here.
+EIGEN_STRONG_INLINE void MatrixTransposeMatrixMultiply(const int ab_rows,
+ const double* a,
+ const int a_cols,
+ const double* b,
+ const int b_cols,
+ double* c,
+ int c_cols) {
+ // Compute c as the sum of ab_rows, rank 1 outer products of the
+ // corresponding rows of a and b.
+ for (int r = 0; r < ab_rows; ++r) {
+ double* c_r = c;
+ for (int i1 = 0; i1 < a_cols; ++i1) {
+ const double a_v = a[i1];
+ for (int i2 = 0; i2 < b_cols; ++i2) {
+ c_r[i2] += a_v * b[i2];
+ }
+ c_r += c_cols;
+ }
+ a += a_cols;
+ b += b_cols;
+ }
+}
+
+} // namespace
+
+// Create the CompressedRowSparseMatrix matrix that will contain the
+// inner product.
+//
+// storage_type controls whether the result matrix contains the upper
+// or the lower triangular part of the product.
+//
+// num_nonzeros is the number of non-zeros in the result matrix.
+CompressedRowSparseMatrix* InnerProductComputer::CreateResultMatrix(
+ const CompressedRowSparseMatrix::StorageType storage_type,
+ const int num_nonzeros) {
+ CompressedRowSparseMatrix* matrix =
+ new CompressedRowSparseMatrix(m_.num_cols(), m_.num_cols(), num_nonzeros);
+ matrix->set_storage_type(storage_type);
+
+ const CompressedRowBlockStructure* bs = m_.block_structure();
+ const std::vector<Block>& blocks = bs->cols;
+ matrix->mutable_row_blocks()->resize(blocks.size());
+ matrix->mutable_col_blocks()->resize(blocks.size());
+ for (int i = 0; i < blocks.size(); ++i) {
+ (*(matrix->mutable_row_blocks()))[i] = blocks[i].size;
+ (*(matrix->mutable_col_blocks()))[i] = blocks[i].size;
+ }
+
+ return matrix;
+}
+
+// Given the set of product terms in the inner product, return the
+// total number of non-zeros in the result and for each row block of
+// the result matrix, compute the number of non-zeros in any one row
+// of the row block.
+int InnerProductComputer::ComputeNonzeros(
+ const std::vector<InnerProductComputer::ProductTerm>& product_terms,
+ std::vector<int>* row_nnz) {
+ const CompressedRowBlockStructure* bs = m_.block_structure();
+ const std::vector<Block>& blocks = bs->cols;
+
+ row_nnz->resize(blocks.size());
+ std::fill(row_nnz->begin(), row_nnz->end(), 0);
+
+ // First product term.
+ (*row_nnz)[product_terms[0].row] = blocks[product_terms[0].col].size;
+ int num_nonzeros =
+ blocks[product_terms[0].row].size * blocks[product_terms[0].col].size;
+
+ // Remaining product terms.
+ for (int i = 1; i < product_terms.size(); ++i) {
+ const ProductTerm& previous = product_terms[i - 1];
+ const ProductTerm& current = product_terms[i];
+
+ // Each (row, col) block counts only once.
+ // This check depends on product sorted on (row, col).
+ if (current.row != previous.row || current.col != previous.col) {
+ (*row_nnz)[current.row] += blocks[current.col].size;
+ num_nonzeros += blocks[current.row].size * blocks[current.col].size;
+ }
+ }
+
+ return num_nonzeros;
+}
+
+InnerProductComputer::InnerProductComputer(const BlockSparseMatrix& m,
+ const int start_row_block,
+ const int end_row_block)
+ : m_(m), start_row_block_(start_row_block), end_row_block_(end_row_block) {}
+
+// Compute the sparsity structure of the product m.transpose() * m
+// and create a CompressedRowSparseMatrix corresponding to it.
+//
+// Also compute the "program" vector, which for every term in the
+// block outer product provides the information for the entry in the
+// values array of the result matrix where it should be accumulated.
+//
+// Since the entries of the program are the same for rows with the
+// same sparsity structure, the program only stores the result for one
+// row per row block. The Compute function reuses this information for
+// each row in the row block.
+//
+// product_storage_type controls the form of the output matrix. It
+// can be LOWER_TRIANGULAR or UPPER_TRIANGULAR.
+InnerProductComputer* InnerProductComputer::Create(
+ const BlockSparseMatrix& m,
+ CompressedRowSparseMatrix::StorageType product_storage_type) {
+ return InnerProductComputer::Create(
+ m, 0, m.block_structure()->rows.size(), product_storage_type);
+}
+
+InnerProductComputer* InnerProductComputer::Create(
+ const BlockSparseMatrix& m,
+ const int start_row_block,
+ const int end_row_block,
+ CompressedRowSparseMatrix::StorageType product_storage_type) {
+ CHECK(product_storage_type == CompressedRowSparseMatrix::LOWER_TRIANGULAR ||
+ product_storage_type == CompressedRowSparseMatrix::UPPER_TRIANGULAR);
+ CHECK_GT(m.num_nonzeros(), 0)
+ << "Congratulations, you found a bug in Ceres. Please report it.";
+ InnerProductComputer* inner_product_computer =
+ new InnerProductComputer(m, start_row_block, end_row_block);
+ inner_product_computer->Init(product_storage_type);
+ return inner_product_computer;
+}
+
+void InnerProductComputer::Init(
+ const CompressedRowSparseMatrix::StorageType product_storage_type) {
+ std::vector<InnerProductComputer::ProductTerm> product_terms;
+ const CompressedRowBlockStructure* bs = m_.block_structure();
+
+ // Give input matrix m in Block Sparse format
+ // (row_block, col_block)
+ // represent each block multiplication
+ // (row_block, col_block1)' X (row_block, col_block2)
+ // by its product term:
+ // (col_block1, col_block2, index)
+ for (int row_block = start_row_block_; row_block < end_row_block_;
+ ++row_block) {
+ const CompressedRow& row = bs->rows[row_block];
+ for (int c1 = 0; c1 < row.cells.size(); ++c1) {
+ const Cell& cell1 = row.cells[c1];
+ int c2_begin, c2_end;
+ if (product_storage_type == CompressedRowSparseMatrix::LOWER_TRIANGULAR) {
+ c2_begin = 0;
+ c2_end = c1 + 1;
+ } else {
+ c2_begin = c1;
+ c2_end = row.cells.size();
+ }
+
+ for (int c2 = c2_begin; c2 < c2_end; ++c2) {
+ const Cell& cell2 = row.cells[c2];
+ product_terms.push_back(InnerProductComputer::ProductTerm(
+ cell1.block_id, cell2.block_id, product_terms.size()));
+ }
+ }
+ }
+
+ std::sort(product_terms.begin(), product_terms.end());
+ ComputeOffsetsAndCreateResultMatrix(product_storage_type, product_terms);
+}
+
+void InnerProductComputer::ComputeOffsetsAndCreateResultMatrix(
+ const CompressedRowSparseMatrix::StorageType product_storage_type,
+ const std::vector<InnerProductComputer::ProductTerm>& product_terms) {
+ const int num_cols = m_.num_cols();
+ const std::vector<Block>& col_blocks = m_.block_structure()->cols;
+
+ std::vector<int> row_block_nnz;
+ const int num_nonzeros = ComputeNonzeros(product_terms, &row_block_nnz);
+
+ result_.reset(CreateResultMatrix(product_storage_type, num_nonzeros));
+
+ // Populate the row non-zero counts in the result matrix.
+ int* crsm_rows = result_->mutable_rows();
+ crsm_rows[0] = 0;
+ for (int i = 0; i < col_blocks.size(); ++i) {
+ for (int j = 0; j < col_blocks[i].size; ++j, ++crsm_rows) {
+ *(crsm_rows + 1) = *crsm_rows + row_block_nnz[i];
+ }
+ }
+
+ // The following macro FILL_CRSM_COL_BLOCK is key to understanding
+ // how this class works.
+ //
+ // It does two things.
+ //
+ // Sets the value for the current term in the result_offsets_ array
+ // and populates the cols array of the result matrix.
+ //
+ // row_block and col_block as the names imply, refer to the row and
+ // column blocks of the current term.
+ //
+ // row_nnz is the number of nonzeros in the result_matrix at the
+ // beginning of the first row of row_block.
+ //
+ // col_nnz is the number of nonzeros in the first row of the row
+ // block that occur before the current column block, i.e. this is
+ // sum of the sizes of all the column blocks in this row block that
+ // came before this column block.
+ //
+ // Given these two numbers and the total number of nonzeros in this
+ // row (nnz_in_row), we can now populate the cols array as follows:
+ //
+ // nnz + j * nnz_in_row is the beginning of the j^th row.
+ //
+ // nnz + j * nnz_in_row + col_nnz is the beginning of the column
+ // block in the j^th row.
+ //
+ // nnz + j * nnz_in_row + col_nnz + k is then the j^th row and the
+ // k^th column of the product block, whose value is
+ //
+ // col_blocks[col_block].position + k, which is the column number of
+ // the k^th column of the current column block.
+#define FILL_CRSM_COL_BLOCK \
+ const int row_block = current->row; \
+ const int col_block = current->col; \
+ const int nnz_in_row = row_block_nnz[row_block]; \
+ int* crsm_cols = result_->mutable_cols(); \
+ result_offsets_[current->index] = nnz + col_nnz; \
+ for (int j = 0; j < col_blocks[row_block].size; ++j) { \
+ for (int k = 0; k < col_blocks[col_block].size; ++k) { \
+ crsm_cols[nnz + j * nnz_in_row + col_nnz + k] = \
+ col_blocks[col_block].position + k; \
+ } \
+ }
+
+ result_offsets_.resize(product_terms.size());
+ int col_nnz = 0;
+ int nnz = 0;
+
+ // Process the first term.
+ const InnerProductComputer::ProductTerm* current = &product_terms[0];
+ FILL_CRSM_COL_BLOCK;
+
+ // Process the rest of the terms.
+ for (int i = 1; i < product_terms.size(); ++i) {
+ current = &product_terms[i];
+ const InnerProductComputer::ProductTerm* previous = &product_terms[i - 1];
+
+ // If the current term is the same as the previous term, then it
+ // stores its product at the same location as the previous term.
+ if (previous->row == current->row && previous->col == current->col) {
+ result_offsets_[current->index] = result_offsets_[previous->index];
+ continue;
+ }
+
+ if (previous->row == current->row) {
+ // if the current and previous terms are in the same row block,
+ // then they differ in the column block, in which case advance
+ // col_nnz by the column size of the prevous term.
+ col_nnz += col_blocks[previous->col].size;
+ } else {
+ // If we have moved to a new row-block , then col_nnz is zero,
+ // and nnz is set to the beginning of the row block.
+ col_nnz = 0;
+ nnz += row_block_nnz[previous->row] * col_blocks[previous->row].size;
+ }
+
+ FILL_CRSM_COL_BLOCK;
+ }
+}
+
+// Use the results_offsets_ array to numerically compute the product
+// m' * m and store it in result_.
+//
+// TODO(sameeragarwal): Multithreading support.
+void InnerProductComputer::Compute() {
+ const double* m_values = m_.values();
+ const CompressedRowBlockStructure* bs = m_.block_structure();
+
+ const CompressedRowSparseMatrix::StorageType storage_type =
+ result_->storage_type();
+ result_->SetZero();
+ double* values = result_->mutable_values();
+ const int* rows = result_->rows();
+ int cursor = 0;
+
+ // Iterate row blocks.
+ for (int r = start_row_block_; r < end_row_block_; ++r) {
+ const CompressedRow& m_row = bs->rows[r];
+ for (int c1 = 0; c1 < m_row.cells.size(); ++c1) {
+ const Cell& cell1 = m_row.cells[c1];
+ const int c1_size = bs->cols[cell1.block_id].size;
+ const int row_nnz = rows[bs->cols[cell1.block_id].position + 1] -
+ rows[bs->cols[cell1.block_id].position];
+
+ int c2_begin, c2_end;
+ if (storage_type == CompressedRowSparseMatrix::LOWER_TRIANGULAR) {
+ c2_begin = 0;
+ c2_end = c1 + 1;
+ } else {
+ c2_begin = c1;
+ c2_end = m_row.cells.size();
+ }
+
+ for (int c2 = c2_begin; c2 < c2_end; ++c2, ++cursor) {
+ const Cell& cell2 = m_row.cells[c2];
+ const int c2_size = bs->cols[cell2.block_id].size;
+ MatrixTransposeMatrixMultiply(m_row.block.size,
+ m_values + cell1.position, c1_size,
+ m_values + cell2.position, c2_size,
+ values + result_offsets_[cursor],
+ row_nnz);
+ }
+ }
+ }
+
+ CHECK_EQ(cursor, result_offsets_.size());
+}
+
+} // namespace internal
+} // namespace ceres
diff --git a/internal/ceres/inner_product_computer.h b/internal/ceres/inner_product_computer.h
new file mode 100644
index 0000000..d206707
--- /dev/null
+++ b/internal/ceres/inner_product_computer.h
@@ -0,0 +1,157 @@
+// Ceres Solver - A fast non-linear least squares minimizer
+// Copyright 2017 Google Inc. All rights reserved.
+// http://ceres-solver.org/
+//
+// Redistribution and use in source and binary forms, with or without
+// modification, are permitted provided that the following conditions are met:
+//
+// * Redistributions of source code must retain the above copyright notice,
+// this list of conditions and the following disclaimer.
+// * Redistributions in binary form must reproduce the above copyright notice,
+// this list of conditions and the following disclaimer in the documentation
+// and/or other materials provided with the distribution.
+// * Neither the name of Google Inc. nor the names of its contributors may be
+// used to endorse or promote products derived from this software without
+// specific prior written permission.
+//
+// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
+// AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
+// ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
+// LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
+// CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
+// SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
+// INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
+// CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
+// ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
+// POSSIBILITY OF SUCH DAMAGE.
+//
+// Author: sameeragarwal@google.com (Sameer Agarwal)
+
+#ifndef CERES_INTERNAL_INNER_PRODUCT_COMPUTER_H_
+#define CERES_INTERNAL_INNER_PRODUCT_COMPUTER_H_
+
+#include <vector>
+
+#include "ceres/block_sparse_matrix.h"
+#include "ceres/compressed_row_sparse_matrix.h"
+#include "ceres/internal/scoped_ptr.h"
+
+namespace ceres {
+namespace internal {
+
+// This class is used to repeatedly compute the inner product
+//
+// result = m' * m
+//
+// where the sparsity structure of m remains constant across calls.
+//
+// Upon creation, the class computes and caches information needed to
+// compute v, and then uses it to efficiently compute the product
+// every time InnerProductComputer::Compute is called.
+//
+// See sparse_normal_cholesky_solver.cc for example usage.
+//
+// Note that the result matrix is a block upper or lower-triangular
+// matrix, i.e., it will contain entries in the upper or lower
+// triangular part of the matrix corresponding to the block that occur
+// along its diagonal.
+//
+// This is not a problem as sparse linear algebra libraries can ignore
+// these entries with ease and the space used is minimal/linear in the
+// size of the matrices.
+class InnerProductComputer {
+ public:
+ // Factory
+ //
+ // m is the input matrix
+ //
+ // Since m' * m is a symmetric matrix, we only compute half of the
+ // matrix and the value of storage_type which must be
+ // UPPER_TRIANGULAR or LOWER_TRIANGULAR determines which half is
+ // computed.
+ //
+ // The user must ensure that the matrix m is valid for the life time
+ // of this object.
+ static InnerProductComputer* Create(
+ const BlockSparseMatrix& m,
+ CompressedRowSparseMatrix::StorageType storage_type);
+
+ // This factory method allows the user control over range of row
+ // blocks of m that should be used to compute the inner product.
+ //
+ // a = m(start_row_block : end_row_block, :);
+ // result = a' * a;
+ static InnerProductComputer* Create(
+ const BlockSparseMatrix& m,
+ int start_row_block,
+ int end_row_block,
+ CompressedRowSparseMatrix::StorageType storage_type);
+
+ // Update result_ to be numerically equal to m' * m.
+ void Compute();
+
+ // Accessors for the result containing the inner product.
+ //
+ // Compute must be called before accessing this result for
+ // the first time.
+ const CompressedRowSparseMatrix& result() const { return *result_; }
+ CompressedRowSparseMatrix* mutable_result() const { return result_.get(); }
+
+ private:
+ // A ProductTerm is a term in the block inner product of a matrix
+ // with itself.
+ struct ProductTerm {
+ ProductTerm(const int row, const int col, const int index)
+ : row(row), col(col), index(index) {}
+
+ bool operator<(const ProductTerm& right) const {
+ if (row == right.row) {
+ if (col == right.col) {
+ return index < right.index;
+ }
+ return col < right.col;
+ }
+ return row < right.row;
+ }
+
+ int row;
+ int col;
+ int index;
+ };
+
+ InnerProductComputer(const BlockSparseMatrix& m,
+ int start_row_block,
+ int end_row_block);
+
+ void Init(CompressedRowSparseMatrix::StorageType storage_type);
+
+ CompressedRowSparseMatrix* CreateResultMatrix(
+ const CompressedRowSparseMatrix::StorageType storage_type,
+ int num_nonzeros);
+
+ int ComputeNonzeros(const std::vector<ProductTerm>& product_terms,
+ std::vector<int>* row_block_nnz);
+
+ void ComputeOffsetsAndCreateResultMatrix(
+ const CompressedRowSparseMatrix::StorageType storage_type,
+ const std::vector<ProductTerm>& product_terms);
+
+ const BlockSparseMatrix& m_;
+ const int start_row_block_;
+ const int end_row_block_;
+ scoped_ptr<CompressedRowSparseMatrix> result_;
+
+ // For each term in the inner product, result_offsets_ contains the
+ // location in the values array of the result_ matrix where it
+ // should be stored.
+ //
+ // This is the principal look up table that allows this class to
+ // compute the inner product fast.
+ std::vector<int> result_offsets_;
+};
+
+} // namespace internal
+} // namespace ceres
+
+#endif // CERES_INTERNAL_INNER_PRODUCT_COMPUTER_H_
diff --git a/internal/ceres/inner_product_computer_test.cc b/internal/ceres/inner_product_computer_test.cc
new file mode 100644
index 0000000..a4c81bc
--- /dev/null
+++ b/internal/ceres/inner_product_computer_test.cc
@@ -0,0 +1,236 @@
+// Ceres Solver - A fast non-linear least squares minimizer
+// Copyright 2017 Google Inc. All rights reserved.
+// http://ceres-solver.org/
+//
+// Redistribution and use in source and binary forms, with or without
+// modification, are permitted provided that the following conditions are met:
+//
+// * Redistributions of source code must retain the above copyright notice,
+// this list of conditions and the following disclaimer.
+// * Redistributions in binary form must reproduce the above copyright notice,
+// this list of conditions and the following disclaimer in the documentation
+// and/or other materials provided with the distribution.
+// * Neither the name of Google Inc. nor the names of its contributors may be
+// used to endorse or promote products derived from this software without
+// specific prior written permission.
+//
+// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
+// AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
+// ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
+// LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
+// CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
+// SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
+// INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
+// CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
+// ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
+// POSSIBILITY OF SUCH DAMAGE.
+//
+// Author: sameeragarwal@google.com (Sameer Agarwal)
+
+#include "ceres/inner_product_computer.h"
+
+#include <numeric>
+#include "ceres/block_sparse_matrix.h"
+#include "ceres/internal/eigen.h"
+#include "ceres/internal/scoped_ptr.h"
+#include "ceres/random.h"
+#include "ceres/triplet_sparse_matrix.h"
+#include "glog/logging.h"
+#include "gtest/gtest.h"
+
+#include "Eigen/SparseCore"
+
+namespace ceres {
+namespace internal {
+
+template <enum Eigen::UpLoType T>
+void CompareTriangularPartOfMatrices(const Matrix& expected,
+ const Matrix& actual) {
+ EXPECT_EQ(actual.rows(), actual.cols());
+ EXPECT_EQ(expected.rows(), expected.cols());
+ EXPECT_EQ(actual.rows(), expected.rows());
+
+ const Matrix expected_t = expected.triangularView<T>();
+ const Matrix actual_t = actual.triangularView<T>();
+
+ // TODO(sameeragarwal): Foo
+ CHECK_LE((expected_t - actual_t).norm() / actual_t.norm(),
+ 100 * std::numeric_limits<double>::epsilon())
+ << "expected: \n"
+ << expected_t << "\nactual: \n"
+ << actual_t;
+}
+
+#define COMPUTE_AND_COMPARE \
+ { \
+ inner_product_computer->Compute(); \
+ CompressedRowSparseMatrix* actual_product_crsm = \
+ inner_product_computer->mutable_result(); \
+ Matrix actual_inner_product = \
+ Eigen::MappedSparseMatrix<double, Eigen::ColMajor>( \
+ actual_product_crsm->num_rows(), \
+ actual_product_crsm->num_rows(), \
+ actual_product_crsm->num_nonzeros(), \
+ actual_product_crsm->mutable_rows(), \
+ actual_product_crsm->mutable_cols(), \
+ actual_product_crsm->mutable_values()); \
+ (actual_product_crsm->storage_type() == \
+ CompressedRowSparseMatrix::LOWER_TRIANGULAR) \
+ ? CompareTriangularPartOfMatrices<Eigen::Upper>( \
+ expected_inner_product, actual_inner_product) \
+ : CompareTriangularPartOfMatrices<Eigen::Lower>( \
+ expected_inner_product, actual_inner_product); \
+ }
+
+
+TEST(InnerProductComputer, NormalOperation) {
+ // "Randomly generated seed."
+ SetRandomState(29823);
+ const int kMaxNumRowBlocks = 10;
+ const int kMaxNumColBlocks = 10;
+ const int kNumTrials = 10;
+
+ // Create a random matrix, compute its outer product using Eigen and
+ // ComputeOuterProduct. Convert both matrices to dense matrices and
+ // compare their upper triangular parts.
+ for (int num_row_blocks = 1; num_row_blocks < kMaxNumRowBlocks;
+ ++num_row_blocks) {
+ for (int num_col_blocks = 1; num_col_blocks < kMaxNumColBlocks;
+ ++num_col_blocks) {
+ for (int trial = 0; trial < kNumTrials; ++trial) {
+ BlockSparseMatrix::RandomMatrixOptions options;
+ options.num_row_blocks = num_row_blocks;
+ options.num_col_blocks = num_col_blocks;
+ options.min_row_block_size = 1;
+ options.max_row_block_size = 5;
+ options.min_col_block_size = 1;
+ options.max_col_block_size = 10;
+ options.block_density = std::max(0.1, RandDouble());
+
+ VLOG(2) << "num row blocks: " << options.num_row_blocks;
+ VLOG(2) << "num col blocks: " << options.num_col_blocks;
+ VLOG(2) << "min row block size: " << options.min_row_block_size;
+ VLOG(2) << "max row block size: " << options.max_row_block_size;
+ VLOG(2) << "min col block size: " << options.min_col_block_size;
+ VLOG(2) << "max col block size: " << options.max_col_block_size;
+ VLOG(2) << "block density: " << options.block_density;
+
+ scoped_ptr<BlockSparseMatrix> random_matrix(
+ BlockSparseMatrix::CreateRandomMatrix(options));
+
+ TripletSparseMatrix tsm(random_matrix->num_rows(),
+ random_matrix->num_cols(),
+ random_matrix->num_nonzeros());
+ random_matrix->ToTripletSparseMatrix(&tsm);
+ std::vector<Eigen::Triplet<double> > triplets;
+ for (int i = 0; i < tsm.num_nonzeros(); ++i) {
+ triplets.push_back(Eigen::Triplet<double>(
+ tsm.rows()[i], tsm.cols()[i], tsm.values()[i]));
+ }
+ Eigen::SparseMatrix<double> eigen_random_matrix(
+ random_matrix->num_rows(), random_matrix->num_cols());
+ eigen_random_matrix.setFromTriplets(triplets.begin(), triplets.end());
+ Matrix expected_inner_product =
+ eigen_random_matrix.transpose() * eigen_random_matrix;
+
+ scoped_ptr<InnerProductComputer> inner_product_computer;
+
+ inner_product_computer.reset(InnerProductComputer::Create(
+ *random_matrix, CompressedRowSparseMatrix::LOWER_TRIANGULAR));
+ COMPUTE_AND_COMPARE;
+ inner_product_computer.reset(InnerProductComputer::Create(
+ *random_matrix, CompressedRowSparseMatrix::UPPER_TRIANGULAR));
+ COMPUTE_AND_COMPARE;
+
+ }
+ }
+ }
+}
+
+
+TEST(InnerProductComputer, SubMatrix) {
+ // "Randomly generated seed."
+ SetRandomState(29823);
+ const int kNumRowBlocks = 10;
+ const int kNumColBlocks = 20;
+ const int kNumTrials = 5;
+
+ // Create a random matrix, compute its outer product using Eigen and
+ // ComputeInnerProductComputer. Convert both matrices to dense matrices and
+ // compare their upper triangular parts.
+ for (int trial = 0; trial < kNumTrials; ++trial) {
+ BlockSparseMatrix::RandomMatrixOptions options;
+ options.num_row_blocks = kNumRowBlocks;
+ options.num_col_blocks = kNumColBlocks;
+ options.min_row_block_size = 1;
+ options.max_row_block_size = 5;
+ options.min_col_block_size = 1;
+ options.max_col_block_size = 10;
+ options.block_density = std::max(0.1, RandDouble());
+
+ VLOG(2) << "num row blocks: " << options.num_row_blocks;
+ VLOG(2) << "num col blocks: " << options.num_col_blocks;
+ VLOG(2) << "min row block size: " << options.min_row_block_size;
+ VLOG(2) << "max row block size: " << options.max_row_block_size;
+ VLOG(2) << "min col block size: " << options.min_col_block_size;
+ VLOG(2) << "max col block size: " << options.max_col_block_size;
+ VLOG(2) << "block density: " << options.block_density;
+
+ scoped_ptr<BlockSparseMatrix> random_matrix(
+ BlockSparseMatrix::CreateRandomMatrix(options));
+
+ const std::vector<CompressedRow>& row_blocks =
+ random_matrix->block_structure()->rows;
+ const int num_row_blocks = row_blocks.size();
+
+ for (int start_row_block = 0; start_row_block < num_row_blocks - 1;
+ ++start_row_block) {
+ for (int end_row_block = start_row_block + 1;
+ end_row_block < num_row_blocks;
+ ++end_row_block) {
+ const int start_row = row_blocks[start_row_block].block.position;
+ const int end_row = row_blocks[end_row_block].block.position;
+
+ TripletSparseMatrix tsm(random_matrix->num_rows(),
+ random_matrix->num_cols(),
+ random_matrix->num_nonzeros());
+ random_matrix->ToTripletSparseMatrix(&tsm);
+ std::vector<Eigen::Triplet<double> > triplets;
+ for (int i = 0; i < tsm.num_nonzeros(); ++i) {
+ if (tsm.rows()[i] >= start_row && tsm.rows()[i] < end_row) {
+ triplets.push_back(Eigen::Triplet<double>(
+ tsm.rows()[i], tsm.cols()[i], tsm.values()[i]));
+ }
+ }
+
+ Eigen::SparseMatrix<double> eigen_random_matrix(
+ random_matrix->num_rows(), random_matrix->num_cols());
+ eigen_random_matrix.setFromTriplets(triplets.begin(), triplets.end());
+
+ Matrix expected_inner_product =
+ eigen_random_matrix.transpose() * eigen_random_matrix;
+
+ scoped_ptr<InnerProductComputer> inner_product_computer;
+ inner_product_computer.reset(InnerProductComputer::Create(
+ *random_matrix,
+ start_row_block,
+ end_row_block,
+ CompressedRowSparseMatrix::LOWER_TRIANGULAR));
+ COMPUTE_AND_COMPARE;
+ inner_product_computer.reset(InnerProductComputer::Create(
+ *random_matrix,
+ start_row_block,
+ end_row_block,
+ CompressedRowSparseMatrix::UPPER_TRIANGULAR));
+ COMPUTE_AND_COMPARE;
+
+ }
+ }
+ }
+}
+
+#undef COMPUTE_AND_COMPARE
+
+} // namespace internal
+} // namespace ceres
diff --git a/jni/Android.mk b/jni/Android.mk
index 20b484d..551cd45 100644
--- a/jni/Android.mk
+++ b/jni/Android.mk
@@ -154,6 +154,7 @@
$(CERES_SRC_PATH)/gradient_problem_solver.cc \
$(CERES_SRC_PATH)/is_close.cc \
$(CERES_SRC_PATH)/implicit_schur_complement.cc \
+ $(CERES_SRC_PATH)/inner_product_computer.cc \
$(CERES_SRC_PATH)/iterative_schur_complement_solver.cc \
$(CERES_SRC_PATH)/lapack.cc \
$(CERES_SRC_PATH)/levenberg_marquardt_strategy.cc \
