| // Ceres Solver - A fast non-linear least squares minimizer |
| // Copyright 2023 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> |
| #include <memory> |
| |
| #include "absl/log/check.h" |
| #include "ceres/small_blas.h" |
| |
| namespace ceres::internal { |
| |
| // 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. |
| std::unique_ptr<CompressedRowSparseMatrix> |
| InnerProductComputer::CreateResultMatrix( |
| const CompressedRowSparseMatrix::StorageType storage_type, |
| const int num_nonzeros) { |
| auto matrix = std::make_unique<CompressedRowSparseMatrix>( |
| m_.num_cols(), m_.num_cols(), num_nonzeros); |
| matrix->set_storage_type(storage_type); |
| const CompressedRowBlockStructure* bs = m_.block_structure(); |
| *matrix->mutable_row_blocks() = bs->cols; |
| *matrix->mutable_col_blocks() = bs->cols; |
| 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); |
| |
| if (product_terms.empty()) { |
| return 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. |
| std::unique_ptr<InnerProductComputer> InnerProductComputer::Create( |
| const BlockSparseMatrix& m, |
| CompressedRowSparseMatrix::StorageType product_storage_type) { |
| return InnerProductComputer::Create( |
| m, 0, m.block_structure()->rows.size(), product_storage_type); |
| } |
| |
| std::unique_ptr<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::StorageType::LOWER_TRIANGULAR || |
| product_storage_type == |
| CompressedRowSparseMatrix::StorageType::UPPER_TRIANGULAR); |
| CHECK_GT(m.num_nonzeros(), 0) |
| << "Congratulations, you found a bug in Ceres. Please report it."; |
| std::unique_ptr<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::StorageType::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.emplace_back( |
| 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 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_ = 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]; |
| } |
| } |
| result_offsets_.resize(product_terms.size()); |
| if (num_nonzeros == 0) { |
| return; |
| } |
| |
| // 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; \ |
| } \ |
| } |
| |
| int col_nnz = 0; |
| int nnz = 0; |
| |
| // Process the first term. |
| const InnerProductComputer::ProductTerm* current = product_terms.data(); |
| 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 previous 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::StorageType::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; |
| // clang-format off |
| MatrixTransposeMatrixMultiply<Eigen::Dynamic, Eigen::Dynamic, |
| Eigen::Dynamic, Eigen::Dynamic, 1>( |
| m_values + cell1.position, |
| m_row.block.size, c1_size, |
| m_values + cell2.position, |
| m_row.block.size, c2_size, |
| values + result_offsets_[cursor], |
| 0, 0, c1_size, row_nnz); |
| // clang-format on |
| } |
| } |
| } |
| |
| CHECK_EQ(cursor, result_offsets_.size()); |
| } |
| |
| } // namespace ceres::internal |