A number of changes to BlockSparseMatrix. 1. Add BlockSparseMatrix::CreateDiagonalMatrix 2. Add BlockSparseMatrix::AppendRows 3. Add BlockSparseMatrix::DeleteRowBlocks 4. Add BlockSparseMatrix::CreateRandomMatrix 5. Add a non-default constructor to Compressedlist. Change-Id: I7cc7656616d059cef4471335f6d5b636807953e6
diff --git a/internal/ceres/block_sparse_matrix.cc b/internal/ceres/block_sparse_matrix.cc index 68d0780..9195e9e 100644 --- a/internal/ceres/block_sparse_matrix.cc +++ b/internal/ceres/block_sparse_matrix.cc
@@ -35,6 +35,7 @@ #include <vector> #include "ceres/block_structure.h" #include "ceres/internal/eigen.h" +#include "ceres/random.h" #include "ceres/small_blas.h" #include "ceres/triplet_sparse_matrix.h" #include "glog/logging.h" @@ -242,5 +243,152 @@ } } +BlockSparseMatrix* BlockSparseMatrix::CreateDiagonalMatrix( + const double* diagonal, const std::vector<Block>& column_blocks) { + // Create the block structure for the diagonal matrix. + CompressedRowBlockStructure* bs = new CompressedRowBlockStructure(); + bs->cols = column_blocks; + int position = 0; + bs->rows.resize(column_blocks.size(), CompressedRow(1)); + for (int i = 0; i < column_blocks.size(); ++i) { + CompressedRow& row = bs->rows[i]; + row.block = column_blocks[i]; + Cell& cell = row.cells[0]; + cell.block_id = i; + cell.position = position; + position += row.block.size * row.block.size; + } + + // Create the BlockSparseMatrix with the given block structure. + BlockSparseMatrix* matrix = new BlockSparseMatrix(bs); + matrix->SetZero(); + + // Fill the values array of the block sparse matrix. + double* values = matrix->mutable_values(); + for (int i = 0; i < column_blocks.size(); ++i) { + const int size = column_blocks[i].size; + for (int j = 0; j < size; ++j) { + // (j + 1) * size is compact way of accessing the (j,j) entry. + values[j * (size + 1)] = diagonal[j]; + } + diagonal += size; + values += size * size; + } + + return matrix; +} + +void BlockSparseMatrix::AppendRows(const BlockSparseMatrix& m) { + const int old_num_nonzeros = num_nonzeros_; + const int old_num_row_blocks = block_structure_->rows.size(); + const CompressedRowBlockStructure* m_bs = m.block_structure(); + block_structure_->rows.resize(old_num_row_blocks + m_bs->rows.size()); + + 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]; + 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()); + 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; + } + } + + double* new_values = new double[num_nonzeros_]; + std::copy(values_.get(), values_.get() + old_num_nonzeros, new_values); + values_.reset(new_values); + + std::copy(m.values(), + m.values() + m.num_nonzeros(), + values_.get() + old_num_nonzeros); +} + +void BlockSparseMatrix::DeleteRowBlocks(const int delta_row_blocks) { + const int num_row_blocks = block_structure_->rows.size(); + int delta_num_nonzeros = 0; + int delta_num_rows = 0; + const std::vector<Block>& column_blocks = block_structure_->cols; + for (int i = 0; i < delta_row_blocks; ++i) { + const CompressedRow& row = block_structure_->rows[num_row_blocks - i - 1]; + delta_num_rows += row.block.size; + 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; + } + } + num_nonzeros_ -= delta_num_nonzeros; + num_rows_ -= delta_num_rows; + block_structure_->rows.resize(num_row_blocks - delta_row_blocks); +} + +BlockSparseMatrix* BlockSparseMatrix::CreateRandomMatrix( + const BlockSparseMatrix::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); + + CompressedRowBlockStructure* bs = new CompressedRowBlockStructure(); + // Generate the col block structure. + int col_block_position = 0; + for (int i = 0; i < options.num_col_blocks; ++i) { + // Generate a random integer in [min_col_block_size, max_col_block_size] + const int delta_block_size = + Uniform(options.max_col_block_size - options.min_col_block_size); + const int col_block_size = options.min_col_block_size + delta_block_size; + bs->cols.push_back(Block(col_block_size, col_block_position)); + col_block_position += col_block_size; + } + + + bool matrix_has_blocks = false; + while (!matrix_has_blocks) { + LOG(INFO) << "clearing"; + bs->rows.clear(); + int row_block_position = 0; + int value_position = 0; + for (int r = 0; r < options.num_row_blocks; ++r) { + + const int delta_block_size = + Uniform(options.max_row_block_size - options.min_row_block_size); + const int row_block_size = options.min_row_block_size + delta_block_size; + bs->rows.push_back(CompressedRow()); + CompressedRow& row = bs->rows.back(); + row.block.size = row_block_size; + row.block.position = row_block_position; + row_block_position += row_block_size; + for (int c = 0; c < options.num_col_blocks; ++c) { + if (RandDouble() > options.block_density) continue; + + row.cells.push_back(Cell()); + Cell& cell = row.cells.back(); + cell.block_id = c; + cell.position = value_position; + value_position += row_block_size * bs->cols[c].size; + matrix_has_blocks = true; + } + } + } + + BlockSparseMatrix* matrix = new BlockSparseMatrix(bs); + double* values = matrix->mutable_values(); + for (int i = 0; i < matrix->num_nonzeros(); ++i) { + values[i] = RandNormal(); + } + + return matrix; +} + } // namespace internal } // namespace ceres
diff --git a/internal/ceres/block_sparse_matrix.h b/internal/ceres/block_sparse_matrix.h index 2f9afb7..1b91a75 100644 --- a/internal/ceres/block_sparse_matrix.h +++ b/internal/ceres/block_sparse_matrix.h
@@ -84,10 +84,52 @@ void ToTripletSparseMatrix(TripletSparseMatrix* matrix) const; const CompressedRowBlockStructure* 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. + void AppendRows(const BlockSparseMatrix& m); + + // Delete the bottom delta_rows_blocks. + void DeleteRowBlocks(int delta_row_blocks); + + static BlockSparseMatrix* CreateDiagonalMatrix( + const double* diagonal, + const std::vector<Block>& column_blocks); + + 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 BlockSparseMatrix whose entries are normally + // distributed and whose structure is determined by + // RandomMatrixOptions. + // + // Caller owns the result. + static BlockSparseMatrix* CreateRandomMatrix( + const RandomMatrixOptions& options); + private: int num_rows_; int num_cols_; - int max_num_nonzeros_; int num_nonzeros_; scoped_array<double> values_; scoped_ptr<CompressedRowBlockStructure> block_structure_;
diff --git a/internal/ceres/block_sparse_matrix_test.cc b/internal/ceres/block_sparse_matrix_test.cc index 226938f..b3d21d0 100644 --- a/internal/ceres/block_sparse_matrix_test.cc +++ b/internal/ceres/block_sparse_matrix_test.cc
@@ -108,5 +108,111 @@ EXPECT_LT((m_a - m_b).norm(), 1e-12); } +TEST_F(BlockSparseMatrixTest, AppendRows) { + scoped_ptr<LinearLeastSquaresProblem> problem( + CreateLinearLeastSquaresProblemFromId(2)); + scoped_ptr<BlockSparseMatrix> m( + down_cast<BlockSparseMatrix*>(problem->A.release())); + A_->AppendRows(*m); + EXPECT_EQ(A_->num_rows(), 2 * m->num_rows()); + EXPECT_EQ(A_->num_cols(), m->num_cols()); + + problem.reset(CreateLinearLeastSquaresProblemFromId(1)); + scoped_ptr<TripletSparseMatrix> m2( + down_cast<TripletSparseMatrix*>(problem->A.release())); + B_->AppendRows(*m2); + + Vector y_a = Vector::Zero(A_->num_rows()); + Vector y_b = Vector::Zero(A_->num_rows()); + for (int i = 0; i < A_->num_cols(); ++i) { + Vector x = Vector::Zero(A_->num_cols()); + x[i] = 1.0; + y_a.setZero(); + y_b.setZero(); + + A_->RightMultiply(x.data(), y_a.data()); + B_->RightMultiply(x.data(), y_b.data()); + EXPECT_LT((y_a - y_b).norm(), 1e-12); + } +} + +TEST_F(BlockSparseMatrixTest, AppendAndDeleteBlockDiagonalMatrix) { + const std::vector<Block>& column_blocks = A_->block_structure()->cols; + const int num_cols = + column_blocks.back().size + column_blocks.back().position; + Vector diagonal(num_cols); + for (int i = 0; i < num_cols; ++i) { + diagonal(i) = 2 * i * i + 1; + } + scoped_ptr<BlockSparseMatrix> appendage( + BlockSparseMatrix::CreateDiagonalMatrix(diagonal.data(), column_blocks)); + + A_->AppendRows(*appendage); + Vector y_a, y_b; + y_a.resize(A_->num_rows()); + y_b.resize(A_->num_rows()); + for (int i = 0; i < A_->num_cols(); ++i) { + Vector x = Vector::Zero(A_->num_cols()); + x[i] = 1.0; + y_a.setZero(); + y_b.setZero(); + + A_->RightMultiply(x.data(), y_a.data()); + B_->RightMultiply(x.data(), y_b.data()); + EXPECT_LT((y_a.head(B_->num_rows()) - y_b.head(B_->num_rows())).norm(), 1e-12); + Vector expected_tail = Vector::Zero(A_->num_cols()); + expected_tail(i) = diagonal(i); + EXPECT_LT((y_a.tail(A_->num_cols()) - expected_tail).norm(), 1e-12); + } + + + A_->DeleteRowBlocks(column_blocks.size()); + EXPECT_EQ(A_->num_rows(), B_->num_rows()); + EXPECT_EQ(A_->num_cols(), B_->num_cols()); + + y_a.resize(A_->num_rows()); + y_b.resize(A_->num_rows()); + for (int i = 0; i < A_->num_cols(); ++i) { + Vector x = Vector::Zero(A_->num_cols()); + x[i] = 1.0; + y_a.setZero(); + y_b.setZero(); + + A_->RightMultiply(x.data(), y_a.data()); + B_->RightMultiply(x.data(), y_b.data()); + EXPECT_LT((y_a - y_b).norm(), 1e-12); + } +} + +TEST(BlockSparseMatrix, CreateDiagonalMatrix) { + std::vector<Block> column_blocks; + column_blocks.push_back(Block(2, 0)); + column_blocks.push_back(Block(1, 2)); + column_blocks.push_back(Block(3, 3)); + const int num_cols = + column_blocks.back().size + column_blocks.back().position; + Vector diagonal(num_cols); + for (int i = 0; i < num_cols; ++i) { + diagonal(i) = 2 * i * i + 1; + } + + scoped_ptr<BlockSparseMatrix> m( + BlockSparseMatrix::CreateDiagonalMatrix(diagonal.data(), column_blocks)); + const CompressedRowBlockStructure* bs = m->block_structure(); + EXPECT_EQ(bs->cols.size(), column_blocks.size()); + for (int i = 0; i < column_blocks.size(); ++i) { + EXPECT_EQ(bs->cols[i].size, column_blocks[i].size); + EXPECT_EQ(bs->cols[i].position, column_blocks[i].position); + } + EXPECT_EQ(m->num_rows(), m->num_cols()); + Vector x = Vector::Ones(num_cols); + Vector y = Vector::Zero(num_cols); + m->RightMultiply(x.data(), y.data()); + for (int i = 0; i < num_cols; ++i) { + EXPECT_NEAR(y[i], diagonal[i], std::numeric_limits<double>::epsilon()); + } +} + + } // namespace internal } // namespace ceres
diff --git a/internal/ceres/block_structure.h b/internal/ceres/block_structure.h index 6e7003a..9d1b5e9 100644 --- a/internal/ceres/block_structure.h +++ b/internal/ceres/block_structure.h
@@ -70,6 +70,11 @@ bool CellLessThan(const Cell& lhs, const Cell& rhs); struct CompressedList { + CompressedList() {} + + // Construct a CompressedList with the cells containing num_cells + // entries. + CompressedList(int num_cells) : cells(num_cells) {} Block block; std::vector<Cell> cells; };