Add dynamic_sparsity option.
The standard sparse normal Cholesky solver assumes a fixed
sparsity pattern which is useful for a large number of problems
presented to Ceres. However, some problems are symbolically dense
but numerically sparse i.e. each residual is a function of a
large number of parameters but at any given state the residual
only depends on a sparse subset of them. For these class of
problems it is faster to re-analyse the sparsity pattern of the
jacobian at each iteration of the non-linear optimisation instead
of including all of the zero entries in the step computation.
The proposed solution adds the dynamic_sparsity option which can
be used with SPARSE_NORMAL_CHOLESKY. A
DynamicCompressedRowSparseMatrix type (which extends
CompressedRowSparseMatrix) has been introduced which allows
dynamic addition and removal of elements. A Finalize method is
provided which then consolidates the matrix so that it can be
used in place of a regular CompressedRowSparseMatrix. An
associated jacobian writer has also been provided.
Changes that were required to make this extension were adding the
SetMaxNumNonZeros method to CompressedRowSparseMatrix and adding
a JacobianFinalizer template parameter to the ProgramEvaluator.
Change-Id: Ia5a8a9523fdae8d5b027bc35e70b4611ec2a8d01
diff --git a/internal/ceres/CMakeLists.txt b/internal/ceres/CMakeLists.txt
index 87faf2b..a9ec1c2 100644
--- a/internal/ceres/CMakeLists.txt
+++ b/internal/ceres/CMakeLists.txt
@@ -59,6 +59,8 @@
dense_sparse_matrix.cc
detect_structure.cc
dogleg_strategy.cc
+ dynamic_compressed_row_jacobian_writer.cc
+ dynamic_compressed_row_sparse_matrix.cc
evaluator.cc
file.cc
gradient_checking_cost_function.cc
@@ -220,6 +222,7 @@
CERES_TEST(covariance)
CERES_TEST(dense_sparse_matrix)
CERES_TEST(dynamic_autodiff_cost_function)
+ CERES_TEST(dynamic_compressed_row_sparse_matrix)
CERES_TEST(dynamic_numeric_diff_cost_function)
CERES_TEST(evaluator)
CERES_TEST(gradient_checker)
diff --git a/internal/ceres/compressed_row_jacobian_writer.cc b/internal/ceres/compressed_row_jacobian_writer.cc
index bbadb77..ed8db14 100644
--- a/internal/ceres/compressed_row_jacobian_writer.cc
+++ b/internal/ceres/compressed_row_jacobian_writer.cc
@@ -40,6 +40,44 @@
namespace ceres {
namespace internal {
+void CompressedRowJacobianWriter::PopulateJacobianRowAndColumnBlockVectors(
+ const Program* program, CompressedRowSparseMatrix* jacobian) {
+ const vector<ParameterBlock*>& parameter_blocks =
+ program->parameter_blocks();
+ vector<int>& col_blocks = *(jacobian->mutable_col_blocks());
+ col_blocks.resize(parameter_blocks.size());
+ for (int i = 0; i < parameter_blocks.size(); ++i) {
+ col_blocks[i] = parameter_blocks[i]->LocalSize();
+ }
+
+ const vector<ResidualBlock*>& residual_blocks =
+ program->residual_blocks();
+ vector<int>& row_blocks = *(jacobian->mutable_row_blocks());
+ row_blocks.resize(residual_blocks.size());
+ for (int i = 0; i < residual_blocks.size(); ++i) {
+ row_blocks[i] = residual_blocks[i]->NumResiduals();
+ }
+}
+
+void CompressedRowJacobianWriter::GetOrderedParameterBlocks(
+ const Program* program,
+ int residual_id,
+ vector<pair<int, int> >* evaluated_jacobian_blocks) {
+ const ResidualBlock* residual_block =
+ program->residual_blocks()[residual_id];
+ const int num_parameter_blocks = residual_block->NumParameterBlocks();
+
+ for (int j = 0; j < num_parameter_blocks; ++j) {
+ const ParameterBlock* parameter_block =
+ residual_block->parameter_blocks()[j];
+ if (!parameter_block->IsConstant()) {
+ evaluated_jacobian_blocks->push_back(
+ make_pair(parameter_block->index(), j));
+ }
+ }
+ sort(evaluated_jacobian_blocks->begin(), evaluated_jacobian_blocks->end());
+}
+
SparseMatrix* CompressedRowJacobianWriter::CreateJacobian() const {
const vector<ResidualBlock*>& residual_blocks =
program_->residual_blocks();
@@ -71,7 +109,7 @@
total_num_effective_parameters,
num_jacobian_nonzeros + total_num_effective_parameters);
- // At this stage, the CompressedSparseMatrix is an invalid state. But this
+ // At this stage, the CompressedRowSparseMatrix is an invalid state. But this
// seems to be the only way to construct it without doing a memory copy.
int* rows = jacobian->mutable_rows();
int* cols = jacobian->mutable_cols();
@@ -132,22 +170,7 @@
}
CHECK_EQ(num_jacobian_nonzeros, rows[total_num_residuals]);
- // Populate the row and column block vectors for use by block
- // oriented ordering algorithms. This is useful when
- // Solver::Options::use_block_amd = true.
- const vector<ParameterBlock*>& parameter_blocks =
- program_->parameter_blocks();
- vector<int>& col_blocks = *(jacobian->mutable_col_blocks());
- col_blocks.resize(parameter_blocks.size());
- for (int i = 0; i < parameter_blocks.size(); ++i) {
- col_blocks[i] = parameter_blocks[i]->LocalSize();
- }
-
- vector<int>& row_blocks = *(jacobian->mutable_row_blocks());
- row_blocks.resize(residual_blocks.size());
- for (int i = 0; i < residual_blocks.size(); ++i) {
- row_blocks[i] = residual_blocks[i]->NumResiduals();
- }
+ PopulateJacobianRowAndColumnBlockVectors(program_, jacobian);
return jacobian;
}
@@ -164,25 +187,10 @@
const ResidualBlock* residual_block =
program_->residual_blocks()[residual_id];
- const int num_parameter_blocks = residual_block->NumParameterBlocks();
const int num_residuals = residual_block->NumResiduals();
- // It is necessary to determine the order of the jacobian blocks before
- // copying them into the CompressedRowSparseMatrix. Just because a cost
- // function uses parameter blocks 1 after 2 in its arguments does not mean
- // that the block 1 occurs before block 2 in the column layout of the
- // jacobian. Thus, determine the order by sorting the jacobian blocks by their
- // position in the state vector.
vector<pair<int, int> > evaluated_jacobian_blocks;
- for (int j = 0; j < num_parameter_blocks; ++j) {
- const ParameterBlock* parameter_block =
- residual_block->parameter_blocks()[j];
- if (!parameter_block->IsConstant()) {
- evaluated_jacobian_blocks.push_back(
- make_pair(parameter_block->index(), j));
- }
- }
- sort(evaluated_jacobian_blocks.begin(), evaluated_jacobian_blocks.end());
+ GetOrderedParameterBlocks(program_, residual_id, &evaluated_jacobian_blocks);
// Where in the current row does the jacobian for a parameter block begin.
int col_pos = 0;
diff --git a/internal/ceres/compressed_row_jacobian_writer.h b/internal/ceres/compressed_row_jacobian_writer.h
index c103165..935e650 100644
--- a/internal/ceres/compressed_row_jacobian_writer.h
+++ b/internal/ceres/compressed_row_jacobian_writer.h
@@ -39,6 +39,7 @@
namespace ceres {
namespace internal {
+class CompressedRowSparseMatrix;
class Program;
class SparseMatrix;
@@ -49,6 +50,36 @@
: program_(program) {
}
+ // `PopulateJacobianRowAndColumnBlockVectors` sets `col_blocks` and
+ // `row_blocks` for a `CompressedRowSparseMatrix`, based on the parameter
+ // block sizes and residual sizes respectively from the program. This is
+ // useful when Solver::Options::use_block_amd = true;
+ //
+ // This function is static so that it is available to other jacobian
+ // writers which use `CompressedRowSparseMatrix` (or derived types).
+ // (Jacobian writers do not fall under any type hierarchy; they only
+ // have to provide an interface as specified in program_evaluator.h).
+ static void PopulateJacobianRowAndColumnBlockVectors(
+ const Program* program,
+ CompressedRowSparseMatrix* jacobian);
+
+ // It is necessary to determine the order of the jacobian blocks before
+ // copying them into a `CompressedRowSparseMatrix` (or derived type).
+ // Just because a cost function uses parameter blocks 1 after 2 in its
+ // arguments does not mean that the block 1 occurs before block 2 in the
+ // column layout of the jacobian. Thus, `GetOrderedParameterBlocks`
+ // determines the order by sorting the jacobian blocks by their position in
+ // the state vector.
+ //
+ // This function is static so that it is available to other jacobian
+ // writers which use `CompressedRowSparseMatrix` (or derived types).
+ // (Jacobian writers do not fall under any type hierarchy; they only
+ // have to provide an interface as specified in program_evaluator.h).
+ static void GetOrderedParameterBlocks(
+ const Program* program,
+ int residual_id,
+ vector<pair<int, int> >* evaluated_jacobian_blocks);
+
// JacobianWriter interface.
// Since the compressed row matrix has different layout than that assumed by
diff --git a/internal/ceres/compressed_row_sparse_matrix.cc b/internal/ceres/compressed_row_sparse_matrix.cc
index bef98d6..7993ed6 100644
--- a/internal/ceres/compressed_row_sparse_matrix.cc
+++ b/internal/ceres/compressed_row_sparse_matrix.cc
@@ -286,6 +286,13 @@
matrix->values.resize(matrix->rows[matrix->num_rows]);
}
+void CompressedRowSparseMatrix::SetMaxNumNonZeros(int num_nonzeros) {
+ CHECK_GE(num_nonzeros, 0);
+
+ cols_.resize(num_nonzeros);
+ values_.resize(num_nonzeros);
+}
+
void CompressedRowSparseMatrix::SolveLowerTriangularInPlace(
double* solution) const {
for (int r = 0; r < num_rows_; ++r) {
diff --git a/internal/ceres/compressed_row_sparse_matrix.h b/internal/ceres/compressed_row_sparse_matrix.h
index 06b8689..a0ba7ee 100644
--- a/internal/ceres/compressed_row_sparse_matrix.h
+++ b/internal/ceres/compressed_row_sparse_matrix.h
@@ -115,6 +115,9 @@
const vector<int>& col_blocks() const { return col_blocks_; }
vector<int>* mutable_col_blocks() { return &col_blocks_; }
+ // 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; }
diff --git a/internal/ceres/dynamic_compressed_row_finalizer.h b/internal/ceres/dynamic_compressed_row_finalizer.h
new file mode 100644
index 0000000..5e6b0d8
--- /dev/null
+++ b/internal/ceres/dynamic_compressed_row_finalizer.h
@@ -0,0 +1,51 @@
+// Ceres Solver - A fast non-linear least squares minimizer
+// Copyright 2014 Google Inc. All rights reserved.
+// http://code.google.com/p/ceres-solver/
+//
+// Redistribution and use in source and binary forms, with or without
+// modification, are permitted provided that the following conditions are met:
+//
+// * Redistributions of source code must retain the above copyright notice,
+// this list of conditions and the following disclaimer.
+// * Redistributions in binary form must reproduce the above copyright notice,
+// this list of conditions and the following disclaimer in the documentation
+// and/or other materials provided with the distribution.
+// * Neither the name of Google Inc. nor the names of its contributors may be
+// used to endorse or promote products derived from this software without
+// specific prior written permission.
+//
+// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
+// AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
+// ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
+// LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
+// CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
+// SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
+// INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
+// CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
+// ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
+// POSSIBILITY OF SUCH DAMAGE.
+//
+// Author: richie.stebbing@gmail.com (Richard Stebbing)
+
+#ifndef CERES_INTERNAL_DYNAMIC_COMPRESED_ROW_FINALIZER_H_
+#define CERES_INTERNAL_DYNAMIC_COMPRESED_ROW_FINALIZER_H_
+
+#include "ceres/casts.h"
+#include "ceres/dynamic_compressed_row_sparse_matrix.h"
+
+namespace ceres {
+namespace internal {
+
+struct DynamicCompressedRowJacobianFinalizer {
+ void operator()(SparseMatrix* base_jacobian, int num_parameters) {
+ DynamicCompressedRowSparseMatrix* jacobian =
+ down_cast<DynamicCompressedRowSparseMatrix*>(base_jacobian);
+ jacobian->Finalize(num_parameters);
+ }
+};
+
+} // namespace internal
+} // namespace ceres
+
+#endif // CERES_INTERNAL_DYNAMIC_COMPRESED_ROW_FINALISER_H_
diff --git a/internal/ceres/dynamic_compressed_row_jacobian_writer.cc b/internal/ceres/dynamic_compressed_row_jacobian_writer.cc
new file mode 100644
index 0000000..2f01617
--- /dev/null
+++ b/internal/ceres/dynamic_compressed_row_jacobian_writer.cc
@@ -0,0 +1,107 @@
+// Ceres Solver - A fast non-linear least squares minimizer
+// Copyright 2014 Google Inc. All rights reserved.
+// http://code.google.com/p/ceres-solver/
+//
+// Redistribution and use in source and binary forms, with or without
+// modification, are permitted provided that the following conditions are met:
+//
+// * Redistributions of source code must retain the above copyright notice,
+// this list of conditions and the following disclaimer.
+// * Redistributions in binary form must reproduce the above copyright notice,
+// this list of conditions and the following disclaimer in the documentation
+// and/or other materials provided with the distribution.
+// * Neither the name of Google Inc. nor the names of its contributors may be
+// used to endorse or promote products derived from this software without
+// specific prior written permission.
+//
+// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
+// AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
+// ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
+// LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
+// CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
+// SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
+// INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
+// CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
+// ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
+// POSSIBILITY OF SUCH DAMAGE.
+//
+// Author: richie.stebbing@gmail.com (Richard Stebbing)
+
+#include "ceres/compressed_row_jacobian_writer.h"
+#include "ceres/dynamic_compressed_row_jacobian_writer.h"
+#include "ceres/casts.h"
+#include "ceres/dynamic_compressed_row_sparse_matrix.h"
+#include "ceres/parameter_block.h"
+#include "ceres/program.h"
+#include "ceres/residual_block.h"
+
+namespace ceres {
+namespace internal {
+
+ScratchEvaluatePreparer*
+DynamicCompressedRowJacobianWriter::CreateEvaluatePreparers(int num_threads) {
+ return ScratchEvaluatePreparer::Create(*program_, num_threads);
+}
+
+SparseMatrix* DynamicCompressedRowJacobianWriter::CreateJacobian() const {
+ // Initialize `jacobian` with zero number of `max_num_nonzeros`.
+ const int num_residuals = program_->NumResiduals();
+ const int num_effective_parameters = program_->NumEffectiveParameters();
+
+ DynamicCompressedRowSparseMatrix* jacobian =
+ new DynamicCompressedRowSparseMatrix(num_residuals,
+ num_effective_parameters,
+ 0);
+
+ CompressedRowJacobianWriter::PopulateJacobianRowAndColumnBlockVectors(
+ program_, jacobian);
+
+ return jacobian;
+}
+
+void DynamicCompressedRowJacobianWriter::Write(int residual_id,
+ int residual_offset,
+ double **jacobians,
+ SparseMatrix* base_jacobian) {
+ DynamicCompressedRowSparseMatrix* jacobian =
+ down_cast<DynamicCompressedRowSparseMatrix*>(base_jacobian);
+
+ // Get the `residual_block` of interest.
+ const ResidualBlock* residual_block =
+ program_->residual_blocks()[residual_id];
+ const int num_residuals = residual_block->NumResiduals();
+
+ vector<pair<int, int> > evaluated_jacobian_blocks;
+ CompressedRowJacobianWriter::GetOrderedParameterBlocks(
+ program_, residual_id, &evaluated_jacobian_blocks);
+
+ // `residual_offset` is the residual row in the global jacobian.
+ // Empty the jacobian rows.
+ jacobian->ClearRows(residual_offset, num_residuals);
+
+ // Iterate over each parameter block.
+ for (int i = 0; i < evaluated_jacobian_blocks.size(); ++i) {
+ const ParameterBlock* parameter_block =
+ program_->parameter_blocks()[evaluated_jacobian_blocks[i].first];
+ const int parameter_block_jacobian_index =
+ evaluated_jacobian_blocks[i].second;
+ const int parameter_block_size = parameter_block->LocalSize();
+
+ // For each parameter block only insert its non-zero entries.
+ for (int r = 0; r < num_residuals; ++r) {
+ for (int c = 0; c < parameter_block_size; ++c) {
+ const double& v = jacobians[parameter_block_jacobian_index][
+ r * parameter_block_size + c];
+ // Only insert non-zero entries.
+ if (v != 0.0) {
+ jacobian->InsertEntry(
+ residual_offset + r, parameter_block->delta_offset() + c, v);
+ }
+ }
+ }
+ }
+}
+
+} // namespace internal
+} // namespace ceres
diff --git a/internal/ceres/dynamic_compressed_row_jacobian_writer.h b/internal/ceres/dynamic_compressed_row_jacobian_writer.h
new file mode 100644
index 0000000..df9581b
--- /dev/null
+++ b/internal/ceres/dynamic_compressed_row_jacobian_writer.h
@@ -0,0 +1,83 @@
+// Ceres Solver - A fast non-linear least squares minimizer
+// Copyright 2014 Google Inc. All rights reserved.
+// http://code.google.com/p/ceres-solver/
+//
+// Redistribution and use in source and binary forms, with or without
+// modification, are permitted provided that the following conditions are met:
+//
+// * Redistributions of source code must retain the above copyright notice,
+// this list of conditions and the following disclaimer.
+// * Redistributions in binary form must reproduce the above copyright notice,
+// this list of conditions and the following disclaimer in the documentation
+// and/or other materials provided with the distribution.
+// * Neither the name of Google Inc. nor the names of its contributors may be
+// used to endorse or promote products derived from this software without
+// specific prior written permission.
+//
+// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
+// AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
+// ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
+// LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
+// CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
+// SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
+// INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
+// CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
+// ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
+// POSSIBILITY OF SUCH DAMAGE.
+//
+// Author: richie.stebbing@gmail.com (Richard Stebbing)
+//
+// A jacobian writer that directly writes to dynamic compressed row sparse
+// matrices.
+
+#ifndef CERES_INTERNAL_DYNAMIC_COMPRESSED_ROW_JACOBIAN_WRITER_H_
+#define CERES_INTERNAL_DYNAMIC_COMPRESSED_ROW_JACOBIAN_WRITER_H_
+
+#include "ceres/evaluator.h"
+#include "ceres/scratch_evaluate_preparer.h"
+
+namespace ceres {
+namespace internal {
+
+class Program;
+class SparseMatrix;
+
+class DynamicCompressedRowJacobianWriter {
+ public:
+ DynamicCompressedRowJacobianWriter(Evaluator::Options /* ignored */,
+ Program* program)
+ : program_(program) {
+ }
+
+ // JacobianWriter interface.
+
+ // The compressed row matrix has different layout than that assumed by
+ // the cost functions. The scratch space is therefore used to store
+ // the jacobians (including zeros) temporarily before only the non-zero
+ // entries are copied over to the larger jacobian in `Write`.
+ ScratchEvaluatePreparer* CreateEvaluatePreparers(int num_threads);
+
+ // Return a `DynamicCompressedRowSparseMatrix` which is filled by
+ // `Write`. Note that `Finalize` must be called to make the
+ // `CompressedRowSparseMatrix` interface valid.
+ SparseMatrix* CreateJacobian() const;
+
+ // Write only the non-zero jacobian entries for a residual block
+ // (specified by `residual_id`) into `base_jacobian`, starting at the row
+ // specifed by `residual_offset`.
+ //
+ // This method is thread-safe over residual blocks (each `residual_id`).
+ void Write(int residual_id,
+ int residual_offset,
+ double **jacobians,
+ SparseMatrix* base_jacobian);
+
+ private:
+ Program* program_;
+};
+
+} // namespace internal
+} // namespace ceres
+
+#endif // CERES_INTERNAL_DYNAMIC_COMPRESSED_ROW_JACOBIAN_WRITER_H_
diff --git a/internal/ceres/dynamic_compressed_row_sparse_matrix.cc b/internal/ceres/dynamic_compressed_row_sparse_matrix.cc
new file mode 100644
index 0000000..f285d52
--- /dev/null
+++ b/internal/ceres/dynamic_compressed_row_sparse_matrix.cc
@@ -0,0 +1,107 @@
+// Ceres Solver - A fast non-linear least squares minimizer
+// Copyright 2014 Google Inc. All rights reserved.
+// http://code.google.com/p/ceres-solver/
+//
+// Redistribution and use in source and binary forms, with or without
+// modification, are permitted provided that the following conditions are met:
+//
+// * Redistributions of source code must retain the above copyright notice,
+// this list of conditions and the following disclaimer.
+// * Redistributions in binary form must reproduce the above copyright notice,
+// this list of conditions and the following disclaimer in the documentation
+// and/or other materials provided with the distribution.
+// * Neither the name of Google Inc. nor the names of its contributors may be
+// used to endorse or promote products derived from this software without
+// specific prior written permission.
+//
+// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
+// AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
+// ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
+// LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
+// CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
+// SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
+// INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
+// CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
+// ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
+// POSSIBILITY OF SUCH DAMAGE.
+//
+// Author: richie.stebbing@gmail.com (Richard Stebbing)
+
+#include <cstring>
+#include "ceres/dynamic_compressed_row_sparse_matrix.h"
+
+namespace ceres {
+namespace internal {
+
+DynamicCompressedRowSparseMatrix::DynamicCompressedRowSparseMatrix(
+ int num_rows,
+ int num_cols,
+ int initial_max_num_nonzeros)
+ : CompressedRowSparseMatrix(num_rows,
+ num_cols,
+ initial_max_num_nonzeros) {
+ dynamic_cols_.resize(num_rows);
+ dynamic_values_.resize(num_rows);
+ }
+
+void DynamicCompressedRowSparseMatrix::InsertEntry(int row,
+ int col,
+ const double& value) {
+ CHECK_GE(row, 0);
+ CHECK_LT(row, num_rows());
+ CHECK_GE(col, 0);
+ CHECK_LT(col, num_cols());
+ dynamic_cols_[row].push_back(col);
+ dynamic_values_[row].push_back(value);
+}
+
+void DynamicCompressedRowSparseMatrix::ClearRows(int row_start,
+ int num_rows) {
+ for (int r = 0; r < num_rows; ++r) {
+ const int i = row_start + r;
+ CHECK_GE(i, 0);
+ CHECK_LT(i, this->num_rows());
+ dynamic_cols_[i].resize(0);
+ dynamic_values_[i].resize(0);
+ }
+}
+
+void DynamicCompressedRowSparseMatrix::Finalize(int num_additional_elements) {
+ // `num_additional_elements` is provided as an argument so that additional
+ // storage can be reserved when it is known by the finalizer.
+ CHECK_GE(num_additional_elements, 0);
+
+ // Count the number of non-zeros and resize `cols_` and `values_`.
+ int num_jacobian_nonzeros = 0;
+ for (int i = 0; i < dynamic_cols_.size(); ++i) {
+ num_jacobian_nonzeros += dynamic_cols_[i].size();
+ }
+
+ SetMaxNumNonZeros(num_jacobian_nonzeros + num_additional_elements);
+
+ // Flatten `dynamic_cols_` into `cols_` and `dynamic_values_`
+ // into `values_`.
+ int index_into_values_and_cols = 0;
+ for (int i = 0; i < num_rows(); ++i) {
+ mutable_rows()[i] = index_into_values_and_cols;
+ const int num_nonzero_columns = dynamic_cols_[i].size();
+ if (num_nonzero_columns > 0) {
+ memcpy(mutable_cols() + index_into_values_and_cols,
+ &dynamic_cols_[i][0],
+ dynamic_cols_[i].size() * sizeof(dynamic_cols_[0][0]));
+ memcpy(mutable_values() + index_into_values_and_cols,
+ &dynamic_values_[i][0],
+ dynamic_values_[i].size() * sizeof(dynamic_values_[0][0]));
+ index_into_values_and_cols += dynamic_cols_[i].size();
+ }
+ }
+ mutable_rows()[num_rows()] = index_into_values_and_cols;
+
+ CHECK_EQ(index_into_values_and_cols, num_jacobian_nonzeros)
+ << "Ceres bug: final index into values_ and cols_ should be equal to "
+ << "the number of jacobian nonzeros. Please contact the developers!";
+}
+
+} // namespace internal
+} // namespace ceres
diff --git a/internal/ceres/dynamic_compressed_row_sparse_matrix.h b/internal/ceres/dynamic_compressed_row_sparse_matrix.h
new file mode 100644
index 0000000..7a89a70
--- /dev/null
+++ b/internal/ceres/dynamic_compressed_row_sparse_matrix.h
@@ -0,0 +1,99 @@
+// Ceres Solver - A fast non-linear least squares minimizer
+// Copyright 2014 Google Inc. All rights reserved.
+// http://code.google.com/p/ceres-solver/
+//
+// Redistribution and use in source and binary forms, with or without
+// modification, are permitted provided that the following conditions are met:
+//
+// * Redistributions of source code must retain the above copyright notice,
+// this list of conditions and the following disclaimer.
+// * Redistributions in binary form must reproduce the above copyright notice,
+// this list of conditions and the following disclaimer in the documentation
+// and/or other materials provided with the distribution.
+// * Neither the name of Google Inc. nor the names of its contributors may be
+// used to endorse or promote products derived from this software without
+// specific prior written permission.
+//
+// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
+// AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
+// ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
+// LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
+// CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
+// SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
+// INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
+// CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
+// ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
+// POSSIBILITY OF SUCH DAMAGE.
+//
+// Author: richie.stebbing@gmail.com (Richard Stebbing)
+//
+// A compressed row sparse matrix that provides an extended interface to
+// allow dynamic insertion of entries. This is provided for the use case
+// where the sparsity structure and number of non-zero entries is dynamic.
+// This flexibility is achieved by using an (internal) scratch space that
+// allows independent insertion of entries into each row (thread-safe).
+// Once insertion is complete, the `Finalize` method must be called to ensure
+// that the underlying `CompressedRowSparseMatrix` is consistent.
+//
+// This should only be used if you really do need a dynamic sparsity pattern.
+
+#ifndef CERES_INTERNAL_DYNAMIC_COMPRESSED_ROW_SPARSE_MATRIX_H_
+#define CERES_INTERNAL_DYNAMIC_COMPRESSED_ROW_SPARSE_MATRIX_H_
+
+#include "ceres/compressed_row_sparse_matrix.h"
+
+namespace ceres {
+namespace internal {
+
+class DynamicCompressedRowSparseMatrix : public CompressedRowSparseMatrix {
+ public:
+ // Set the number of rows and columns for the underlyig
+ // `CompressedRowSparseMatrix` and set the initial number of maximum non-zero
+ // entries. Note that following the insertion of entries, when `Finalize`
+ // is called the number of non-zeros is determined and all internal
+ // structures are adjusted as required. If you know the upper limit on the
+ // number of non-zeros, then passing this value here can prevent future
+ // memory reallocations which may improve performance. Otherwise, if no
+ // upper limit is available a value of 0 is sufficient.
+ //
+ // Typical usage of this class is to define a new instance with a given
+ // number of rows, columns and maximum number of non-zero elements
+ // (if available). Next, entries are inserted at row and column positions
+ // using `InsertEntry`. Finally, once all elements have been inserted,
+ // `Finalize` must be called to make the underlying
+ // `CompressedRowSparseMatrix` consistent.
+ DynamicCompressedRowSparseMatrix(int num_rows,
+ int num_cols,
+ int initial_max_num_nonzeros);
+
+ // Insert an entry at a given row and column position. This method is
+ // thread-safe across rows i.e. different threads can insert values
+ // simultaneously into different rows. It should be emphasised that this
+ // method always inserts a new entry and does not check for existing
+ // entries at the specified row and column position. Duplicate entries
+ // for a given row and column position will result in undefined
+ // behavior.
+ void InsertEntry(int row, int col, const double& value);
+
+ // Clear all entries for rows, starting from row index `row_start`
+ // and proceeding for `num_rows`.
+ void ClearRows(int row_start, int num_rows);
+
+ // Make the underlying internal `CompressedRowSparseMatrix` data structures
+ // consistent. Additional space for non-zero entries in the
+ // `CompressedRowSparseMatrix` can be reserved by specifying
+ // `num_additional_elements`. This is useful when it is known that rows will
+ // be appended to the `CompressedRowSparseMatrix` (e.g. appending a diagonal
+ // matrix to the jacobian) as it prevents need for future reallocation.
+ void Finalize(int num_additional_elements);
+
+ private:
+ vector<vector<int> > dynamic_cols_;
+ vector<vector<double> > dynamic_values_;
+};
+
+} // namespace internal
+} // namespace ceres
+
+#endif // CERES_INTERNAL_DYNAMIC_COMPRESSED_ROW_SPARSE_MATRIX_H_
diff --git a/internal/ceres/dynamic_compressed_row_sparse_matrix_test.cc b/internal/ceres/dynamic_compressed_row_sparse_matrix_test.cc
new file mode 100644
index 0000000..03bfcb6
--- /dev/null
+++ b/internal/ceres/dynamic_compressed_row_sparse_matrix_test.cc
@@ -0,0 +1,217 @@
+// Ceres Solver - A fast non-linear least squares minimizer
+// Copyright 2014 Google Inc. All rights reserved.
+// http://code.google.com/p/ceres-solver/
+//
+// Redistribution and use in source and binary forms, with or without
+// modification, are permitted provided that the following conditions are met:
+//
+// * Redistributions of source code must retain the above copyright notice,
+// this list of conditions and the following disclaimer.
+// * Redistributions in binary form must reproduce the above copyright notice,
+// this list of conditions and the following disclaimer in the documentation
+// and/or other materials provided with the distribution.
+// * Neither the name of Google Inc. nor the names of its contributors may be
+// used to endorse or promote products derived from this software without
+// specific prior written permission.
+//
+// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
+// AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
+// ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
+// LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
+// CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
+// SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
+// INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
+// CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
+// ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
+// POSSIBILITY OF SUCH DAMAGE.
+//
+// Author: richie.stebbing@gmail.com (Richard Stebbing)
+
+#include "ceres/dynamic_compressed_row_sparse_matrix.h"
+
+#include "ceres/casts.h"
+#include "ceres/compressed_row_sparse_matrix.h"
+#include "ceres/casts.h"
+#include "ceres/internal/eigen.h"
+#include "ceres/internal/scoped_ptr.h"
+#include "ceres/linear_least_squares_problems.h"
+#include "ceres/triplet_sparse_matrix.h"
+#include "gtest/gtest.h"
+
+namespace ceres {
+namespace internal {
+
+class DynamicCompressedRowSparseMatrixTest : public ::testing::Test {
+ protected:
+ virtual void SetUp() {
+ num_rows = 7;
+ num_cols = 4;
+
+ // The number of additional elements reserved when `Finalize` is called
+ // should have no effect on the number of rows, columns or nonzeros.
+ // Set this to some nonzero value to be sure.
+ num_additional_elements = 13;
+
+ expected_num_nonzeros = num_rows * num_cols - min(num_rows, num_cols);
+
+ InitialiseDenseReference();
+ InitialiseSparseMatrixReferences();
+
+ dcrsm.reset(new DynamicCompressedRowSparseMatrix(num_rows,
+ num_cols,
+ 0));
+ }
+
+ void Finalize() {
+ dcrsm->Finalize(num_additional_elements);
+ }
+
+ void InitialiseDenseReference() {
+ dense.resize(num_rows, num_cols);
+ dense.setZero();
+ int num_nonzeros = 0;
+ for (int i = 0; i < (num_rows * num_cols); ++i) {
+ const int r = i / num_cols, c = i % num_cols;
+ if (r != c) {
+ dense(r, c) = i + 1;
+ ++num_nonzeros;
+ }
+ }
+ ASSERT_EQ(num_nonzeros, expected_num_nonzeros);
+ }
+
+ void InitialiseSparseMatrixReferences() {
+ std::vector<int> rows, cols;
+ std::vector<double> values;
+ for (int i = 0; i < (num_rows * num_cols); ++i) {
+ const int r = i / num_cols, c = i % num_cols;
+ if (r != c) {
+ rows.push_back(r);
+ cols.push_back(c);
+ values.push_back(i + 1);
+ }
+ }
+ ASSERT_EQ(values.size(), expected_num_nonzeros);
+
+ tsm.reset(new TripletSparseMatrix(num_rows,
+ num_cols,
+ expected_num_nonzeros));
+ std::copy(rows.begin(), rows.end(), tsm->mutable_rows());
+ std::copy(cols.begin(), cols.end(), tsm->mutable_cols());
+ std::copy(values.begin(), values.end(), tsm->mutable_values());
+ tsm->set_num_nonzeros(values.size());
+
+ Matrix dense_from_tsm;
+ tsm->ToDenseMatrix(&dense_from_tsm);
+ ASSERT_TRUE((dense.array() == dense_from_tsm.array()).all());
+
+ crsm.reset(new CompressedRowSparseMatrix(*tsm));
+ Matrix dense_from_crsm;
+ crsm->ToDenseMatrix(&dense_from_crsm);
+ ASSERT_TRUE((dense.array() == dense_from_crsm.array()).all());
+ }
+
+ void InsertNonZeroEntriesFromDenseReference() {
+ for (int r = 0; r < num_rows; ++r) {
+ for (int c = 0; c < num_cols; ++c) {
+ const double& v = dense(r, c);
+ if (v != 0.0) {
+ dcrsm->InsertEntry(r, c, v);
+ }
+ }
+ }
+ }
+
+ void ExpectEmpty() {
+ EXPECT_EQ(dcrsm->num_rows(), num_rows);
+ EXPECT_EQ(dcrsm->num_cols(), num_cols);
+ EXPECT_EQ(dcrsm->num_nonzeros(), 0);
+
+ Matrix dense_from_dcrsm;
+ dcrsm->ToDenseMatrix(&dense_from_dcrsm);
+ EXPECT_EQ(dense_from_dcrsm.rows(), num_rows);
+ EXPECT_EQ(dense_from_dcrsm.cols(), num_cols);
+ EXPECT_TRUE((dense_from_dcrsm.array() == 0.0).all());
+ }
+
+ void ExpectEqualToDenseReference() {
+ Matrix dense_from_dcrsm;
+ dcrsm->ToDenseMatrix(&dense_from_dcrsm);
+ EXPECT_TRUE((dense.array() == dense_from_dcrsm.array()).all());
+ }
+
+ void ExpectEqualToCompressedRowSparseMatrixReference() {
+ typedef Eigen::Map<const Eigen::VectorXi> ConstIntVectorRef;
+
+ ConstIntVectorRef crsm_rows(crsm->rows(), crsm->num_rows() + 1);
+ ConstIntVectorRef dcrsm_rows(dcrsm->rows(), dcrsm->num_rows() + 1);
+ EXPECT_TRUE((crsm_rows.array() == dcrsm_rows.array()).all());
+
+ ConstIntVectorRef crsm_cols(crsm->cols(), crsm->num_nonzeros());
+ ConstIntVectorRef dcrsm_cols(dcrsm->cols(), dcrsm->num_nonzeros());
+ EXPECT_TRUE((crsm_cols.array() == dcrsm_cols.array()).all());
+
+ ConstVectorRef crsm_values(crsm->values(), crsm->num_nonzeros());
+ ConstVectorRef dcrsm_values(dcrsm->values(), dcrsm->num_nonzeros());
+ EXPECT_TRUE((crsm_values.array() == dcrsm_values.array()).all());
+ }
+
+ int num_rows;
+ int num_cols;
+
+ int num_additional_elements;
+
+ int expected_num_nonzeros;
+
+ Matrix dense;
+ scoped_ptr<TripletSparseMatrix> tsm;
+ scoped_ptr<CompressedRowSparseMatrix> crsm;
+
+ scoped_ptr<DynamicCompressedRowSparseMatrix> dcrsm;
+};
+
+TEST_F(DynamicCompressedRowSparseMatrixTest, Initialization) {
+ ExpectEmpty();
+
+ Finalize();
+ ExpectEmpty();
+}
+
+TEST_F(DynamicCompressedRowSparseMatrixTest, InsertEntryAndFinalize) {
+ InsertNonZeroEntriesFromDenseReference();
+ ExpectEmpty();
+
+ Finalize();
+ ExpectEqualToDenseReference();
+ ExpectEqualToCompressedRowSparseMatrixReference();
+}
+
+TEST_F(DynamicCompressedRowSparseMatrixTest, ClearRows) {
+ InsertNonZeroEntriesFromDenseReference();
+ Finalize();
+ ExpectEqualToDenseReference();
+ ExpectEqualToCompressedRowSparseMatrixReference();
+
+ dcrsm->ClearRows(0, 0);
+ Finalize();
+ ExpectEqualToDenseReference();
+ ExpectEqualToCompressedRowSparseMatrixReference();
+
+ dcrsm->ClearRows(0, num_rows);
+ ExpectEqualToCompressedRowSparseMatrixReference();
+
+ Finalize();
+ ExpectEmpty();
+
+ InsertNonZeroEntriesFromDenseReference();
+ dcrsm->ClearRows(1, 2);
+ Finalize();
+ dense.block(1, 0, 2, num_cols).setZero();
+ ExpectEqualToDenseReference();
+
+ InitialiseDenseReference();
+}
+
+} // namespace internal
+} // namespace ceres
diff --git a/internal/ceres/evaluator.cc b/internal/ceres/evaluator.cc
index 31a4176..c94c62c 100644
--- a/internal/ceres/evaluator.cc
+++ b/internal/ceres/evaluator.cc
@@ -35,6 +35,8 @@
#include "ceres/compressed_row_sparse_matrix.h"
#include "ceres/crs_matrix.h"
#include "ceres/dense_jacobian_writer.h"
+#include "ceres/dynamic_compressed_row_finalizer.h"
+#include "ceres/dynamic_compressed_row_jacobian_writer.h"
#include "ceres/evaluator.h"
#include "ceres/internal/port.h"
#include "ceres/program_evaluator.h"
@@ -63,9 +65,17 @@
BlockJacobianWriter>(options,
program);
case SPARSE_NORMAL_CHOLESKY:
- return new ProgramEvaluator<ScratchEvaluatePreparer,
- CompressedRowJacobianWriter>(options,
- program);
+ if (options.dynamic_sparsity) {
+ return new ProgramEvaluator<ScratchEvaluatePreparer,
+ DynamicCompressedRowJacobianWriter,
+ DynamicCompressedRowJacobianFinalizer>(
+ options, program);
+ } else {
+ return new ProgramEvaluator<ScratchEvaluatePreparer,
+ CompressedRowJacobianWriter>(options,
+ program);
+ }
+
default:
*error = "Invalid Linear Solver Type. Unable to create evaluator.";
return NULL;
diff --git a/internal/ceres/evaluator.h b/internal/ceres/evaluator.h
index 3d25462..8fc60b8 100644
--- a/internal/ceres/evaluator.h
+++ b/internal/ceres/evaluator.h
@@ -61,11 +61,13 @@
Options()
: num_threads(1),
num_eliminate_blocks(-1),
- linear_solver_type(DENSE_QR) {}
+ linear_solver_type(DENSE_QR),
+ dynamic_sparsity(false) {}
int num_threads;
int num_eliminate_blocks;
LinearSolverType linear_solver_type;
+ bool dynamic_sparsity;
};
static Evaluator* Create(const Options& options,
diff --git a/internal/ceres/evaluator_test.cc b/internal/ceres/evaluator_test.cc
index ea24504..c0de3fc 100644
--- a/internal/ceres/evaluator_test.cc
+++ b/internal/ceres/evaluator_test.cc
@@ -44,6 +44,7 @@
#include "ceres/program.h"
#include "ceres/sized_cost_function.h"
#include "ceres/sparse_matrix.h"
+#include "ceres/stringprintf.h"
#include "ceres/types.h"
#include "gtest/gtest.h"
@@ -91,18 +92,42 @@
}
};
+struct EvaluatorTestOptions {
+ EvaluatorTestOptions(LinearSolverType linear_solver_type,
+ int num_eliminate_blocks,
+ bool dynamic_sparsity = false)
+ : linear_solver_type(linear_solver_type),
+ num_eliminate_blocks(num_eliminate_blocks),
+ dynamic_sparsity(dynamic_sparsity) {}
+
+ LinearSolverType linear_solver_type;
+ int num_eliminate_blocks;
+ bool dynamic_sparsity;
+};
+
struct EvaluatorTest
- : public ::testing::TestWithParam<pair<LinearSolverType, int> > {
+ : public ::testing::TestWithParam<EvaluatorTestOptions> {
Evaluator* CreateEvaluator(Program* program) {
// This program is straight from the ProblemImpl, and so has no index/offset
// yet; compute it here as required by the evalutor implementations.
program->SetParameterOffsetsAndIndex();
- VLOG(1) << "Creating evaluator with type: " << GetParam().first
- << " and num_eliminate_blocks: " << GetParam().second;
+ if (VLOG_IS_ON(1)) {
+ string report;
+ StringAppendF(&report, "Creating evaluator with type: %d",
+ GetParam().linear_solver_type);
+ if (GetParam().linear_solver_type == SPARSE_NORMAL_CHOLESKY) {
+ StringAppendF(&report, ", dynamic_sparsity: %d",
+ GetParam().dynamic_sparsity);
+ }
+ StringAppendF(&report, " and num_eliminate_blocks: %d",
+ GetParam().num_eliminate_blocks);
+ VLOG(1) << report;
+ }
Evaluator::Options options;
- options.linear_solver_type = GetParam().first;
- options.num_eliminate_blocks = GetParam().second;
+ options.linear_solver_type = GetParam().linear_solver_type;
+ options.num_eliminate_blocks = GetParam().num_eliminate_blocks;
+ options.dynamic_sparsity = GetParam().dynamic_sparsity;
string error;
return Evaluator::Create(options, program, &error);
}
@@ -517,23 +542,25 @@
INSTANTIATE_TEST_CASE_P(
LinearSolvers,
EvaluatorTest,
- ::testing::Values(make_pair(DENSE_QR, 0),
- make_pair(DENSE_SCHUR, 0),
- make_pair(DENSE_SCHUR, 1),
- make_pair(DENSE_SCHUR, 2),
- make_pair(DENSE_SCHUR, 3),
- make_pair(DENSE_SCHUR, 4),
- make_pair(SPARSE_SCHUR, 0),
- make_pair(SPARSE_SCHUR, 1),
- make_pair(SPARSE_SCHUR, 2),
- make_pair(SPARSE_SCHUR, 3),
- make_pair(SPARSE_SCHUR, 4),
- make_pair(ITERATIVE_SCHUR, 0),
- make_pair(ITERATIVE_SCHUR, 1),
- make_pair(ITERATIVE_SCHUR, 2),
- make_pair(ITERATIVE_SCHUR, 3),
- make_pair(ITERATIVE_SCHUR, 4),
- make_pair(SPARSE_NORMAL_CHOLESKY, 0)));
+ ::testing::Values(
+ EvaluatorTestOptions(DENSE_QR, 0),
+ EvaluatorTestOptions(DENSE_SCHUR, 0),
+ EvaluatorTestOptions(DENSE_SCHUR, 1),
+ EvaluatorTestOptions(DENSE_SCHUR, 2),
+ EvaluatorTestOptions(DENSE_SCHUR, 3),
+ EvaluatorTestOptions(DENSE_SCHUR, 4),
+ EvaluatorTestOptions(SPARSE_SCHUR, 0),
+ EvaluatorTestOptions(SPARSE_SCHUR, 1),
+ EvaluatorTestOptions(SPARSE_SCHUR, 2),
+ EvaluatorTestOptions(SPARSE_SCHUR, 3),
+ EvaluatorTestOptions(SPARSE_SCHUR, 4),
+ EvaluatorTestOptions(ITERATIVE_SCHUR, 0),
+ EvaluatorTestOptions(ITERATIVE_SCHUR, 1),
+ EvaluatorTestOptions(ITERATIVE_SCHUR, 2),
+ EvaluatorTestOptions(ITERATIVE_SCHUR, 3),
+ EvaluatorTestOptions(ITERATIVE_SCHUR, 4),
+ EvaluatorTestOptions(SPARSE_NORMAL_CHOLESKY, 0, false),
+ EvaluatorTestOptions(SPARSE_NORMAL_CHOLESKY, 0, true)));
// Simple cost function used to check if the evaluator is sensitive to
// state changes.
diff --git a/internal/ceres/linear_solver.h b/internal/ceres/linear_solver.h
index 8737c7c..f091bc5 100644
--- a/internal/ceres/linear_solver.h
+++ b/internal/ceres/linear_solver.h
@@ -98,6 +98,7 @@
dense_linear_algebra_library_type(EIGEN),
sparse_linear_algebra_library_type(SUITE_SPARSE),
use_postordering(false),
+ dynamic_sparsity(false),
min_num_iterations(1),
max_num_iterations(1),
num_threads(1),
@@ -115,6 +116,7 @@
// See solver.h for information about this flag.
bool use_postordering;
+ bool dynamic_sparsity;
// Number of internal iterations that the solver uses. This
// parameter only makes sense for iterative solvers like CG.
diff --git a/internal/ceres/program_evaluator.h b/internal/ceres/program_evaluator.h
index 8aa2a39..a06dc8c 100644
--- a/internal/ceres/program_evaluator.h
+++ b/internal/ceres/program_evaluator.h
@@ -97,7 +97,13 @@
namespace ceres {
namespace internal {
-template<typename EvaluatePreparer, typename JacobianWriter>
+struct NullJacobianFinalizer {
+ void operator()(SparseMatrix* jacobian, int num_parameters) {}
+};
+
+template<typename EvaluatePreparer,
+ typename JacobianWriter,
+ typename JacobianFinalizer = NullJacobianFinalizer>
class ProgramEvaluator : public Evaluator {
public:
ProgramEvaluator(const Evaluator::Options &options, Program* program)
@@ -244,9 +250,10 @@
}
if (!abort) {
+ const int num_parameters = program_->NumEffectiveParameters();
+
// Sum the cost and gradient (if requested) from each thread.
(*cost) = 0.0;
- int num_parameters = program_->NumEffectiveParameters();
if (gradient != NULL) {
VectorRef(gradient, num_parameters).setZero();
}
@@ -257,6 +264,15 @@
VectorRef(evaluate_scratch_[i].gradient.get(), num_parameters);
}
}
+
+ // Finalize the Jacobian if it is available.
+ // `num_parameters` is passed to the finalizer so that additional
+ // storage can be reserved for additional diagonal elements if
+ // necessary.
+ if (jacobian != NULL) {
+ JacobianFinalizer f;
+ f(jacobian, num_parameters);
+ }
}
return !abort;
}
diff --git a/internal/ceres/solver_impl.cc b/internal/ceres/solver_impl.cc
index e85ffbf..6ae2c90 100644
--- a/internal/ceres/solver_impl.cc
+++ b/internal/ceres/solver_impl.cc
@@ -1183,7 +1183,8 @@
return transformed_program.release();
}
- if (options->linear_solver_type == SPARSE_NORMAL_CHOLESKY) {
+ if (options->linear_solver_type == SPARSE_NORMAL_CHOLESKY &&
+ !options->dynamic_sparsity) {
if (!ReorderProgramForSparseNormalCholesky(
options->sparse_linear_algebra_library_type,
linear_solver_ordering,
@@ -1305,6 +1306,7 @@
linear_solver_options.dense_linear_algebra_library_type =
options->dense_linear_algebra_library_type;
linear_solver_options.use_postordering = options->use_postordering;
+ linear_solver_options.dynamic_sparsity = options->dynamic_sparsity;
// Ignore user's postordering preferences and force it to be true if
// cholmod_camd is not available. This ensures that the linear
@@ -1457,6 +1459,7 @@
->second.size())
: 0;
evaluator_options.num_threads = options.num_threads;
+ evaluator_options.dynamic_sparsity = options.dynamic_sparsity;
return Evaluator::Create(evaluator_options, program, error);
}
diff --git a/internal/ceres/sparse_normal_cholesky_solver.cc b/internal/ceres/sparse_normal_cholesky_solver.cc
index 880adf7..07537e3 100644
--- a/internal/ceres/sparse_normal_cholesky_solver.cc
+++ b/internal/ceres/sparse_normal_cholesky_solver.cc
@@ -56,13 +56,13 @@
options_(options) {
}
-SparseNormalCholeskySolver::~SparseNormalCholeskySolver() {
+void SparseNormalCholeskySolver::FreeFactorization() {
#ifndef CERES_NO_SUITESPARSE
if (factor_ != NULL) {
ss_.Free(factor_);
factor_ = NULL;
}
-#endif
+#endif // CERES_NO_SUITESPARSE
#ifndef CERES_NO_CXSPARSE
if (cxsparse_factor_ != NULL) {
@@ -72,6 +72,10 @@
#endif // CERES_NO_CXSPARSE
}
+SparseNormalCholeskySolver::~SparseNormalCholeskySolver() {
+ FreeFactorization();
+}
+
LinearSolver::Summary SparseNormalCholeskySolver::SolveImpl(
CompressedRowSparseMatrix* A,
const double* b,
@@ -150,13 +154,20 @@
event_logger.AddEvent("Setup");
// Compute symbolic factorization if not available.
+ if (options_.dynamic_sparsity) {
+ FreeFactorization();
+ }
if (cxsparse_factor_ == NULL) {
if (options_.use_postordering) {
cxsparse_factor_ = cxsparse_.BlockAnalyzeCholesky(AtA,
A->col_blocks(),
A->col_blocks());
} else {
- cxsparse_factor_ = cxsparse_.AnalyzeCholeskyWithNaturalOrdering(AtA);
+ if (options_.dynamic_sparsity) {
+ cxsparse_factor_ = cxsparse_.AnalyzeCholesky(AtA);
+ } else {
+ cxsparse_factor_ = cxsparse_.AnalyzeCholeskyWithNaturalOrdering(AtA);
+ }
}
}
event_logger.AddEvent("Analysis");
@@ -169,6 +180,7 @@
summary.termination_type = LINEAR_SOLVER_FAILURE;
}
event_logger.AddEvent("Solve");
+
return summary;
}
#else
@@ -198,6 +210,9 @@
cholmod_sparse lhs = ss_.CreateSparseMatrixTransposeView(A);
event_logger.AddEvent("Setup");
+ if (options_.dynamic_sparsity) {
+ FreeFactorization();
+ }
if (factor_ == NULL) {
if (options_.use_postordering) {
factor_ = ss_.BlockAnalyzeCholesky(&lhs,
@@ -205,7 +220,11 @@
A->row_blocks(),
&summary.message);
} else {
- factor_ = ss_.AnalyzeCholeskyWithNaturalOrdering(&lhs, &summary.message);
+ if (options_.dynamic_sparsity) {
+ factor_ = ss_.AnalyzeCholesky(&lhs, &summary.message);
+ } else {
+ factor_ = ss_.AnalyzeCholeskyWithNaturalOrdering(&lhs, &summary.message);
+ }
}
}
event_logger.AddEvent("Analysis");
diff --git a/internal/ceres/sparse_normal_cholesky_solver.h b/internal/ceres/sparse_normal_cholesky_solver.h
index adca08a..ea91534 100644
--- a/internal/ceres/sparse_normal_cholesky_solver.h
+++ b/internal/ceres/sparse_normal_cholesky_solver.h
@@ -71,6 +71,8 @@
const LinearSolver::PerSolveOptions& options,
double* rhs_and_solution);
+ void FreeFactorization();
+
SuiteSparse ss_;
// Cached factorization
cholmod_factor* factor_;
diff --git a/internal/ceres/unsymmetric_linear_solver_test.cc b/internal/ceres/unsymmetric_linear_solver_test.cc
index 624db68..412e611 100644
--- a/internal/ceres/unsymmetric_linear_solver_test.cc
+++ b/internal/ceres/unsymmetric_linear_solver_test.cc
@@ -172,6 +172,15 @@
options.use_postordering = true;
TestSolver(options);
}
+
+TEST_F(UnsymmetricLinearSolverTest,
+ SparseNormalCholeskyUsingSuiteSparseDynamicSparsity) {
+ LinearSolver::Options options;
+ options.sparse_linear_algebra_library_type = SUITE_SPARSE;
+ options.type = SPARSE_NORMAL_CHOLESKY;
+ options.dynamic_sparsity = true;
+ TestSolver(options);
+}
#endif
#ifndef CERES_NO_CXSPARSE
@@ -192,6 +201,15 @@
options.use_postordering = true;
TestSolver(options);
}
+
+TEST_F(UnsymmetricLinearSolverTest,
+ SparseNormalCholeskyUsingCXSparseDynamicSparsity) {
+ LinearSolver::Options options;
+ options.sparse_linear_algebra_library_type = CX_SPARSE;
+ options.type = SPARSE_NORMAL_CHOLESKY;
+ options.dynamic_sparsity = true;
+ TestSolver(options);
+}
#endif
} // namespace internal
