Matrix generation cleanup
1. Convert a CompressedRowSparseMatrix constructor which
takes a TripletSparseMatrix as input into a factory method
which allows the input to be transposed.
2. Move the random matrix creation routine for CompressedRowSparseMatrix
from being a standalone function to a static method.
3. Add a corresponding random matrix generation static method to
TripletSparseMatrix.
4. Add a new constructor to TripletSparseMatrix, which takes as input
the row, col and values arrays.
Change-Id: Iec7b184646818f432a5e6822bea3b2f3128a82aa
diff --git a/internal/ceres/compressed_row_sparse_matrix.cc b/internal/ceres/compressed_row_sparse_matrix.cc
index 9c438f6..ccb85bc 100644
--- a/internal/ceres/compressed_row_sparse_matrix.cc
+++ b/internal/ceres/compressed_row_sparse_matrix.cc
@@ -116,6 +116,43 @@
}
transpose_rows[0] = 0;
}
+
+void AddRandomBlock(const int num_rows,
+ const int num_cols,
+ const int row_block_begin,
+ const int col_block_begin,
+ std::vector<int>* rows,
+ std::vector<int>* cols,
+ std::vector<double>* values) {
+ for (int r = 0; r < num_rows; ++r) {
+ for (int c = 0; c < num_cols; ++c) {
+ rows->push_back(row_block_begin + r);
+ cols->push_back(col_block_begin + c);
+ values->push_back(RandNormal());
+ }
+ }
+}
+
+void AddRandomSymmetricBlock(const int num_rows,
+ const int row_block_begin,
+ std::vector<int>* rows,
+ std::vector<int>* cols,
+ std::vector<double>* values) {
+ for (int r = 0; r < num_rows; ++r) {
+ for (int c = r; c < num_rows; ++c) {
+ const double v = RandNormal();
+ rows->push_back(row_block_begin + r);
+ cols->push_back(row_block_begin + c);
+ values->push_back(v);
+ if (c != r) {
+ cols->push_back(row_block_begin + r);
+ rows->push_back(row_block_begin + c);
+ values->push_back(v);
+ }
+ }
+ }
+}
+
} // namespace
// This constructor gives you a semi-initialized CompressedRowSparseMatrix.
@@ -136,47 +173,70 @@
cols_.size() * sizeof(double); // NOLINT
}
-CompressedRowSparseMatrix::CompressedRowSparseMatrix(
- const TripletSparseMatrix& m) {
- num_rows_ = m.num_rows();
- num_cols_ = m.num_cols();
- storage_type_ = UNSYMMETRIC;
+CompressedRowSparseMatrix* CompressedRowSparseMatrix::FromTripletSparseMatrix(
+ const TripletSparseMatrix& input) {
+ return CompressedRowSparseMatrix::FromTripletSparseMatrix(input, false);
+}
- rows_.resize(num_rows_ + 1, 0);
- cols_.resize(m.num_nonzeros(), 0);
- values_.resize(m.max_num_nonzeros(), 0.0);
+CompressedRowSparseMatrix*
+CompressedRowSparseMatrix::FromTripletSparseMatrixTransposed(
+ const TripletSparseMatrix& input) {
+ return CompressedRowSparseMatrix::FromTripletSparseMatrix(input, true);
+}
- // index is the list of indices into the TripletSparseMatrix m.
- vector<int> index(m.num_nonzeros(), 0);
- for (int i = 0; i < m.num_nonzeros(); ++i) {
+CompressedRowSparseMatrix* CompressedRowSparseMatrix::FromTripletSparseMatrix(
+ const TripletSparseMatrix& input, bool transpose) {
+ int num_rows = input.num_rows();
+ int num_cols = input.num_cols();
+ const int* rows = input.rows();
+ const int* cols = input.cols();
+ const double* values = input.values();
+
+ if (transpose) {
+ std::swap(num_rows, num_cols);
+ std::swap(rows, cols);
+ }
+
+ // index is the list of indices into the TripletSparseMatrix input.
+ vector<int> index(input.num_nonzeros(), 0);
+ for (int i = 0; i < input.num_nonzeros(); ++i) {
index[i] = i;
}
// Sort index such that the entries of m are ordered by row and ties
// are broken by column.
- sort(index.begin(), index.end(), RowColLessThan(m.rows(), m.cols()));
+ std::sort(index.begin(), index.end(), RowColLessThan(rows, cols));
- VLOG(1) << "# of rows: " << num_rows_ << " # of columns: " << num_cols_
- << " max_num_nonzeros: " << cols_.size() << ". Allocating "
- << ((num_rows_ + 1) * sizeof(int) + // NOLINT
- cols_.size() * sizeof(int) + // NOLINT
- cols_.size() * sizeof(double)); // NOLINT
+ VLOG(1) << "# of rows: " << num_rows << " # of columns: " << num_cols
+ << " num_nonzeros: " << input.num_nonzeros() << ". Allocating "
+ << ((num_rows + 1) * sizeof(int) + // NOLINT
+ input.num_nonzeros() * sizeof(int) + // NOLINT
+ input.num_nonzeros() * sizeof(double)); // NOLINT
+
+ CompressedRowSparseMatrix* output =
+ new CompressedRowSparseMatrix(num_rows, num_cols, input.num_nonzeros());
// Copy the contents of the cols and values array in the order given
// by index and count the number of entries in each row.
- for (int i = 0; i < m.num_nonzeros(); ++i) {
+ int* output_rows = output->mutable_rows();
+ int* output_cols = output->mutable_cols();
+ double* output_values = output->mutable_values();
+
+ output_rows[0] = 0;
+ for (int i = 0; i < index.size(); ++i) {
const int idx = index[i];
- ++rows_[m.rows()[idx] + 1];
- cols_[i] = m.cols()[idx];
- values_[i] = m.values()[idx];
+ ++output_rows[rows[idx] + 1];
+ output_cols[i] = cols[idx];
+ output_values[i] = values[idx];
}
// Find the cumulative sum of the row counts.
- for (int i = 1; i < num_rows_ + 1; ++i) {
- rows_[i] += rows_[i - 1];
+ for (int i = 1; i < num_rows + 1; ++i) {
+ output_rows[i] += output_rows[i - 1];
}
- CHECK_EQ(num_nonzeros(), m.num_nonzeros());
+ CHECK_EQ(output->num_nonzeros(), input.num_nonzeros());
+ return output;
}
CompressedRowSparseMatrix::CompressedRowSparseMatrix(const double* diagonal,
@@ -585,8 +645,8 @@
CHECK_GT(product.size(), 0);
vector<int> row_nnz;
- CompressedRowSparseMatrix* matrix =
- CreateOuterProductMatrix(num_cols, storage_type, blocks, product, &row_nnz);
+ CompressedRowSparseMatrix* matrix = CreateOuterProductMatrix(
+ num_cols, storage_type, blocks, product, &row_nnz);
const vector<int>& block_offsets = matrix->block_offsets();
@@ -697,10 +757,10 @@
CompressedRowSparseMatrix*
CompressedRowSparseMatrix::CreateOuterProductMatrixAndProgram(
- const CompressedRowSparseMatrix& m,
- const CompressedRowSparseMatrix::StorageType storage_type,
- vector<int>* program) {
- CHECK(storage_type == LOWER_TRIANGULAR || storage_type == UPPER_TRIANGULAR);
+ 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.";
@@ -778,8 +838,8 @@
const CompressedRowSparseMatrix& m,
const vector<int>& program,
CompressedRowSparseMatrix* result) {
- CHECK(result->storage_type() == LOWER_TRIANGULAR ||
- result->storage_type() == UPPER_TRIANGULAR);
+ CHECK(result->storage_type() == LOWER_TRIANGULAR ||
+ result->storage_type() == UPPER_TRIANGULAR);
result->SetZero();
double* values = result->mutable_values();
const int* rows = result->rows();
@@ -796,7 +856,6 @@
#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++]) {
@@ -853,8 +912,19 @@
CHECK_EQ(cursor, program.size());
}
-CompressedRowSparseMatrix* CreateRandomCompressedRowSparseMatrix(
- const RandomMatrixOptions& options) {
+CompressedRowSparseMatrix* CompressedRowSparseMatrix::CreateRandomMatrix(
+ const CompressedRowSparseMatrix::RandomMatrixOptions& options) {
+ CHECK_GT(options.num_row_blocks, 0);
+ CHECK_GT(options.min_row_block_size, 0);
+ CHECK_GT(options.max_row_block_size, 0);
+ CHECK_LE(options.min_row_block_size, options.max_row_block_size);
+ CHECK_GT(options.num_col_blocks, 0);
+ CHECK_GT(options.min_col_block_size, 0);
+ CHECK_GT(options.max_col_block_size, 0);
+ CHECK_LE(options.min_col_block_size, options.max_col_block_size);
+ CHECK_GT(options.block_density, 0.0);
+ CHECK_LE(options.block_density, 1.0);
+
vector<int> row_blocks;
vector<int> col_blocks;
@@ -889,22 +959,26 @@
// not what the user wants, so do the matrix generation till we have
// at least one non-zero entry.
while (tsm_values.empty()) {
- int row_block_begin = 0;
crsb_rows.clear();
crsb_cols.clear();
+ tsm_rows.clear();
+ tsm_cols.clear();
+ tsm_values.clear();
+
+ int row_block_begin = 0;
for (int r = 0; r < options.num_row_blocks; ++r) {
int col_block_begin = 0;
crsb_rows.push_back(crsb_cols.size());
for (int c = 0; c < options.num_col_blocks; ++c) {
// Randomly determine if this block is present or not.
if (RandDouble() <= options.block_density) {
- for (int i = 0; i < row_blocks[r]; ++i) {
- for (int j = 0; j < col_blocks[c]; ++j) {
- tsm_rows.push_back(row_block_begin + i);
- tsm_cols.push_back(col_block_begin + j);
- tsm_values.push_back(RandNormal());
- }
- }
+ AddRandomBlock(row_blocks[r],
+ col_blocks[c],
+ row_block_begin,
+ col_block_begin,
+ &tsm_rows,
+ &tsm_cols,
+ &tsm_values);
// Add the block to the block sparse structure.
crsb_cols.push_back(c);
}
@@ -917,21 +991,17 @@
const int num_rows = std::accumulate(row_blocks.begin(), row_blocks.end(), 0);
const int num_cols = std::accumulate(col_blocks.begin(), col_blocks.end(), 0);
- const int num_nonzeros = tsm_values.size();
-
- // Create a TripletSparseMatrix
- TripletSparseMatrix tsm(num_rows, num_cols, num_nonzeros);
- std::copy(tsm_rows.begin(), tsm_rows.end(), tsm.mutable_rows());
- std::copy(tsm_cols.begin(), tsm_cols.end(), tsm.mutable_cols());
- std::copy(tsm_values.begin(), tsm_values.end(), tsm.mutable_values());
- tsm.set_num_nonzeros(num_nonzeros);
-
- // Convert the TripletSparseMatrix to a CompressedRowSparseMatrix.
- CompressedRowSparseMatrix* matrix = new CompressedRowSparseMatrix(tsm);
+ const bool kDoNotTranspose = false;
+ CompressedRowSparseMatrix* matrix =
+ CompressedRowSparseMatrix::FromTripletSparseMatrix(
+ TripletSparseMatrix(
+ num_rows, num_cols, tsm_rows, tsm_cols, tsm_values),
+ kDoNotTranspose);
(*matrix->mutable_row_blocks()) = row_blocks;
(*matrix->mutable_col_blocks()) = col_blocks;
(*matrix->mutable_crsb_rows()) = crsb_rows;
(*matrix->mutable_crsb_cols()) = crsb_cols;
+ matrix->set_storage_type(CompressedRowSparseMatrix::UNSYMMETRIC);
return matrix;
}
diff --git a/internal/ceres/compressed_row_sparse_matrix.h b/internal/ceres/compressed_row_sparse_matrix.h
index bf67737..1e26f7c 100644
--- a/internal/ceres/compressed_row_sparse_matrix.h
+++ b/internal/ceres/compressed_row_sparse_matrix.h
@@ -50,18 +50,33 @@
public:
enum StorageType {
UNSYMMETRIC,
+ // Matrix is assumed to be symmetric but only the lower triangular
+ // part of the matrix is stored.
LOWER_TRIANGULAR,
+ // Matrix is assumed to be symmetric but only the upper triangular
+ // part of the matrix is stored.
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
- // objects.
+ // Create a matrix with the same content as the TripletSparseMatrix
+ // input. We assume that input does not have any repeated
+ // entries.
//
- // We assume that m does not have any repeated entries.
- explicit CompressedRowSparseMatrix(const TripletSparseMatrix& m);
+ // The storage type of the matrix is set to UNSYMMETRIC.
+ //
+ // Caller owns the result.
+ static CompressedRowSparseMatrix* FromTripletSparseMatrix(
+ const TripletSparseMatrix& input);
+
+ // Create a matrix with the same content as the TripletSparseMatrix
+ // input transposed. We assume that input does not have any repeated
+ // entries.
+ //
+ // The storage type of the matrix is set to UNSYMMETRIC.
+ //
+ // Caller owns the result.
+ static CompressedRowSparseMatrix* FromTripletSparseMatrixTransposed(
+ const TripletSparseMatrix& input);
// Use this constructor only if you know what you are doing. This
// creates a "blank" matrix with the appropriate amount of memory
@@ -74,17 +89,20 @@
// manually, instead of going via the indirect route of first
// constructing a TripletSparseMatrix, which leads to more than
// double the peak memory usage.
+ //
+ // The storage type is set to UNSYMMETRIC.
CompressedRowSparseMatrix(int num_rows,
int num_cols,
int max_num_nonzeros);
// Build a square sparse diagonal matrix with num_rows rows and
// columns. The diagonal m(i,i) = diagonal(i);
+ //
+ // The storage type is set to UNSYMMETRIC
CompressedRowSparseMatrix(const double* diagonal, int num_rows);
- virtual ~CompressedRowSparseMatrix();
-
// SparseMatrix interface.
+ virtual ~CompressedRowSparseMatrix();
virtual void SetZero();
virtual void RightMultiply(const double* x, double* y) const;
virtual void LeftMultiply(const double* x, double* y) const;
@@ -145,10 +163,60 @@
const std::vector<int>& crsb_cols() const { return crsb_cols_; }
std::vector<int>* mutable_crsb_cols() { return &crsb_cols_; }
+ // Create a block diagonal CompressedRowSparseMatrix with the given
+ // block structure. The individual blocks are assumed to be laid out
+ // contiguously in the diagonal array, one block at a time.
+ //
+ // Caller owns the result.
static CompressedRowSparseMatrix* CreateBlockDiagonalMatrix(
const double* diagonal,
const std::vector<int>& blocks);
+ // Options struct to control the generation of random block sparse
+ // matrices in compressed row sparse format.
+ //
+ // The random matrix generation proceeds as follows.
+ //
+ // First the row and column block structure is determined by
+ // generating random row and column block sizes that lie within the
+ // given bounds.
+ //
+ // Then we walk the block structure of the resulting matrix, and with
+ // probability block_density detemine whether they are structurally
+ // zero or not. If the answer is no, then we generate entries for the
+ // block which are distributed normally.
+ struct RandomMatrixOptions {
+ RandomMatrixOptions()
+ : num_row_blocks(0),
+ min_row_block_size(0),
+ max_row_block_size(0),
+ num_col_blocks(0),
+ min_col_block_size(0),
+ max_col_block_size(0),
+ block_density(0.0) {
+ }
+
+ int num_row_blocks;
+ int min_row_block_size;
+ int max_row_block_size;
+ int num_col_blocks;
+ int min_col_block_size;
+ int max_col_block_size;
+
+ // 0 < block_density <= 1 is the probability of a block being
+ // present in the matrix. A given random matrix will not have
+ // precisely this density.
+ double block_density;
+ };
+
+ // Create a random CompressedRowSparseMatrix whose entries are
+ // normally distributed and whose structure is determined by
+ // RandomMatrixOptions.
+ //
+ // Caller owns the result.
+ static CompressedRowSparseMatrix* CreateRandomMatrix(
+ const RandomMatrixOptions& options);
+
// Compute the sparsity structure of the product m.transpose() * m
// and create a CompressedRowSparseMatrix corresponding to it.
//
@@ -180,6 +248,10 @@
CompressedRowSparseMatrix* result);
private:
+
+ static CompressedRowSparseMatrix* FromTripletSparseMatrix(
+ const TripletSparseMatrix& input, bool transpose);
+
int num_rows_;
int num_cols_;
std::vector<int> rows_;
@@ -206,45 +278,8 @@
// carry with it compressed row sparse block information.
std::vector<int> crsb_rows_;
std::vector<int> crsb_cols_;
-
- CERES_DISALLOW_COPY_AND_ASSIGN(CompressedRowSparseMatrix);
};
-// Options struct to control the generation of random block sparse
-// matrices in compressed row sparse format.
-//
-// The random matrix generation proceeds as follows.
-//
-// First the row and column block structure is determined by
-// generating random row and column block sizes that lie within the
-// given bounds.
-//
-// Then we walk the block structure of the resulting matrix, and with
-// probability block_density detemine whether they are structurally
-// zero or not. If the answer is no, then we generate entries for the
-// block which are distributed normally.
-struct RandomMatrixOptions {
- int num_row_blocks;
- int min_row_block_size;
- int max_row_block_size;
- int num_col_blocks;
- int min_col_block_size;
- int max_col_block_size;
-
- // 0 <= block_density <= 1 is the probability of a block being
- // present in the matrix. A given random matrix will not have
- // precisely this density.
- double block_density;
-};
-
-// Create a random CompressedRowSparseMatrix whose entries are
-// normally distributed and whose structure is determined by
-// RandomMatrixOptions.
-//
-// Caller owns the result.
-CompressedRowSparseMatrix* CreateRandomCompressedRowSparseMatrix(
- const RandomMatrixOptions& options);
-
} // namespace internal
} // namespace ceres
diff --git a/internal/ceres/compressed_row_sparse_matrix_test.cc b/internal/ceres/compressed_row_sparse_matrix_test.cc
index 66e2390..bb8456e 100644
--- a/internal/ceres/compressed_row_sparse_matrix_test.cc
+++ b/internal/ceres/compressed_row_sparse_matrix_test.cc
@@ -70,7 +70,7 @@
}
class CompressedRowSparseMatrixTest : public ::testing::Test {
- protected :
+ protected:
virtual void SetUp() {
scoped_ptr<LinearLeastSquaresProblem> problem(
CreateLinearLeastSquaresProblemFromId(1));
@@ -78,7 +78,7 @@
CHECK_NOTNULL(problem.get());
tsm.reset(down_cast<TripletSparseMatrix*>(problem->A.release()));
- crsm.reset(new CompressedRowSparseMatrix(*tsm));
+ crsm.reset(CompressedRowSparseMatrix::FromTripletSparseMatrix(*tsm));
num_rows = tsm->num_rows();
num_cols = tsm->num_cols();
@@ -93,9 +93,11 @@
// With all blocks of size 1, crsb_rows and crsb_cols are equivalent to
// rows and cols.
- std::copy(crsm->rows(), crsm->rows() + crsm->num_rows() + 1,
+ std::copy(crsm->rows(),
+ crsm->rows() + crsm->num_rows() + 1,
std::back_inserter(*crsm->mutable_crsb_rows()));
- std::copy(crsm->cols(), crsm->cols() + crsm->num_nonzeros(),
+ std::copy(crsm->cols(),
+ crsm->cols() + crsm->num_nonzeros(),
std::back_inserter(*crsm->mutable_crsb_cols()));
}
@@ -172,9 +174,10 @@
tsm_appendage.Resize(i, num_cols);
tsm->AppendRows(tsm_appendage);
- CompressedRowSparseMatrix crsm_appendage(tsm_appendage);
- crsm->AppendRows(crsm_appendage);
+ scoped_ptr<CompressedRowSparseMatrix> crsm_appendage(
+ CompressedRowSparseMatrix::FromTripletSparseMatrix(tsm_appendage));
+ crsm->AppendRows(*crsm_appendage);
CompareMatrices(tsm.get(), crsm.get());
}
}
@@ -219,9 +222,9 @@
EXPECT_EQ(expected_col_blocks, crsm->col_blocks());
EXPECT_EQ(crsm->crsb_cols().size(),
- pre_crsb_cols.size() + row_and_column_blocks.size());
+ pre_crsb_cols.size() + row_and_column_blocks.size());
EXPECT_EQ(crsm->crsb_rows().size(),
- pre_crsb_rows.size() + row_and_column_blocks.size());
+ pre_crsb_rows.size() + row_and_column_blocks.size());
for (int i = 0; i < row_and_column_blocks.size(); ++i) {
EXPECT_EQ(crsm->crsb_rows()[i + pre_crsb_rows.size()],
pre_crsb_rows.back() + i + 1);
@@ -234,7 +237,6 @@
EXPECT_EQ(crsm->crsb_rows(), pre_crsb_rows);
EXPECT_EQ(crsm->crsb_cols(), pre_crsb_cols);
-
}
TEST_F(CompressedRowSparseMatrixTest, ToDenseMatrix) {
@@ -278,8 +280,8 @@
}
scoped_ptr<CompressedRowSparseMatrix> matrix(
- CompressedRowSparseMatrix::CreateBlockDiagonalMatrix(
- diagonal.data(), blocks));
+ CompressedRowSparseMatrix::CreateBlockDiagonalMatrix(diagonal.data(),
+ blocks));
EXPECT_EQ(matrix->num_rows(), 5);
EXPECT_EQ(matrix->num_cols(), 5);
@@ -351,7 +353,6 @@
cols[5] = 2;
cols[6] = 5;
-
rows[2] = 7;
cols[7] = 0;
cols[8] = 1;
@@ -394,21 +395,19 @@
TEST(CompressedRowSparseMatrix, ComputeOuterProduct) {
// "Randomly generated seed."
SetRandomState(29823);
- int kMaxNumRowBlocks = 10;
- int kMaxNumColBlocks = 10;
- int kNumTrials = 10;
+ const int kMaxNumRowBlocks = 10;
+ const int kMaxNumColBlocks = 10;
+ const int kNumTrials = 10;
// Create a random matrix, compute its outer product using Eigen and
// ComputeOuterProduct. Convert both matrices to dense matrices and
// compare their upper triangular parts.
- for (int num_row_blocks = 1;
- num_row_blocks < kMaxNumRowBlocks;
+ for (int num_row_blocks = 1; num_row_blocks < kMaxNumRowBlocks;
++num_row_blocks) {
- for (int num_col_blocks = 1;
- num_col_blocks < kMaxNumColBlocks;
+ for (int num_col_blocks = 1; num_col_blocks < kMaxNumColBlocks;
++num_col_blocks) {
for (int trial = 0; trial < kNumTrials; ++trial) {
- RandomMatrixOptions options;
+ CompressedRowSparseMatrix::RandomMatrixOptions options;
options.num_row_blocks = num_row_blocks;
options.num_col_blocks = num_col_blocks;
options.min_row_block_size = 1;
@@ -426,7 +425,7 @@
VLOG(2) << "block density: " << options.block_density;
scoped_ptr<CompressedRowSparseMatrix> random_matrix(
- CreateRandomCompressedRowSparseMatrix(options));
+ CompressedRowSparseMatrix::CreateRandomMatrix(options));
Eigen::MappedSparseMatrix<double, Eigen::RowMajor> mapped_random_matrix(
random_matrix->num_rows(),
@@ -480,5 +479,60 @@
}
}
+TEST(CompressedRowSparseMatrix, FromTripletSparseMatrix) {
+ TripletSparseMatrix::RandomMatrixOptions options;
+ options.num_rows = 5;
+ options.num_cols = 7;
+ options.density = 0.5;
+
+ const int kNumTrials = 10;
+ for (int i = 0; i < kNumTrials; ++i) {
+ scoped_ptr<TripletSparseMatrix> tsm(
+ TripletSparseMatrix::CreateRandomMatrix(options));
+ scoped_ptr<CompressedRowSparseMatrix> crsm(
+ CompressedRowSparseMatrix::FromTripletSparseMatrix(*tsm));
+
+ Matrix expected;
+ tsm->ToDenseMatrix(&expected);
+ Matrix actual;
+ crsm->ToDenseMatrix(&actual);
+ EXPECT_NEAR((expected - actual).norm() / actual.norm(),
+ 0.0,
+ std::numeric_limits<double>::epsilon())
+ << "\nexpected: \n"
+ << expected << "\nactual: \n"
+ << actual;
+ }
+}
+
+TEST(CompressedRowSparseMatrix, FromTripletSparseMatrixTransposed) {
+ TripletSparseMatrix::RandomMatrixOptions options;
+ options.num_rows = 5;
+ options.num_cols = 7;
+ options.density = 0.5;
+
+ const int kNumTrials = 10;
+ for (int i = 0; i < kNumTrials; ++i) {
+ scoped_ptr<TripletSparseMatrix> tsm(
+ TripletSparseMatrix::CreateRandomMatrix(options));
+ scoped_ptr<CompressedRowSparseMatrix> crsm(
+ CompressedRowSparseMatrix::FromTripletSparseMatrixTransposed(*tsm));
+
+ Matrix tmp;
+ tsm->ToDenseMatrix(&tmp);
+ Matrix expected = tmp.transpose();
+ Matrix actual;
+ crsm->ToDenseMatrix(&actual);
+ EXPECT_NEAR((expected - actual).norm() / actual.norm(),
+ 0.0,
+ std::numeric_limits<double>::epsilon())
+ << "\nexpected: \n"
+ << expected << "\nactual: \n"
+ << actual;
+ }
+}
+
+// TODO(sameeragarwal) Add tests for the random matrix creation methods.
+
} // namespace internal
} // namespace ceres
diff --git a/internal/ceres/dynamic_compressed_row_sparse_matrix_test.cc b/internal/ceres/dynamic_compressed_row_sparse_matrix_test.cc
index d27385a..4030142 100644
--- a/internal/ceres/dynamic_compressed_row_sparse_matrix_test.cc
+++ b/internal/ceres/dynamic_compressed_row_sparse_matrix_test.cc
@@ -108,7 +108,7 @@
tsm->ToDenseMatrix(&dense_from_tsm);
ASSERT_TRUE((dense.array() == dense_from_tsm.array()).all());
- crsm.reset(new CompressedRowSparseMatrix(*tsm));
+ crsm.reset(CompressedRowSparseMatrix::FromTripletSparseMatrix(*tsm));
Matrix dense_from_crsm;
crsm->ToDenseMatrix(&dense_from_crsm);
ASSERT_TRUE((dense.array() == dense_from_crsm.array()).all());
diff --git a/internal/ceres/schur_complement_solver.cc b/internal/ceres/schur_complement_solver.cc
index 81e0d46..26c0e89 100644
--- a/internal/ceres/schur_complement_solver.cc
+++ b/internal/ceres/schur_complement_solver.cc
@@ -547,15 +547,16 @@
}
// This is an upper triangular matrix.
- CompressedRowSparseMatrix crsm(*tsm);
+ scoped_ptr<CompressedRowSparseMatrix> crsm(
+ CompressedRowSparseMatrix::FromTripletSparseMatrix(*tsm));
// Map this to a column major, lower triangular matrix.
Eigen::MappedSparseMatrix<double, Eigen::ColMajor> eigen_lhs(
- crsm.num_rows(),
- crsm.num_rows(),
- crsm.num_nonzeros(),
- crsm.mutable_rows(),
- crsm.mutable_cols(),
- crsm.mutable_values());
+ crsm->num_rows(),
+ crsm->num_rows(),
+ crsm->num_nonzeros(),
+ crsm->mutable_rows(),
+ crsm->mutable_cols(),
+ crsm->mutable_values());
event_logger.AddEvent("ToCompressedRowSparseMatrix");
// Compute symbolic factorization if one does not exist.
diff --git a/internal/ceres/triplet_sparse_matrix.cc b/internal/ceres/triplet_sparse_matrix.cc
index 8df405c..a7df060 100644
--- a/internal/ceres/triplet_sparse_matrix.cc
+++ b/internal/ceres/triplet_sparse_matrix.cc
@@ -35,6 +35,7 @@
#include "ceres/internal/eigen.h"
#include "ceres/internal/port.h"
#include "ceres/internal/scoped_ptr.h"
+#include "ceres/random.h"
#include "ceres/types.h"
#include "glog/logging.h"
@@ -69,6 +70,27 @@
AllocateMemory();
}
+TripletSparseMatrix::TripletSparseMatrix(const int num_rows,
+ const int num_cols,
+ const std::vector<int>& rows,
+ const std::vector<int>& cols,
+ const std::vector<double>& values)
+ : num_rows_(num_rows),
+ num_cols_(num_cols),
+ max_num_nonzeros_(values.size()),
+ num_nonzeros_(values.size()),
+ rows_(NULL),
+ cols_(NULL),
+ values_(NULL) {
+ // All the sizes should at least be zero
+ CHECK_GE(num_rows, 0);
+ CHECK_GE(num_cols, 0);
+ AllocateMemory();
+ std::copy(rows.begin(), rows.end(), rows_.get());
+ std::copy(cols.begin(), cols.end(), cols_.get());
+ std::copy(values.begin(), values.end(), values_.get());
+}
+
TripletSparseMatrix::TripletSparseMatrix(const TripletSparseMatrix& orig)
: SparseMatrix(),
num_rows_(orig.num_rows_),
@@ -260,5 +282,37 @@
}
}
+TripletSparseMatrix* TripletSparseMatrix::CreateRandomMatrix(
+ const TripletSparseMatrix::RandomMatrixOptions& options) {
+ CHECK_GT(options.num_rows, 0);
+ CHECK_GT(options.num_cols, 0);
+ CHECK_GT(options.density, 0.0);
+ CHECK_LE(options.density, 1.0);
+
+ std::vector<int> rows;
+ std::vector<int> cols;
+ std::vector<double> values;
+ while (rows.empty()) {
+ rows.clear();
+ cols.clear();
+ values.clear();
+ for (int r = 0; r < options.num_rows; ++r) {
+ for (int c = 0; c < options.num_cols; ++c) {
+ if (RandDouble() <= options.density) {
+ rows.push_back(r);
+ cols.push_back(c);
+ values.push_back(RandNormal());
+ }
+ }
+ }
+ if (!rows.empty()) {
+ break;
+ }
+ }
+
+ return new TripletSparseMatrix(
+ options.num_rows, options.num_cols, rows, cols, values);
+}
+
} // namespace internal
} // namespace ceres
diff --git a/internal/ceres/triplet_sparse_matrix.h b/internal/ceres/triplet_sparse_matrix.h
index f3f5370..1d82e14 100644
--- a/internal/ceres/triplet_sparse_matrix.h
+++ b/internal/ceres/triplet_sparse_matrix.h
@@ -31,6 +31,7 @@
#ifndef CERES_INTERNAL_TRIPLET_SPARSE_MATRIX_H_
#define CERES_INTERNAL_TRIPLET_SPARSE_MATRIX_H_
+#include <vector>
#include "ceres/sparse_matrix.h"
#include "ceres/internal/eigen.h"
#include "ceres/internal/scoped_ptr.h"
@@ -47,6 +48,12 @@
public:
TripletSparseMatrix();
TripletSparseMatrix(int num_rows, int num_cols, int max_num_nonzeros);
+ TripletSparseMatrix(int num_rows,
+ int num_cols,
+ const std::vector<int>& rows,
+ const std::vector<int>& cols,
+ const std::vector<double>& values);
+
explicit TripletSparseMatrix(const TripletSparseMatrix& orig);
TripletSparseMatrix& operator=(const TripletSparseMatrix& rhs);
@@ -105,6 +112,25 @@
static TripletSparseMatrix* CreateSparseDiagonalMatrix(const double* values,
int num_rows);
+ // Options struct to control the generation of random
+ // TripletSparseMatrix objects.
+ struct RandomMatrixOptions {
+ int num_rows;
+ int num_cols;
+ // 0 < density <= 1 is the probability of an entry being
+ // structurally non-zero. A given random matrix will not have
+ // precisely this density.
+ double density;
+ };
+
+ // Create a random CompressedRowSparseMatrix whose entries are
+ // normally distributed and whose structure is determined by
+ // RandomMatrixOptions.
+ //
+ // Caller owns the result.
+ static TripletSparseMatrix* CreateRandomMatrix(
+ const TripletSparseMatrix::RandomMatrixOptions& options);
+
private:
void AllocateMemory();
void CopyData(const TripletSparseMatrix& orig);
diff --git a/internal/ceres/unsymmetric_linear_solver_test.cc b/internal/ceres/unsymmetric_linear_solver_test.cc
index 95797c5..a670f00 100644
--- a/internal/ceres/unsymmetric_linear_solver_test.cc
+++ b/internal/ceres/unsymmetric_linear_solver_test.cc
@@ -69,7 +69,8 @@
options.type == DENSE_NORMAL_CHOLESKY) {
transformed_A.reset(new DenseSparseMatrix(*A_));
} else if (options.type == SPARSE_NORMAL_CHOLESKY) {
- CompressedRowSparseMatrix* crsm = new CompressedRowSparseMatrix(*A_);
+ CompressedRowSparseMatrix* crsm =
+ CompressedRowSparseMatrix::FromTripletSparseMatrix(*A_);
// Add row/column blocks structure.
for (int i = 0; i < A_->num_rows(); ++i) {
crsm->mutable_row_blocks()->push_back(1);
