Refactor BlockJacobiPreconditioner
1. Rename BlockJacobiPreconditionet to BlockSparseJacobiPreconditioner.
2. Add CompressedRowSparseJacobiPreconditioner which is a
block Jacobi preconditioner for CompressedRowSparseMatrix objects.
3. Re-write the tests to be more comprehensive.
Change-Id: Icbc91f9ad2cefaad593c11397f8cdcf805d7e118
diff --git a/internal/ceres/block_jacobi_preconditioner.cc b/internal/ceres/block_jacobi_preconditioner.cc
index fdba2b8..779693c 100644
--- a/internal/ceres/block_jacobi_preconditioner.cc
+++ b/internal/ceres/block_jacobi_preconditioner.cc
@@ -30,6 +30,7 @@
#include "ceres/block_jacobi_preconditioner.h"
+#include "Eigen/Dense"
#include "ceres/block_random_access_diagonal_matrix.h"
#include "ceres/block_sparse_matrix.h"
#include "ceres/block_structure.h"
@@ -38,7 +39,7 @@
namespace ceres::internal {
-BlockJacobiPreconditioner::BlockJacobiPreconditioner(
+BlockSparseJacobiPreconditioner::BlockSparseJacobiPreconditioner(
const BlockSparseMatrix& A) {
const CompressedRowBlockStructure* bs = A.block_structure();
std::vector<int> blocks(bs->cols.size());
@@ -49,10 +50,10 @@
m_ = std::make_unique<BlockRandomAccessDiagonalMatrix>(blocks);
}
-BlockJacobiPreconditioner::~BlockJacobiPreconditioner() = default;
+BlockSparseJacobiPreconditioner::~BlockSparseJacobiPreconditioner() = default;
-bool BlockJacobiPreconditioner::UpdateImpl(const BlockSparseMatrix& A,
- const double* D) {
+bool BlockSparseJacobiPreconditioner::UpdateImpl(const BlockSparseMatrix& A,
+ const double* D) {
const CompressedRowBlockStructure* bs = A.block_structure();
const double* values = A.values();
m_->SetZero();
@@ -90,9 +91,91 @@
return true;
}
-void BlockJacobiPreconditioner::RightMultiplyAndAccumulate(const double* x,
- double* y) const {
- m_->RightMultiplyAndAccumulate(x, y);
+BlockCRSJacobiPreconditioner::BlockCRSJacobiPreconditioner(
+ const CompressedRowSparseMatrix& A) {
+ auto& col_blocks = A.col_blocks();
+
+ // Compute the number of non-zeros in the preconditioner. This is needed so
+ // that we can construct the CompressedRowSparseMatrix.
+ int m_nnz = 0;
+ for (int col_block_size : col_blocks) {
+ m_nnz += col_block_size * col_block_size;
+ }
+
+ m_ = std::make_unique<CompressedRowSparseMatrix>(
+ A.num_cols(), A.num_cols(), m_nnz);
+
+ const int num_col_blocks = col_blocks.size();
+
+ // Populate the sparsity structure of the preconditioner matrix.
+ int* m_cols = m_->mutable_cols();
+ int* m_rows = m_->mutable_rows();
+ m_rows[0] = 0;
+ for (int i = 0, col = 0, idx = 0; i < num_col_blocks; ++i) {
+ // For each column block populate a diagonal block in the preconditioner.
+ // Not that the because of the way the CompressedRowSparseMatrix format
+ // works, the entire diagonal block is laid out contiguously in memory as a
+ // row-major matrix. We will use this when updating the block.
+ const int col_block_size = col_blocks[i];
+ for (int j = 0; j < col_block_size; ++j) {
+ for (int k = 0; k < col_block_size; ++k, ++idx) {
+ m_cols[idx] = col + k;
+ }
+ m_rows[col + j + 1] = idx;
+ }
+ col += col_block_size;
+ }
+
+ CHECK_EQ(m_rows[A.num_cols()], m_nnz);
+}
+
+BlockCRSJacobiPreconditioner::~BlockCRSJacobiPreconditioner() = default;
+
+bool BlockCRSJacobiPreconditioner::UpdateImpl(
+ const CompressedRowSparseMatrix& A, const double* D) {
+ const auto& col_blocks = A.col_blocks();
+ const int num_col_blocks = col_blocks.size();
+
+ const int* a_rows = A.rows();
+ const int* a_cols = A.cols();
+ const double* a_values = A.values();
+
+ m_->SetZero();
+ double* m_values = m_->mutable_values();
+ const int* m_rows = m_->rows();
+
+ const int num_rows = A.num_rows();
+ // The following loop can likely be optimized by exploiting the fact that each
+ // row block has exactly the same sparsity structure.
+ for (int r = 0; r < num_rows; ++r) {
+ int idx = a_rows[r];
+ while (idx < a_rows[r + 1]) {
+ const int col = a_cols[idx];
+ const int col_block_size = m_rows[col + 1] - m_rows[col];
+ // We make use of the fact that the entire diagonal block is stored
+ // contiguously in memory as a row-major matrix.
+ MatrixRef m(m_values + m_rows[col], col_block_size, col_block_size);
+ ConstVectorRef b(a_values + idx, col_block_size);
+ m.selfadjointView<Eigen::Upper>().rankUpdate(b);
+ idx += col_block_size;
+ }
+ }
+
+ for (int i = 0, col = 0; i < num_col_blocks; ++i) {
+ const int col_block_size = m_rows[col + 1] - m_rows[col];
+ MatrixRef m(m_values + m_rows[col], col_block_size, col_block_size);
+
+ if (D != nullptr) {
+ m.diagonal() +=
+ ConstVectorRef(D + col, col_block_size).array().square().matrix();
+ }
+
+ m = m.selfadjointView<Eigen::Upper>().llt().solve(
+ Matrix::Identity(col_block_size, col_block_size));
+ col += col_block_size;
+ }
+
+ return true;
}
} // namespace ceres::internal
diff --git a/internal/ceres/block_jacobi_preconditioner.h b/internal/ceres/block_jacobi_preconditioner.h
index 7728eb9..265b7d6 100644
--- a/internal/ceres/block_jacobi_preconditioner.h
+++ b/internal/ceres/block_jacobi_preconditioner.h
@@ -41,30 +41,22 @@
namespace ceres::internal {
class BlockSparseMatrix;
-struct CompressedRowBlockStructure;
+class CompressedRowSparseMatrix;
// A block Jacobi preconditioner. This is intended for use with
-// conjugate gradients, or other iterative symmetric solvers. To use
-// the preconditioner, create one by passing a BlockSparseMatrix "A"
-// to the constructor. This fixes the sparsity pattern to the pattern
-// of the matrix A^TA.
-//
-// Before each use of the preconditioner in a solve with conjugate gradients,
-// update the matrix by running Update(A, D). The values of the matrix A are
-// inspected to construct the preconditioner. The vector D is applied as the
-// D^TD diagonal term.
-class CERES_NO_EXPORT BlockJacobiPreconditioner
+// conjugate gradients, or other iterative symmetric solvers.
+
+// This version of the preconditioner is for use with BlockSparseMatrix
+// Jacobians.
+class CERES_NO_EXPORT BlockSparseJacobiPreconditioner
: public BlockSparseMatrixPreconditioner {
public:
// A must remain valid while the BlockJacobiPreconditioner is.
- explicit BlockJacobiPreconditioner(const BlockSparseMatrix& A);
- BlockJacobiPreconditioner(const BlockJacobiPreconditioner&) = delete;
- void operator=(const BlockJacobiPreconditioner&) = delete;
-
- ~BlockJacobiPreconditioner() override;
-
- // Preconditioner interface
- void RightMultiplyAndAccumulate(const double* x, double* y) const final;
+ explicit BlockSparseJacobiPreconditioner(const BlockSparseMatrix& A);
+ ~BlockSparseJacobiPreconditioner() override;
+ void RightMultiplyAndAccumulate(const double* x, double* y) const final {
+ return m_->RightMultiplyAndAccumulate(x, y);
+ }
int num_rows() const final { return m_->num_rows(); }
int num_cols() const final { return m_->num_rows(); }
const BlockRandomAccessDiagonalMatrix& matrix() const { return *m_; }
@@ -75,6 +67,27 @@
std::unique_ptr<BlockRandomAccessDiagonalMatrix> m_;
};
+// This version of the preconditioner is for use with CompressedRowSparseMatrix
+// Jacobians.
+class CERES_NO_EXPORT BlockCRSJacobiPreconditioner
+ : public CompressedRowSparseMatrixPreconditioner {
+ public:
+ // A must remain valid while the BlockJacobiPreconditioner is.
+ explicit BlockCRSJacobiPreconditioner(const CompressedRowSparseMatrix& A);
+ ~BlockCRSJacobiPreconditioner() override;
+ void RightMultiplyAndAccumulate(const double* x, double* y) const final {
+ m_->RightMultiplyAndAccumulate(x, y);
+ }
+ int num_rows() const final { return m_->num_rows(); }
+ int num_cols() const final { return m_->num_rows(); }
+ const CompressedRowSparseMatrix& matrix() const { return *m_; }
+
+ private:
+ bool UpdateImpl(const CompressedRowSparseMatrix& A, const double* D) final;
+
+ std::unique_ptr<CompressedRowSparseMatrix> m_;
+};
+
} // namespace ceres::internal
#include "ceres/internal/reenable_warnings.h"
diff --git a/internal/ceres/block_jacobi_preconditioner_test.cc b/internal/ceres/block_jacobi_preconditioner_test.cc
index ab69148..6e35387 100644
--- a/internal/ceres/block_jacobi_preconditioner_test.cc
+++ b/internal/ceres/block_jacobi_preconditioner_test.cc
@@ -41,58 +41,108 @@
namespace ceres::internal {
-class BlockJacobiPreconditionerTest : public ::testing::Test {
- protected:
- void SetUpFromProblemId(int problem_id) {
- std::unique_ptr<LinearLeastSquaresProblem> problem =
- CreateLinearLeastSquaresProblemFromId(problem_id);
+TEST(BlockSparseJacobiPreconditioner, _) {
+ constexpr int kNumtrials = 10;
+ BlockSparseMatrix::RandomMatrixOptions options;
+ options.num_col_blocks = 3;
+ options.min_col_block_size = 1;
+ options.max_col_block_size = 3;
- CHECK(problem != nullptr);
- A.reset(down_cast<BlockSparseMatrix*>(problem->A.release()));
- D = std::move(problem->D);
+ options.num_row_blocks = 5;
+ options.min_row_block_size = 1;
+ options.max_row_block_size = 4;
+ options.block_density = 0.25;
+ std::mt19937 prng;
- Matrix dense_a;
- A->ToDenseMatrix(&dense_a);
- dense_ata = dense_a.transpose() * dense_a;
- dense_ata += VectorRef(D.get(), A->num_cols())
- .array()
- .square()
- .matrix()
- .asDiagonal();
- }
+ for (int trial = 0; trial < kNumtrials; ++trial) {
+ auto jacobian = BlockSparseMatrix::CreateRandomMatrix(options, prng);
+ Vector diagonal = Vector::Ones(jacobian->num_cols());
+ Matrix dense_jacobian;
+ jacobian->ToDenseMatrix(&dense_jacobian);
+ Matrix hessian = dense_jacobian.transpose() * dense_jacobian;
+ hessian.diagonal() += diagonal.array().square().matrix();
- void VerifyDiagonalBlocks(const int problem_id) {
- SetUpFromProblemId(problem_id);
+ BlockSparseJacobiPreconditioner pre(*jacobian);
+ pre.Update(*jacobian, diagonal.data());
- BlockJacobiPreconditioner pre(*A);
- pre.Update(*A, D.get());
+ // The const_cast is needed to be able to call GetCell.
auto* m = const_cast<BlockRandomAccessDiagonalMatrix*>(&pre.matrix());
- EXPECT_EQ(m->num_rows(), A->num_cols());
- EXPECT_EQ(m->num_cols(), A->num_cols());
+ EXPECT_EQ(m->num_rows(), jacobian->num_cols());
+ EXPECT_EQ(m->num_cols(), jacobian->num_cols());
- const CompressedRowBlockStructure* bs = A->block_structure();
+ const CompressedRowBlockStructure* bs = jacobian->block_structure();
for (int i = 0; i < bs->cols.size(); ++i) {
const int block_size = bs->cols[i].size;
int r, c, row_stride, col_stride;
CellInfo* cell_info = m->GetCell(i, i, &r, &c, &row_stride, &col_stride);
- MatrixRef m(cell_info->values, row_stride, col_stride);
- Matrix actual_block_inverse = m.block(r, c, block_size, block_size);
- Matrix expected_block = dense_ata.block(
+ Matrix actual_block_inverse =
+ MatrixRef(cell_info->values, row_stride, col_stride)
+ .block(r, c, block_size, block_size);
+ Matrix expected_block = hessian.block(
bs->cols[i].position, bs->cols[i].position, block_size, block_size);
const double residual = (actual_block_inverse * expected_block -
Matrix::Identity(block_size, block_size))
.norm();
EXPECT_NEAR(residual, 0.0, 1e-12) << "Block: " << i;
}
+ options.num_col_blocks++;
+ options.num_row_blocks++;
}
+}
- std::unique_ptr<BlockSparseMatrix> A;
- std::unique_ptr<double[]> D;
- Matrix dense_ata;
-};
+TEST(CompressedRowSparseJacobiPreconditioner, _) {
+ constexpr int kNumtrials = 10;
+ CompressedRowSparseMatrix::RandomMatrixOptions options;
+ options.num_col_blocks = 3;
+ options.min_col_block_size = 1;
+ options.max_col_block_size = 3;
-TEST_F(BlockJacobiPreconditionerTest, SmallProblem) { VerifyDiagonalBlocks(2); }
+ options.num_row_blocks = 5;
+ options.min_row_block_size = 1;
+ options.max_row_block_size = 4;
+ options.block_density = 0.25;
+ std::mt19937 prng;
-TEST_F(BlockJacobiPreconditionerTest, LargeProblem) { VerifyDiagonalBlocks(3); }
+ for (int trial = 0; trial < kNumtrials; ++trial) {
+ auto jacobian =
+ CompressedRowSparseMatrix::CreateRandomMatrix(options, prng);
+ Vector diagonal = Vector::Ones(jacobian->num_cols());
+
+ Matrix dense_jacobian;
+ jacobian->ToDenseMatrix(&dense_jacobian);
+ Matrix hessian = dense_jacobian.transpose() * dense_jacobian;
+ hessian.diagonal() += diagonal.array().square().matrix();
+
+ BlockCRSJacobiPreconditioner pre(*jacobian);
+ pre.Update(*jacobian, diagonal.data());
+ auto& m = pre.matrix();
+
+ EXPECT_EQ(m.num_rows(), jacobian->num_cols());
+ EXPECT_EQ(m.num_cols(), jacobian->num_cols());
+
+ const auto& col_blocks = jacobian->col_blocks();
+ for (int i = 0, col = 0; i < col_blocks.size(); ++i) {
+ const int block_size = col_blocks[i];
+ int idx = m.rows()[col];
+ for (int j = 0; j < block_size; ++j) {
+ EXPECT_EQ(m.rows()[col + j + 1] - m.rows()[col + j], block_size);
+ for (int k = 0; k < block_size; ++k, ++idx) {
+ EXPECT_EQ(m.cols()[idx], col + k);
+ }
+ }
+
+ ConstMatrixRef actual_block_inverse(
+ m.values() + m.rows()[col], block_size, block_size);
+ Matrix expected_block = hessian.block(col, col, block_size, block_size);
+ const double residual = (actual_block_inverse * expected_block -
+ Matrix::Identity(block_size, block_size))
+ .norm();
+ EXPECT_NEAR(residual, 0.0, 1e-12) << "Block: " << i;
+ col += block_size;
+ }
+ options.num_col_blocks++;
+ options.num_row_blocks++;
+ }
+}
} // namespace ceres::internal
diff --git a/internal/ceres/cgnr_solver.cc b/internal/ceres/cgnr_solver.cc
index f4f9e0f..3452e64 100644
--- a/internal/ceres/cgnr_solver.cc
+++ b/internal/ceres/cgnr_solver.cc
@@ -137,7 +137,7 @@
EventLogger event_logger("CgnrSolver::Solve");
if (!preconditioner_) {
if (options_.preconditioner_type == JACOBI) {
- preconditioner_ = std::make_unique<BlockJacobiPreconditioner>(*A);
+ preconditioner_ = std::make_unique<BlockSparseJacobiPreconditioner>(*A);
} else if (options_.preconditioner_type == SUBSET) {
Preconditioner::Options preconditioner_options;
preconditioner_options.type = SUBSET;
