Updates to sparse block matrix structures to support new sparse linear solvers. * Add methods to convert TripletSparseMatrix and BlockSparseMatrix to CRSMatrix structure. * Added tests for conversion of TripletSparseMatrix and BlockSparseMatrix to CRSMatrix structure. * Added documentation on the BlockSparseMatrix structure. Change-Id: I020cfa91c301567ceeb39ff2064183c5d88c9ed5
diff --git a/internal/ceres/block_sparse_matrix.cc b/internal/ceres/block_sparse_matrix.cc index 9ba054f..22224a3 100644 --- a/internal/ceres/block_sparse_matrix.cc +++ b/internal/ceres/block_sparse_matrix.cc
@@ -36,6 +36,7 @@ #include <vector> #include "ceres/block_structure.h" +#include "ceres/crs_matrix.h" #include "ceres/internal/eigen.h" #include "ceres/random.h" #include "ceres/small_blas.h" @@ -165,6 +166,42 @@ } } +void BlockSparseMatrix::ToCRSMatrix(CRSMatrix* crs_matrix) const { + CHECK(crs_matrix != nullptr); + crs_matrix->num_rows = num_rows_; + crs_matrix->num_cols = num_cols_; + vector<int>& rows = crs_matrix->rows; + vector<int>& cols = crs_matrix->cols; + vector<double>& values = crs_matrix->values; + rows.clear(); + cols.clear(); + values.clear(); + rows.reserve(num_nonzeros_ + 1); + cols.reserve(num_nonzeros_); + values.reserve(num_nonzeros_); + + for (const auto& row_block : block_structure_->rows) { + int row_block_size = row_block.block.size; + const vector<Cell>& cells = row_block.cells; + const int row_start = row_block.block.position; + for (int r = 0; r < row_block_size; ++r) { + rows.push_back(values.size()); + for (const auto& cell : cells) { + int col_block_id = cell.block_id; + int col_block_size = block_structure_->cols[col_block_id].size; + int col_start = block_structure_->cols[col_block_id].position; + const MatrixRef m( + values_.get() + cell.position, row_block_size, col_block_size); + for (int c = 0; c < col_block_size; ++c) { + cols.push_back(c + col_start); + values.push_back(m(r, c)); + } + } + } + } + rows.push_back(values.size()); +} + void BlockSparseMatrix::ToDenseMatrix(Matrix* dense_matrix) const { CHECK(dense_matrix != nullptr);
diff --git a/internal/ceres/block_sparse_matrix.h b/internal/ceres/block_sparse_matrix.h index 987a435..da6b641 100644 --- a/internal/ceres/block_sparse_matrix.h +++ b/internal/ceres/block_sparse_matrix.h
@@ -37,6 +37,7 @@ #include <memory> #include "ceres/block_structure.h" +#include "ceres/crs_matrix.h" #include "ceres/internal/disable_warnings.h" #include "ceres/internal/eigen.h" #include "ceres/internal/export.h" @@ -74,6 +75,7 @@ void LeftMultiply(const double* x, double* y) const final; void SquaredColumnNorm(double* x) const final; void ScaleColumns(const double* scale) final; + void ToCRSMatrix(CRSMatrix* matrix) const; void ToDenseMatrix(Matrix* dense_matrix) const final; void ToTextFile(FILE* file) const final;
diff --git a/internal/ceres/block_sparse_matrix_test.cc b/internal/ceres/block_sparse_matrix_test.cc index d8cf987..7b64a7c 100644 --- a/internal/ceres/block_sparse_matrix_test.cc +++ b/internal/ceres/block_sparse_matrix_test.cc
@@ -34,6 +34,7 @@ #include <string> #include "ceres/casts.h" +#include "ceres/crs_matrix.h" #include "ceres/internal/eigen.h" #include "ceres/linear_least_squares_problems.h" #include "ceres/triplet_sparse_matrix.h" @@ -43,6 +44,92 @@ namespace ceres { namespace internal { +namespace { +template<typename T> +void CheckVectorEq(const std::vector<T>& a, const std::vector<T>& b) { + EXPECT_EQ(a.size(), b.size()); + for (int i = 0; i < a.size(); ++i) { + EXPECT_EQ(a[i], b[i]); + } +} + +std::unique_ptr<BlockSparseMatrix> CreateTestMatrixFromId(int id) { + if (id == 0) { + // Create the following block sparse matrix: + // [ 1 2 0 0 0 0 ] + // [ 3 4 0 0 0 0 ] + // [ 0 0 5 6 7 0 ] + // [ 0 0 8 9 10 0 ] + CompressedRowBlockStructure* bs = new CompressedRowBlockStructure; + bs->cols = { + // Block size 2, position 0. + Block(2, 0), + // Block size 3, position 2. + Block(3, 2), + // Block size 1, position 5. + Block(1, 5), + }; + bs->rows = { + CompressedRow(1), + CompressedRow(1) + }; + bs->rows[0].block = Block(2, 0); + bs->rows[0].cells = { Cell(0, 0) }; + + bs->rows[1].block = Block(2, 2); + bs->rows[1].cells = { Cell(1, 4) }; + std::unique_ptr<BlockSparseMatrix> m = + std::make_unique<BlockSparseMatrix>(bs); + EXPECT_NE(m, nullptr); + EXPECT_EQ(m->num_rows(), 4); + EXPECT_EQ(m->num_cols(), 6); + EXPECT_EQ(m->num_nonzeros(), 10); + double* values = m->mutable_values(); + for (int i = 0; i < 10; ++i) { + values[i] = i + 1; + } + return m; + } else if (id == 1) { + // Create the following block sparse matrix: + // [ 1 2 0 5 6 0 ] + // [ 3 4 0 7 8 0 ] + // [ 0 0 9 0 0 0 ] + CompressedRowBlockStructure* bs = new CompressedRowBlockStructure; + bs->cols = { + // Block size 2, position 0. + Block(2, 0), + // Block size 1, position 2. + Block(1, 2), + // Block size 2, position 3. + Block(2, 3), + // Block size 1, position 5. + Block(1, 5), + }; + bs->rows = { + CompressedRow(2), + CompressedRow(1) + }; + bs->rows[0].block = Block(2, 0); + bs->rows[0].cells = { Cell(0, 0), Cell(2, 4) }; + + bs->rows[1].block = Block(1, 2); + bs->rows[1].cells = { Cell(1, 8) }; + std::unique_ptr<BlockSparseMatrix> m = + std::make_unique<BlockSparseMatrix>(bs); + EXPECT_NE(m, nullptr); + EXPECT_EQ(m->num_rows(), 3); + EXPECT_EQ(m->num_cols(), 6); + EXPECT_EQ(m->num_nonzeros(), 9); + double* values = m->mutable_values(); + for (int i = 0; i < 9; ++i) { + values[i] = i + 1; + } + return m; + } + return nullptr; +} +} // namespace + class BlockSparseMatrixTest : public ::testing::Test { protected: void SetUp() final { @@ -214,5 +301,59 @@ } } +TEST(BlockSparseMatrix, ToDenseMatrix) { + { + std::unique_ptr<BlockSparseMatrix> m = CreateTestMatrixFromId(0); + Matrix m_dense; + m->ToDenseMatrix(&m_dense); + EXPECT_EQ(m_dense.rows(), 4); + EXPECT_EQ(m_dense.cols(), 6); + Matrix m_expected(4, 6); + m_expected << 1, 2, 0, 0, 0, 0, + 3, 4, 0, 0, 0, 0, + 0, 0, 5, 6, 7, 0, + 0, 0, 8, 9, 10, 0; + EXPECT_EQ(m_dense, m_expected); + } + + { + std::unique_ptr<BlockSparseMatrix> m = CreateTestMatrixFromId(1); + Matrix m_dense; + m->ToDenseMatrix(&m_dense); + EXPECT_EQ(m_dense.rows(), 3); + EXPECT_EQ(m_dense.cols(), 6); + Matrix m_expected(3, 6); + m_expected << 1, 2, 0, 5, 6, 0, + 3, 4, 0, 7, 8, 0, + 0, 0, 9, 0, 0, 0; + EXPECT_EQ(m_dense, m_expected); + } +} + +TEST(BlockSparseMatrix, ToCRSMatrix) { + { + std::unique_ptr<BlockSparseMatrix> m = CreateTestMatrixFromId(0); + CRSMatrix m_crs; + m->ToCRSMatrix(&m_crs); + std::vector<int> rows_expected = {0, 2, 4, 7, 10}; + std::vector<int> cols_expected = {0, 1, 0, 1, 2, 3, 4, 2, 3, 4}; + std::vector<double> values_expected = {1, 2, 3, 4, 5, 6, 7, 8, 9, 10}; + CheckVectorEq(rows_expected, m_crs.rows); + CheckVectorEq(cols_expected, m_crs.cols); + CheckVectorEq(values_expected, m_crs.values); + } + { + std::unique_ptr<BlockSparseMatrix> m = CreateTestMatrixFromId(1); + CRSMatrix m_crs; + m->ToCRSMatrix(&m_crs); + std::vector<int> rows_expected = {0, 4, 8, 9}; + std::vector<int> cols_expected = {0, 1, 3, 4, 0, 1, 3, 4, 2}; + std::vector<double> values_expected = {1, 2, 5, 6, 3, 4, 7, 8, 9}; + CheckVectorEq(rows_expected, m_crs.rows); + CheckVectorEq(cols_expected, m_crs.cols); + CheckVectorEq(values_expected, m_crs.values); + } +} + } // namespace internal } // namespace ceres
diff --git a/internal/ceres/block_structure.h b/internal/ceres/block_structure.h index 52b7b2f..2ba6b32 100644 --- a/internal/ceres/block_structure.h +++ b/internal/ceres/block_structure.h
@@ -82,6 +82,82 @@ using CompressedRow = CompressedList; using CompressedColumn = CompressedList; +// CompressedRowBlockStructure specifies the storage structure of a row block +// sparse matrix. +// +// Consider the following matrix A: +// A = [A_11 A_12 ... +// A_21 A_22 ... +// ... +// A_m1 A_m2 ... ] +// +// A row block sparse matrix is a matrix where the following properties hold: +// 1. The number of rows in every block A_ij and A_ik are the same. +// 2. The number of columns in every block A_ij and A_kj are the same. +// 3. The number of rows in A_ij and A_kj may be different (i != k). +// 4. The number of columns in A_ij and A_ik may be different (j != k). +// 5. Any block A_ij may be all 0s, in which case the block is not stored. +// +// The structure of the matrix is stored as follows: +// +// The `rows' array contains the following information for each row block: +// - rows[i].block.size: The number of rows in each block A_ij in the row block. +// - rows[i].block.position: The starting row in the full matrix A of the +// row block i. +// - rows[i].cells[j].block_id: The index into the `cols' array corresponding to +// the non-zero blocks A_ij. +// - rows[i].cells[j].position: The index in the `values' array for the contents +// of block A_ij. +// +// The `cols' array contains the following information for block: +// - cols[.].size: The number of columns spanned by the block. +// - cols[.].position: The starting column in the full matrix A of the block. +// +// +// Example of a row block sparse matrix: +// block_id: | 0 |1|2 |3 | +// rows[0]: [ 1 2 0 3 4 0 ] +// [ 5 6 0 7 8 0 ] +// rows[1]: [ 0 0 9 0 0 0 ] +// +// This matrix is stored as follows: +// +// There are four column blocks: +// cols[0].size = 2 +// cols[0].position = 0 +// cols[1].size = 1 +// cols[1].position = 2 +// cols[2].size = 2 +// cols[2].position = 3 +// cols[3].size = 1 +// cols[3].position = 5 + +// The first row block spans two rows, starting at row 0: +// rows[0].block.size = 2 // This row block spans two rows. +// rows[0].block.position = 0 // It starts at row 0. +// rows[0] has two cells, at column blocks 0 and 2: +// rows[0].cells[0].block_id = 0 // This cell is in column block 0. +// rows[0].cells[0].position = 0 // See below for an explanation of this. +// rows[0].cells[1].block_id = 2 // This cell is in column block 2. +// rows[0].cells[1].position = 4 // See below for an explanation of this. +// +// The second row block spans two rows, starting at row 2: +// rows[1].block.size = 1 // This row block spans one row. +// rows[1].block.position = 2 // It starts at row 2. +// rows[1] has one cell at column block 1: +// rows[1].cells[0].block_id = 1 // This cell is in column block 1. +// rows[1].cells[0].position = 8 // See below for an explanation of this. +// +// The values in each blocks are stored contiguously in row major order. +// However, there is no unique way to order the blocks -- it is usually +// optimized to promote cache coherent access, e.g. ordering it so that +// Jacobian blocks of parameters of the same type are stored nearby. +// This is one possible way to store the values of the blocks in a values array: +// values = { 1, 2, 5, 6, 3, 4, 7, 8, 9 } +// | | | | // The three blocks. +// ^ rows[0].cells[0].position = 0 +// ^ rows[0].cells[1].position = 4 +// ^ rows[1].cells[0].position = 8 struct CERES_NO_EXPORT CompressedRowBlockStructure { std::vector<Block> cols; std::vector<CompressedRow> rows;
diff --git a/internal/ceres/triplet_sparse_matrix.cc b/internal/ceres/triplet_sparse_matrix.cc index c1c7ddc..86a444e 100644 --- a/internal/ceres/triplet_sparse_matrix.cc +++ b/internal/ceres/triplet_sparse_matrix.cc
@@ -33,6 +33,8 @@ #include <algorithm> #include <memory> +#include "ceres/compressed_row_sparse_matrix.h" +#include "ceres/crs_matrix.h" #include "ceres/internal/eigen.h" #include "ceres/internal/export.h" #include "ceres/random.h" @@ -194,6 +196,11 @@ } } +void TripletSparseMatrix::ToCRSMatrix(CRSMatrix* crs_matrix) const { + CompressedRowSparseMatrix::FromTripletSparseMatrix(*this)->ToCRSMatrix( + crs_matrix); +} + void TripletSparseMatrix::ToDenseMatrix(Matrix* dense_matrix) const { dense_matrix->resize(num_rows_, num_cols_); dense_matrix->setZero(); @@ -275,6 +282,31 @@ } } +std::unique_ptr<TripletSparseMatrix> +TripletSparseMatrix::CreateFromTextFile(FILE* file) { + CHECK(file != nullptr); + int num_rows = 0; + int num_cols = 0; + std::vector<int> rows; + std::vector<int> cols; + std::vector<double> values; + while (true) { + int row, col; + double value; + if (fscanf(file, "%d %d %lf", &row, &col, &value) != 3) { + break; + } + rows.push_back(row); + cols.push_back(col); + values.push_back(value); + num_rows = std::max(num_rows, row + 1); + num_cols = std::max(num_cols, col + 1); + } + VLOG(1) << "Read " << rows.size() << " nonzeros from file."; + return std::make_unique<TripletSparseMatrix>( + num_rows, num_cols, rows, cols, values); +} + std::unique_ptr<TripletSparseMatrix> TripletSparseMatrix::CreateRandomMatrix( const TripletSparseMatrix::RandomMatrixOptions& options) { CHECK_GT(options.num_rows, 0);
diff --git a/internal/ceres/triplet_sparse_matrix.h b/internal/ceres/triplet_sparse_matrix.h index 2cf4290..a2ba1f1 100644 --- a/internal/ceres/triplet_sparse_matrix.h +++ b/internal/ceres/triplet_sparse_matrix.h
@@ -34,6 +34,7 @@ #include <memory> #include <vector> +#include "ceres/crs_matrix.h" #include "ceres/internal/disable_warnings.h" #include "ceres/internal/eigen.h" #include "ceres/internal/export.h" @@ -68,6 +69,7 @@ void LeftMultiply(const double* x, double* y) const final; void SquaredColumnNorm(double* x) const final; void ScaleColumns(const double* scale) final; + void ToCRSMatrix(CRSMatrix* matrix) const; void ToDenseMatrix(Matrix* dense_matrix) const final; void ToTextFile(FILE* file) const final; // clang-format off @@ -135,6 +137,8 @@ static std::unique_ptr<TripletSparseMatrix> CreateRandomMatrix( const TripletSparseMatrix::RandomMatrixOptions& options); + // Load a triplet sparse matrix from a text file. + static std::unique_ptr<TripletSparseMatrix> CreateFromTextFile(FILE* file); private: void AllocateMemory(); void CopyData(const TripletSparseMatrix& orig);
diff --git a/internal/ceres/triplet_sparse_matrix_test.cc b/internal/ceres/triplet_sparse_matrix_test.cc index e96c994..577d959 100644 --- a/internal/ceres/triplet_sparse_matrix_test.cc +++ b/internal/ceres/triplet_sparse_matrix_test.cc
@@ -28,6 +28,7 @@ // // Author: sameeragarwal@google.com (Sameer Agarwal) +#include "ceres/crs_matrix.h" #include "ceres/triplet_sparse_matrix.h" #include <memory> @@ -343,4 +344,42 @@ } } +TEST(TripletSparseMatrix, ToCRSMatrix) { + // Test matrix: + // [1, 2, 0, 5, 6, 0, + // 3, 4, 0, 7, 8, 0, + // 0, 0, 9, 0, 0, 0] + TripletSparseMatrix m(3, + 6, + {0, 0, 0, 0, 1, 1, 1, 1, 2}, + {0, 1, 3, 4, 0, 1, 3, 4, 2}, + {1, 2, 3, 4, 5, 6, 7, 8, 9}); + CRSMatrix m_crs; + m.ToCRSMatrix(&m_crs); + EXPECT_EQ(m_crs.num_rows, 3); + EXPECT_EQ(m_crs.num_cols, 6); + + EXPECT_EQ(m_crs.rows.size(), 4); + EXPECT_EQ(m_crs.rows[0], 0); + EXPECT_EQ(m_crs.rows[1], 4); + EXPECT_EQ(m_crs.rows[2], 8); + EXPECT_EQ(m_crs.rows[3], 9); + + EXPECT_EQ(m_crs.cols.size(), 9); + EXPECT_EQ(m_crs.cols[0], 0); + EXPECT_EQ(m_crs.cols[1], 1); + EXPECT_EQ(m_crs.cols[2], 3); + EXPECT_EQ(m_crs.cols[3], 4); + EXPECT_EQ(m_crs.cols[4], 0); + EXPECT_EQ(m_crs.cols[5], 1); + EXPECT_EQ(m_crs.cols[6], 3); + EXPECT_EQ(m_crs.cols[7], 4); + EXPECT_EQ(m_crs.cols[8], 2); + + EXPECT_EQ(m_crs.values.size(), 9); + for (int i = 0; i < 9; ++i) { + EXPECT_EQ(m_crs.values[i], i + 1); + } +} + } // namespace ceres::internal