Add a storage type to CompressedRowSparseMatrix
By adding an enum to CompressedRowSparseMatrix, which indicates
whether the matrix is unsymmetric, upper or lower triangular
we are able to improve the readability and fix some minor
bugs in the way some matrix manipulation code was being
called.
Thank to William Rucklidge for this suggestion.
Change-Id: I355c90d11cd5d31f5a25741b0bda4fc4583e9095
diff --git a/internal/ceres/compressed_row_sparse_matrix.cc b/internal/ceres/compressed_row_sparse_matrix.cc
index 68af5f8..dceb3bc 100644
--- a/internal/ceres/compressed_row_sparse_matrix.cc
+++ b/internal/ceres/compressed_row_sparse_matrix.cc
@@ -109,6 +109,7 @@
int max_num_nonzeros) {
num_rows_ = num_rows;
num_cols_ = num_cols;
+ storage_type_ = UNSYMMETRIC;
rows_.resize(num_rows + 1, 0);
cols_.resize(max_num_nonzeros, 0);
values_.resize(max_num_nonzeros, 0.0);
@@ -124,6 +125,7 @@
const TripletSparseMatrix& m) {
num_rows_ = m.num_rows();
num_cols_ = m.num_cols();
+ storage_type_ = UNSYMMETRIC;
rows_.resize(num_rows_ + 1, 0);
cols_.resize(m.num_nonzeros(), 0);
@@ -168,6 +170,7 @@
num_rows_ = num_rows;
num_cols_ = num_rows;
+ storage_type_ = UNSYMMETRIC;
rows_.resize(num_rows + 1);
cols_.resize(num_rows);
values_.resize(num_rows);
@@ -448,6 +451,20 @@
CompressedRowSparseMatrix* transpose =
new CompressedRowSparseMatrix(num_cols_, num_rows_, num_nonzeros());
+ switch (storage_type_) {
+ case UNSYMMETRIC:
+ transpose->set_storage_type(UNSYMMETRIC);
+ break;
+ case LOWER_TRIANGULAR:
+ transpose->set_storage_type(UPPER_TRIANGULAR);
+ break;
+ case UPPER_TRIANGULAR:
+ transpose->set_storage_type(LOWER_TRIANGULAR);
+ break;
+ default:
+ LOG(FATAL) << "Unknown storage type: " << storage_type_;
+ };
+
TransposeForCompressedRowSparseStructure(num_rows(),
num_cols(),
num_nonzeros(),
@@ -524,6 +541,7 @@
// the outerproduct matrix, which is used in outer product computation.
CompressedRowSparseMatrix* CreateOuterProductMatrix(
const int num_cols,
+ const CompressedRowSparseMatrix::StorageType storage_type,
const vector<int>& blocks,
const vector<ProductTerm>& product,
vector<int>* row_nnz) {
@@ -546,6 +564,7 @@
CompressedRowSparseMatrix* matrix =
new CompressedRowSparseMatrix(num_cols, num_cols, num_nonzeros);
+ matrix->set_storage_type(storage_type);
// Compute block offsets for outer product matrix, which is used
// in ComputeOuterProduct.
@@ -561,6 +580,7 @@
CompressedRowSparseMatrix* CompressAndFillProgram(
const int num_cols,
+ const CompressedRowSparseMatrix::StorageType storage_type,
const vector<int>& blocks,
const vector<ProductTerm>& product,
vector<int>* program) {
@@ -568,7 +588,7 @@
vector<int> row_nnz;
CompressedRowSparseMatrix* matrix =
- CreateOuterProductMatrix(num_cols, blocks, product, &row_nnz);
+ CreateOuterProductMatrix(num_cols, storage_type, blocks, product, &row_nnz);
const vector<int>& block_offsets = matrix->block_offsets();
@@ -679,7 +699,10 @@
CompressedRowSparseMatrix*
CompressedRowSparseMatrix::CreateOuterProductMatrixAndProgram(
- const CompressedRowSparseMatrix& m, const int stype, vector<int>* program) {
+ const CompressedRowSparseMatrix& m,
+ const CompressedRowSparseMatrix::StorageType storage_type,
+ vector<int>* program) {
+ CHECK(storage_type == LOWER_TRIANGULAR || storage_type == UPPER_TRIANGULAR);
CHECK_NOTNULL(program)->clear();
CHECK_GT(m.num_nonzeros(), 0)
<< "Congratulations, you found a bug in Ceres. Please report it.";
@@ -701,7 +724,7 @@
for (int row_block = 1; row_block < crsb_rows.size(); ++row_block) {
for (int idx1 = crsb_rows[row_block - 1]; idx1 < crsb_rows[row_block];
++idx1) {
- if (stype > 0) { // Lower triangular matrix.
+ if (storage_type == LOWER_TRIANGULAR) {
for (int idx2 = crsb_rows[row_block - 1]; idx2 <= idx1; ++idx2) {
product.push_back(
ProductTerm(crsb_cols[idx1], crsb_cols[idx2], product.size()));
@@ -716,7 +739,8 @@
}
sort(product.begin(), product.end());
- return CompressAndFillProgram(m.num_cols(), col_blocks, product, program);
+ return CompressAndFillProgram(
+ m.num_cols(), storage_type, col_blocks, product, program);
}
// Give input matrix m in Compressed Row Sparse Block format
@@ -754,9 +778,10 @@
// So there is no special handling for diagonal blocks.
void CompressedRowSparseMatrix::ComputeOuterProduct(
const CompressedRowSparseMatrix& m,
- const int stype,
const vector<int>& program,
CompressedRowSparseMatrix* result) {
+ CHECK(result->storage_type() == LOWER_TRIANGULAR ||
+ result->storage_type() == UPPER_TRIANGULAR);
result->SetZero();
double* values = result->mutable_values();
const int* rows = result->rows();
@@ -769,10 +794,11 @@
const vector<int>& col_blocks = m.col_blocks();
const vector<int>& crsb_rows = m.crsb_rows();
const vector<int>& crsb_cols = m.crsb_cols();
-
+ const StorageType storage_type = result->storage_type();
#define COL_BLOCK1 (crsb_cols[idx1])
#define COL_BLOCK2 (crsb_cols[idx2])
+
// Iterate row blocks.
for (int row_block = 0, m_row_begin = 0; row_block < row_blocks.size();
m_row_begin += row_blocks[row_block++]) {
@@ -790,8 +816,7 @@
// outerproduct matrix compressed row.
const int row_begin = block_offsets[COL_BLOCK1];
const int row_nnz = rows[row_begin + 1] - rows[row_begin];
-
- if (stype > 0) { // Lower triangular matrix.
+ if (storage_type == LOWER_TRIANGULAR) {
for (int idx2 = crsb_rows[row_block], m_col_nnz2 = 0; idx2 <= idx1;
m_col_nnz2 += col_blocks[COL_BLOCK2], ++idx2, ++cursor) {
int col_nnz = program[cursor];
@@ -805,7 +830,7 @@
row_nnz,
values + rows[row_begin] + col_nnz);
}
- } else { // Upper triangular matrix.
+ } else {
for (int idx2 = idx1, m_col_nnz2 = m_col_nnz1;
idx2 < crsb_rows[row_block + 1];
m_col_nnz2 += col_blocks[COL_BLOCK2], ++idx2, ++cursor) {
diff --git a/internal/ceres/compressed_row_sparse_matrix.h b/internal/ceres/compressed_row_sparse_matrix.h
index 0bc5d01..9759d13 100644
--- a/internal/ceres/compressed_row_sparse_matrix.h
+++ b/internal/ceres/compressed_row_sparse_matrix.h
@@ -48,6 +48,13 @@
class CompressedRowSparseMatrix : public SparseMatrix {
public:
+ enum StorageType {
+ UNSYMMETRIC,
+ LOWER_TRIANGULAR,
+ UPPER_TRIANGULAR
+ };
+
+
// Build a matrix with the same content as the TripletSparseMatrix
// m. TripletSparseMatrix objects are easier to construct
// incrementally, so we use them to initialize SparseMatrix
@@ -102,6 +109,18 @@
void ToCRSMatrix(CRSMatrix* matrix) const;
+ void SolveLowerTriangularInPlace(double* solution) const;
+ void SolveLowerTriangularTransposeInPlace(double* solution) const;
+
+ CompressedRowSparseMatrix* Transpose() const;
+
+ // Destructive array resizing method.
+ void SetMaxNumNonZeros(int num_nonzeros);
+
+ // Non-destructive array resizing method.
+ void set_num_rows(const int num_rows) { num_rows_ = num_rows; }
+ void set_num_cols(const int num_cols) { num_cols_ = num_cols; }
+
// Low level access methods that expose the structure of the matrix.
const int* cols() const { return &cols_[0]; }
int* mutable_cols() { return &cols_[0]; }
@@ -109,6 +128,11 @@
const int* rows() const { return &rows_[0]; }
int* mutable_rows() { return &rows_[0]; }
+ const StorageType& storage_type() const { return storage_type_; }
+ void set_storage_type(const StorageType& storage_type) {
+ storage_type_ = storage_type;
+ }
+
const std::vector<int>& row_blocks() const { return row_blocks_; }
std::vector<int>* mutable_row_blocks() { return &row_blocks_; }
@@ -124,18 +148,6 @@
const std::vector<int>& crsb_cols() const { return crsb_cols_; }
std::vector<int>* mutable_crsb_cols() { return &crsb_cols_; }
- // Destructive array resizing method.
- void SetMaxNumNonZeros(int num_nonzeros);
-
- // Non-destructive array resizing method.
- void set_num_rows(const int num_rows) { num_rows_ = num_rows; }
- void set_num_cols(const int num_cols) { num_cols_ = num_cols; }
-
- void SolveLowerTriangularInPlace(double* solution) const;
- void SolveLowerTriangularTransposeInPlace(double* solution) const;
-
- CompressedRowSparseMatrix* Transpose() const;
-
static CompressedRowSparseMatrix* CreateBlockDiagonalMatrix(
const double* diagonal,
const std::vector<int>& blocks);
@@ -155,12 +167,11 @@
// one row per row block. The ComputeOuterProduct function reuses
// this information for each row in the row block.
//
- // stype controls the result matrix in upper or lower triangular matrix
- // stype = 1: lower triangular matrix
- // stype = -1: higher triangular matrix
+ // storage_type controls the form of the output matrix. It can be
+ // LOWER_TRIANGULAR or UPPER_TRIANGULAR.
static CompressedRowSparseMatrix* CreateOuterProductMatrixAndProgram(
const CompressedRowSparseMatrix& m,
- const int stype,
+ const StorageType storage_type,
std::vector<int>* program);
// Compute the values array for the expression m.transpose() * m,
@@ -168,7 +179,6 @@
// created using the CreateOuterProductMatrixAndProgram function
// above.
static void ComputeOuterProduct(const CompressedRowSparseMatrix& m,
- const int stype,
const std::vector<int>& program,
CompressedRowSparseMatrix* result);
@@ -178,6 +188,7 @@
std::vector<int> rows_;
std::vector<int> cols_;
std::vector<double> values_;
+ StorageType storage_type_;
// If the matrix has an underlying block structure, then it can also
// carry with it row and column block sizes. This is auxilliary and
@@ -187,11 +198,11 @@
std::vector<int> row_blocks_;
std::vector<int> col_blocks_;
- // For outerproduct matrix (J' * J), we pre-compute its block offsets
- // information here for fast outerproduct computation in block unit.
- // Since the outerproduct matrix is symmetric, we do not need to
- // distinguish row or col block. In another word, this is the prefix
- // sum of row_blocks_/col_blocks_.
+ // For outer product matrix (J' * J), we pre-compute its block
+ // offsets information here for fast outer product computation in
+ // block unit. Since the outer product matrix is symmetric, we do
+ // not need to distinguish row or col block. In another word, this
+ // is the prefix sum of row_blocks_/col_blocks_.
std::vector<int> block_offsets_;
// If the matrix has an underlying block structure, then it can also
diff --git a/internal/ceres/compressed_row_sparse_matrix_test.cc b/internal/ceres/compressed_row_sparse_matrix_test.cc
index b6e732c..cc5ed08 100644
--- a/internal/ceres/compressed_row_sparse_matrix_test.cc
+++ b/internal/ceres/compressed_row_sparse_matrix_test.cc
@@ -494,16 +494,16 @@
Matrix expected_outer_product =
mapped_random_matrix.transpose() * mapped_random_matrix;
- // Use compressed row lower triangular matrix.
- const int stype = 1;
+ // Use compressed row lower triangular matrix, which will then
+ // get mapped to a compressed column upper triangular matrix.
vector<int> program;
scoped_ptr<CompressedRowSparseMatrix> outer_product(
CompressedRowSparseMatrix::CreateOuterProductMatrixAndProgram(
- *random_matrix, stype, &program));
- CompressedRowSparseMatrix::ComputeOuterProduct(*random_matrix,
- stype,
- program,
- outer_product.get());
+ *random_matrix,
+ CompressedRowSparseMatrix::LOWER_TRIANGULAR,
+ &program));
+ CompressedRowSparseMatrix::ComputeOuterProduct(
+ *random_matrix, program, outer_product.get());
Matrix actual_outer_product =
Eigen::MappedSparseMatrix<double, Eigen::ColMajor>(
diff --git a/internal/ceres/sparse_normal_cholesky_solver.cc b/internal/ceres/sparse_normal_cholesky_solver.cc
index 8401202..64fd080 100644
--- a/internal/ceres/sparse_normal_cholesky_solver.cc
+++ b/internal/ceres/sparse_normal_cholesky_solver.cc
@@ -275,21 +275,18 @@
&event_logger);
}
- // Compute outerproduct to compressed row lower triangular matrix.
- // Eigen SimplicialLDLT default uses lower triangular part of matrix.
- // This can change to upper triangular matrix if specifying
- // Eigen::SimplicialLDLT< _MatrixType, _UpLo, _Ordering >
- // with _UpLo = Upper.
- const int stype = 1;
-
+ // Compute outer product as a compressed row lower triangular
+ // matrix, because after mapping to a column major matrix, this will
+ // become a compressed column upper triangular matrix. Which is the
+ // default that Eigen uses.
if (outer_product_.get() == NULL) {
outer_product_.reset(
CompressedRowSparseMatrix::CreateOuterProductMatrixAndProgram(
- *A, stype, &pattern_));
+ *A, CompressedRowSparseMatrix::LOWER_TRIANGULAR, &pattern_));
}
CompressedRowSparseMatrix::ComputeOuterProduct(
- *A, stype, pattern_, outer_product_.get());
+ *A, pattern_, outer_product_.get());
// Map to an upper triangular column major matrix.
//
@@ -369,23 +366,21 @@
summary.termination_type = LINEAR_SOLVER_SUCCESS;
summary.message = "Success.";
- // Compute outerproduct to compressed row lower triangular matrix.
- // CXSparse Cholesky factorization uses lower triangular part of the matrix.
- const int stype = 1;
- // Compute the normal equations. J'J delta = J'f and solve them
- // using a sparse Cholesky factorization. Notice that we explicitly
- // compute the normal equations before they can be factorized.
+ // Compute outer product as a compressed row lower triangular
+ // matrix, which would mapped to a compressed column upper
+ // triangular matrix, which is the representation used by CXSparse's
+ // sparse Cholesky factorization.
if (outer_product_.get() == NULL) {
outer_product_.reset(
CompressedRowSparseMatrix::CreateOuterProductMatrixAndProgram(
- *A, stype, &pattern_));
+ *A, CompressedRowSparseMatrix::LOWER_TRIANGULAR, &pattern_));
}
CompressedRowSparseMatrix::ComputeOuterProduct(
- *A, stype, pattern_, outer_product_.get());
- cs_di lhs =
- cxsparse_.CreateSparseMatrixTransposeView(outer_product_.get());
+ *A, pattern_, outer_product_.get());
+
+ cs_di lhs = cxsparse_.CreateSparseMatrixTransposeView(outer_product_.get());
event_logger.AddEvent("Setup");
@@ -439,26 +434,22 @@
summary.num_iterations = 1;
summary.message = "Success.";
- // Compute outerproduct to compressed row upper triangular matrix.
- // This is the fastest option for the our default natural ordering
- // (see comment in cholmod_factorize.c:205 in SuiteSparse).
- const int stype = -1;
-
- // Compute the normal equations. J'J delta = J'f and solve them
- // using a sparse Cholesky factorization. Notice that we explicitly
- // compute the normal equations before they can be factorized.
+ // Compute outer product to compressed row upper triangular matrix,
+ // this will be mapped to a compressed-column lower triangular
+ // matrix, which is the fastest option for our default natural
+ // ordering (see comment in cholmod_factorize.c:205 in SuiteSparse).
if (outer_product_.get() == NULL) {
outer_product_.reset(
CompressedRowSparseMatrix::CreateOuterProductMatrixAndProgram(
- *A, stype, &pattern_));
+ *A, CompressedRowSparseMatrix::UPPER_TRIANGULAR, &pattern_));
}
CompressedRowSparseMatrix::ComputeOuterProduct(
- *A, stype, pattern_, outer_product_.get());
+ *A, pattern_, outer_product_.get());
const int num_cols = A->num_cols();
cholmod_sparse lhs =
- ss_.CreateSparseMatrixTransposeView(outer_product_.get(), stype);
+ ss_.CreateSparseMatrixTransposeView(outer_product_.get());
event_logger.AddEvent("Setup");
if (options_.dynamic_sparsity) {
diff --git a/internal/ceres/suitesparse.cc b/internal/ceres/suitesparse.cc
index b136ce2..1a0927b 100644
--- a/internal/ceres/suitesparse.cc
+++ b/internal/ceres/suitesparse.cc
@@ -96,7 +96,7 @@
}
cholmod_sparse SuiteSparse::CreateSparseMatrixTransposeView(
- CompressedRowSparseMatrix* A, const int stype) {
+ CompressedRowSparseMatrix* A) {
cholmod_sparse m;
m.nrow = A->num_cols();
m.ncol = A->num_rows();
@@ -106,7 +106,15 @@
m.i = reinterpret_cast<void*>(A->mutable_cols());
m.x = reinterpret_cast<void*>(A->mutable_values());
m.z = NULL;
- m.stype = stype;
+
+ if (A->storage_type() == CompressedRowSparseMatrix::LOWER_TRIANGULAR) {
+ m.stype = 1;
+ } else if (A->storage_type() == CompressedRowSparseMatrix::UPPER_TRIANGULAR) {
+ m.stype = -1;
+ } else {
+ m.stype = 0;
+ }
+
m.itype = CHOLMOD_INT;
m.xtype = CHOLMOD_REAL;
m.dtype = CHOLMOD_DOUBLE;
diff --git a/internal/ceres/suitesparse.h b/internal/ceres/suitesparse.h
index 380c16b..380d76e 100644
--- a/internal/ceres/suitesparse.h
+++ b/internal/ceres/suitesparse.h
@@ -96,11 +96,8 @@
// Create a cholmod_sparse wrapper around the contents of A. This is
// a shallow object, which refers to the contents of A and does not
- // use the SuiteSparse machinery to allocate memory. This can
- // create a choldmod_sparse structure with different stype (upper
- // triangular: stype= -1 or lower triangular: stype = 1).
- cholmod_sparse CreateSparseMatrixTransposeView(CompressedRowSparseMatrix* A,
- const int stype);
+ // use the SuiteSparse machinery to allocate memory.
+ cholmod_sparse CreateSparseMatrixTransposeView(CompressedRowSparseMatrix* A);
// Given a vector x, build a cholmod_dense vector of size out_size
// with the first in_size entries copied from x. If x is NULL, then
