AddBlockStructureTranspose to BlockSparseMatrix Add structure of transposed matrix to BlockSparseMatrix Number of non-zero values per row block and cumulative non-zero values count are maintained for transposed structure Change-Id: Icf38bb7a734ca695c788579eece1c92d36d78e54
diff --git a/internal/ceres/block_sparse_matrix.cc b/internal/ceres/block_sparse_matrix.cc index 0dd3b20..0e05693 100644 --- a/internal/ceres/block_sparse_matrix.cc +++ b/internal/ceres/block_sparse_matrix.cc
@@ -46,6 +46,19 @@ namespace ceres::internal { +namespace { +void ComputeCumulativeNumberOfNonZeros(std::vector<CompressedList>& rows) { + if (!rows.size()) { + return; + } + rows[0].cumulative_nnz = rows[0].nnz; + for (int c = 1; c < rows.size(); ++c) { + const int curr_nnz = rows[c].nnz; + rows[c].cumulative_nnz = curr_nnz + rows[c - 1].cumulative_nnz; + } +} +} // namespace + BlockSparseMatrix::BlockSparseMatrix( CompressedRowBlockStructure* block_structure) : num_rows_(0), @@ -83,6 +96,12 @@ CHECK(values_ != nullptr); } +void BlockSparseMatrix::AddTransposeBlockStructure() { + if (transpose_block_structure_ == nullptr) { + transpose_block_structure_ = CreateTranspose(*block_structure_); + } +} + void BlockSparseMatrix::SetZero() { std::fill(values_.get(), values_.get() + num_nonzeros_, 0.0); } @@ -133,7 +152,6 @@ double* y) const { CHECK(x != nullptr); CHECK(y != nullptr); - for (int i = 0; i < block_structure_->rows.size(); ++i) { int row_block_pos = block_structure_->rows[i].block.position; int row_block_size = block_structure_->rows[i].block.size; @@ -267,6 +285,13 @@ return block_structure_.get(); } +// Return a pointer to the block structure of matrix transpose. We continue to +// hold ownership of the object though. +const CompressedRowBlockStructure* +BlockSparseMatrix::transpose_block_structure() const { + return transpose_block_structure_.get(); +} + void BlockSparseMatrix::ToTextFile(FILE* file) const { CHECK(file != nullptr); for (int i = 0; i < block_structure_->rows.size(); ++i) { @@ -337,16 +362,28 @@ for (int i = 0; i < m_bs->rows.size(); ++i) { const CompressedRow& m_row = m_bs->rows[i]; - CompressedRow& row = block_structure_->rows[old_num_row_blocks + i]; + const int row_block_id = old_num_row_blocks + i; + CompressedRow& row = block_structure_->rows[row_block_id]; row.block.size = m_row.block.size; row.block.position = num_rows_; num_rows_ += m_row.block.size; row.cells.resize(m_row.cells.size()); + if (transpose_block_structure_) { + transpose_block_structure_->cols.emplace_back(row.block); + } for (int c = 0; c < m_row.cells.size(); ++c) { const int block_id = m_row.cells[c].block_id; row.cells[c].block_id = block_id; row.cells[c].position = num_nonzeros_; - num_nonzeros_ += m_row.block.size * m_bs->cols[block_id].size; + + const int cell_nnz = m_row.block.size * m_bs->cols[block_id].size; + if (transpose_block_structure_) { + transpose_block_structure_->rows[block_id].cells.emplace_back( + row_block_id, num_nonzeros_); + transpose_block_structure_->rows[block_id].nnz += cell_nnz; + } + + num_nonzeros_ += cell_nnz; } } @@ -360,10 +397,16 @@ std::copy(m.values(), m.values() + m.num_nonzeros(), values_.get() + old_num_nonzeros); + + if (transpose_block_structure_ == nullptr) { + return; + } + ComputeCumulativeNumberOfNonZeros(transpose_block_structure_->rows); } void BlockSparseMatrix::DeleteRowBlocks(const int delta_row_blocks) { const int num_row_blocks = block_structure_->rows.size(); + const int new_num_row_blocks = num_row_blocks - delta_row_blocks; int delta_num_nonzeros = 0; int delta_num_rows = 0; const std::vector<Block>& column_blocks = block_structure_->cols; @@ -373,11 +416,34 @@ for (int c = 0; c < row.cells.size(); ++c) { const Cell& cell = row.cells[c]; delta_num_nonzeros += row.block.size * column_blocks[cell.block_id].size; + + if (transpose_block_structure_) { + auto& col_cells = transpose_block_structure_->rows[cell.block_id].cells; + while (!col_cells.empty() && + col_cells.back().block_id >= new_num_row_blocks) { + const int del_block_id = col_cells.back().block_id; + const int del_block_rows = + block_structure_->rows[del_block_id].block.size; + const int del_block_cols = column_blocks[cell.block_id].size; + const int del_cell_nnz = del_block_rows * del_block_cols; + transpose_block_structure_->rows[cell.block_id].nnz -= del_cell_nnz; + col_cells.pop_back(); + } + } } } num_nonzeros_ -= delta_num_nonzeros; num_rows_ -= delta_num_rows; - block_structure_->rows.resize(num_row_blocks - delta_row_blocks); + block_structure_->rows.resize(new_num_row_blocks); + + if (transpose_block_structure_ == nullptr) { + return; + } + for (int i = 0; i < delta_row_blocks; ++i) { + transpose_block_structure_->cols.pop_back(); + } + + ComputeCumulativeNumberOfNonZeros(transpose_block_structure_->rows); } std::unique_ptr<BlockSparseMatrix> BlockSparseMatrix::CreateRandomMatrix( @@ -449,4 +515,30 @@ return matrix; } +std::unique_ptr<CompressedRowBlockStructure> CreateTranspose( + const CompressedRowBlockStructure& bs) { + auto transpose = std::make_unique<CompressedRowBlockStructure>(); + + transpose->rows.resize(bs.cols.size()); + for (int i = 0; i < bs.cols.size(); ++i) { + transpose->rows[i].block = bs.cols[i]; + transpose->rows[i].nnz = 0; + } + + transpose->cols.resize(bs.rows.size()); + for (int i = 0; i < bs.rows.size(); ++i) { + auto& row = bs.rows[i]; + transpose->cols[i] = row.block; + + const int nrows = row.block.size; + for (auto& cell : row.cells) { + transpose->rows[cell.block_id].cells.emplace_back(i, cell.position); + const int ncols = transpose->rows[cell.block_id].block.size; + transpose->rows[cell.block_id].nnz += nrows * ncols; + } + } + ComputeCumulativeNumberOfNonZeros(transpose->rows); + return transpose; +} + } // namespace ceres::internal
diff --git a/internal/ceres/block_sparse_matrix.h b/internal/ceres/block_sparse_matrix.h index 3a8dac9..b4abdc7 100644 --- a/internal/ceres/block_sparse_matrix.h +++ b/internal/ceres/block_sparse_matrix.h
@@ -84,6 +84,8 @@ void ToDenseMatrix(Matrix* dense_matrix) const final; void ToTextFile(FILE* file) const final; + void AddTransposeBlockStructure(); + // clang-format off int num_rows() const final { return num_rows_; } int num_cols() const final { return num_cols_; } @@ -94,6 +96,7 @@ void ToTripletSparseMatrix(TripletSparseMatrix* matrix) const; const CompressedRowBlockStructure* block_structure() const; + const CompressedRowBlockStructure* transpose_block_structure() const; // Append the contents of m to the bottom of this matrix. m must // have the same column blocks structure as this matrix. @@ -137,6 +140,7 @@ int max_num_nonzeros_; std::unique_ptr<double[]> values_; std::unique_ptr<CompressedRowBlockStructure> block_structure_; + std::unique_ptr<CompressedRowBlockStructure> transpose_block_structure_; }; // A number of algorithms like the SchurEliminator do not need @@ -164,6 +168,9 @@ const double* values_; }; +std::unique_ptr<CompressedRowBlockStructure> CreateTranspose( + const CompressedRowBlockStructure& bs); + } // namespace ceres::internal #include "ceres/internal/reenable_warnings.h"
diff --git a/internal/ceres/block_sparse_matrix_test.cc b/internal/ceres/block_sparse_matrix_test.cc index 3f5c78a..62bb0f4 100644 --- a/internal/ceres/block_sparse_matrix_test.cc +++ b/internal/ceres/block_sparse_matrix_test.cc
@@ -227,6 +227,85 @@ } } +TEST_F(BlockSparseMatrixTest, AppendDeleteRowsTransposedStructure) { + auto problem = CreateLinearLeastSquaresProblemFromId(2); + std::unique_ptr<BlockSparseMatrix> m( + down_cast<BlockSparseMatrix*>(problem->A.release())); + + A_->AddTransposeBlockStructure(); + auto block_structure = A_->block_structure(); + + // Several AppendRows and DeleteRowBlocks operations are applied to matrix, + // with regular and transpose block structures being compared after each + // operation. + // + // Non-negative values encode number of row blocks to remove + // -1 encodes appending matrix m + const int num_row_blocks_to_delete[] = {0, -1, 1, -1, 8, -1, 10}; + for (auto& t : num_row_blocks_to_delete) { + if (t == -1) { + A_->AppendRows(*m); + } else if (t > 0) { + CHECK_GE(block_structure->rows.size(), t); + A_->DeleteRowBlocks(t); + } + + auto block_structure = A_->block_structure(); + auto transpose_block_structure = A_->transpose_block_structure(); + ASSERT_NE(block_structure, nullptr); + ASSERT_NE(transpose_block_structure, nullptr); + + EXPECT_EQ(block_structure->rows.size(), + transpose_block_structure->cols.size()); + EXPECT_EQ(block_structure->cols.size(), + transpose_block_structure->rows.size()); + + std::vector<int> nnz_col(transpose_block_structure->rows.size()); + for (int i = 0; i < block_structure->cols.size(); ++i) { + EXPECT_EQ(block_structure->cols[i].position, + transpose_block_structure->rows[i].block.position); + const int col_size = transpose_block_structure->rows[i].block.size; + EXPECT_EQ(block_structure->cols[i].size, col_size); + + for (auto& col_cell : transpose_block_structure->rows[i].cells) { + int matches = 0; + const int row_block_id = col_cell.block_id; + nnz_col[i] += + col_size * transpose_block_structure->cols[row_block_id].size; + for (auto& row_cell : block_structure->rows[row_block_id].cells) { + if (row_cell.block_id != i) continue; + EXPECT_EQ(row_cell.position, col_cell.position); + ++matches; + } + EXPECT_EQ(matches, 1); + } + EXPECT_EQ(nnz_col[i], transpose_block_structure->rows[i].nnz); + if (i > 0) { + nnz_col[i] += nnz_col[i - 1]; + } + EXPECT_EQ(nnz_col[i], transpose_block_structure->rows[i].cumulative_nnz); + } + for (int i = 0; i < block_structure->rows.size(); ++i) { + EXPECT_EQ(block_structure->rows[i].block.position, + transpose_block_structure->cols[i].position); + EXPECT_EQ(block_structure->rows[i].block.size, + transpose_block_structure->cols[i].size); + + for (auto& row_cell : block_structure->rows[i].cells) { + int matches = 0; + const int col_block_id = row_cell.block_id; + for (auto& col_cell : + transpose_block_structure->rows[col_block_id].cells) { + if (col_cell.block_id != i) continue; + EXPECT_EQ(col_cell.position, row_cell.position); + ++matches; + } + EXPECT_EQ(matches, 1); + } + } + } +} + TEST_F(BlockSparseMatrixTest, AppendAndDeleteBlockDiagonalMatrix) { const std::vector<Block>& column_blocks = A_->block_structure()->cols; const int num_cols = @@ -368,5 +447,46 @@ } } +TEST(BlockSparseMatrix, CreateTranspose) { + constexpr int kNumtrials = 10; + BlockSparseMatrix::RandomMatrixOptions options; + options.num_col_blocks = 10; + options.min_col_block_size = 1; + options.max_col_block_size = 3; + + options.num_row_blocks = 20; + options.min_row_block_size = 1; + options.max_row_block_size = 4; + options.block_density = 0.25; + std::mt19937 prng; + + for (int trial = 0; trial < kNumtrials; ++trial) { + auto a = BlockSparseMatrix::CreateRandomMatrix(options, prng); + + auto ap_bs = std::make_unique<CompressedRowBlockStructure>(); + *ap_bs = *a->block_structure(); + BlockSparseMatrix ap(ap_bs.release()); + std::copy_n(a->values(), a->num_nonzeros(), ap.mutable_values()); + ap.AddTransposeBlockStructure(); + + Vector x = Vector::Random(a->num_cols()); + Vector y = Vector::Random(a->num_rows()); + Vector a_x = Vector::Zero(a->num_rows()); + Vector a_t_y = Vector::Zero(a->num_cols()); + Vector ap_x = Vector::Zero(a->num_rows()); + Vector ap_t_y = Vector::Zero(a->num_cols()); + a->RightMultiplyAndAccumulate(x.data(), a_x.data()); + ap.RightMultiplyAndAccumulate(x.data(), ap_x.data()); + EXPECT_NEAR((a_x - ap_x).norm() / a_x.norm(), + 0.0, + std::numeric_limits<double>::epsilon()); + a->LeftMultiplyAndAccumulate(y.data(), a_t_y.data()); + ap.LeftMultiplyAndAccumulate(y.data(), ap_t_y.data()); + EXPECT_NEAR((a_t_y - ap_t_y).norm() / a_t_y.norm(), + 0.0, + std::numeric_limits<double>::epsilon()); + } +} + } // namespace internal } // namespace ceres
diff --git a/internal/ceres/block_structure.h b/internal/ceres/block_structure.h index 9ee10af..aaa4080 100644 --- a/internal/ceres/block_structure.h +++ b/internal/ceres/block_structure.h
@@ -81,6 +81,11 @@ explicit CompressedList(int num_cells) noexcept : cells(num_cells) {} Block block; std::vector<Cell> cells; + // Number of non-zeros in cells of this row block + int nnz{-1}; + // Number of non-zeros in cells of this and every preceeding row block in + // block-sparse matrix + int cumulative_nnz{-1}; }; using CompressedRow = CompressedList;