An implementation of SubsetPreconditioner.
The key idea being, use some subset of the rows of the Jacobian
as the preconditioner.
This CL only implements the preconditioner assuming that the row
selection has already been done. How the rows are selected will be
left to the user based on their knowledge of the problem.
A follow up CL will hook this preconditioner into the rest of the
solver.
Change-Id: I3e18dc57811116534e9ddf35d7b154bcce496d3b
diff --git a/BUILD b/BUILD
index a82d2cf..504a23e 100644
--- a/BUILD
+++ b/BUILD
@@ -141,6 +141,7 @@
"solver",
"sparse_cholesky",
"sparse_normal_cholesky_solver",
+ "subset_preconditioner",
"system",
"tiny_solver_autodiff_function",
"tiny_solver_cost_function_adapter",
diff --git a/bazel/ceres.bzl b/bazel/ceres.bzl
index 833c937..d9d7502 100644
--- a/bazel/ceres.bzl
+++ b/bazel/ceres.bzl
@@ -111,6 +111,7 @@
"sparse_normal_cholesky_solver.cc",
"split.cc",
"stringprintf.cc",
+ "subset_preconditioner.cc",
"suitesparse.cc",
"thread_token_provider.cc",
"triplet_sparse_matrix.cc",
diff --git a/internal/ceres/CMakeLists.txt b/internal/ceres/CMakeLists.txt
index 884fd57..91ffc11 100644
--- a/internal/ceres/CMakeLists.txt
+++ b/internal/ceres/CMakeLists.txt
@@ -111,6 +111,7 @@
sparse_matrix.cc
sparse_cholesky.cc
sparse_normal_cholesky_solver.cc
+ subset_preconditioner.cc
split.cc
stringprintf.cc
suitesparse.cc
@@ -359,6 +360,7 @@
ceres_test(solver)
ceres_test(sparse_cholesky)
ceres_test(sparse_normal_cholesky_solver)
+ ceres_test(subset_preconditioner)
ceres_test(system)
ceres_test(tiny_solver)
ceres_test(tiny_solver_autodiff_function)
diff --git a/internal/ceres/block_sparse_matrix.cc b/internal/ceres/block_sparse_matrix.cc
index 214daef..c814a72 100644
--- a/internal/ceres/block_sparse_matrix.cc
+++ b/internal/ceres/block_sparse_matrix.cc
@@ -280,9 +280,12 @@
}
void BlockSparseMatrix::AppendRows(const BlockSparseMatrix& m) {
+ CHECK_EQ(m.num_cols(), num_cols());
+ const CompressedRowBlockStructure* m_bs = m.block_structure();
+ CHECK_EQ(m_bs->cols.size(), block_structure_->cols.size());
+
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) {
@@ -336,29 +339,33 @@
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;
- }
+ if (options.col_blocks.empty()) {
+ 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);
+ // 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;
+ }
+ } else {
+ bs->cols = options.col_blocks;
+ }
bool matrix_has_blocks = false;
while (!matrix_has_blocks) {
- LOG(INFO) << "clearing";
+ VLOG(1) << "Clearing";
bs->rows.clear();
int row_block_position = 0;
int value_position = 0;
@@ -372,7 +379,7 @@
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) {
+ for (int c = 0; c < bs->cols.size(); ++c) {
if (RandDouble() > options.block_density) continue;
row.cells.push_back(Cell());
diff --git a/internal/ceres/block_sparse_matrix.h b/internal/ceres/block_sparse_matrix.h
index 17491bf..abde2a6 100644
--- a/internal/ceres/block_sparse_matrix.h
+++ b/internal/ceres/block_sparse_matrix.h
@@ -117,6 +117,11 @@
// present in the matrix. A given random matrix will not have
// precisely this density.
double block_density;
+
+ // If col_blocks is non-empty, then the generated random matrix
+ // has this block structure and the column related options in this
+ // struct are ignored.
+ std::vector<Block> col_blocks;
};
// Create a random BlockSparseMatrix whose entries are normally
diff --git a/internal/ceres/preconditioner.h b/internal/ceres/preconditioner.h
index a248eae..5d293b8 100644
--- a/internal/ceres/preconditioner.h
+++ b/internal/ceres/preconditioner.h
@@ -51,6 +51,8 @@
: type(JACOBI),
visibility_clustering_type(CANONICAL_VIEWS),
sparse_linear_algebra_library_type(SUITE_SPARSE),
+ use_postordering(false),
+ subset_preconditioner_start_row_block(-1),
num_threads(1),
row_block_size(Eigen::Dynamic),
e_block_size(Eigen::Dynamic),
@@ -61,6 +63,21 @@
VisibilityClusteringType visibility_clustering_type;
SparseLinearAlgebraLibraryType sparse_linear_algebra_library_type;
+ // When using the subset preconditioner, all row blocks starting
+ // from this row block are used to construct the preconditioner.
+ //
+ // i.e., the Jacobian matrix A is horizonatally partitioned as
+ //
+ // A = [P]
+ // [Q]
+ //
+ // where P has subset_preconditioner_start_row_block row blocks,
+ // and the preconditioner is the inverse of the matrix Q'Q.
+ int subset_preconditioner_start_row_block;
+
+ // See solver.h for information about these flags.
+ bool use_postordering;
+
// If possible, how many threads the preconditioner can use.
int num_threads;
diff --git a/internal/ceres/subset_preconditioner.cc b/internal/ceres/subset_preconditioner.cc
new file mode 100644
index 0000000..e970b91
--- /dev/null
+++ b/internal/ceres/subset_preconditioner.cc
@@ -0,0 +1,112 @@
+// Ceres Solver - A fast non-linear least squares minimizer
+// Copyright 2017 Google Inc. All rights reserved.
+// http://ceres-solver.org/
+//
+// 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: sameeragarwal@google.com (Sameer Agarwal)
+
+#include "ceres/subset_preconditioner.h"
+
+#include <string>
+#include "ceres/compressed_row_sparse_matrix.h"
+#include "ceres/inner_product_computer.h"
+#include "ceres/internal/scoped_ptr.h"
+#include "ceres/linear_solver.h"
+#include "ceres/sparse_cholesky.h"
+#include "ceres/types.h"
+
+namespace ceres {
+namespace internal {
+
+SubsetPreconditioner::SubsetPreconditioner(
+ const Preconditioner::Options& options, const BlockSparseMatrix& A)
+ : options_(options), num_cols_(A.num_cols()) {
+ sparse_cholesky_.reset(
+ SparseCholesky::Create(options_.sparse_linear_algebra_library_type,
+ options_.use_postordering ? AMD : NATURAL));
+ CHECK_GE(options_.subset_preconditioner_start_row_block, 0);
+}
+
+SubsetPreconditioner::~SubsetPreconditioner() {}
+
+void SubsetPreconditioner::RightMultiply(const double* x, double* y) const {
+ CHECK_NOTNULL(x);
+ CHECK_NOTNULL(y);
+ std::string message;
+ sparse_cholesky_->Solve(x, y, &message);
+}
+
+bool SubsetPreconditioner::UpdateImpl(const BlockSparseMatrix& A,
+ const double* D) {
+ BlockSparseMatrix* m = const_cast<BlockSparseMatrix*>(&A);
+ const CompressedRowBlockStructure* bs = m->block_structure();
+
+ // A = [P]
+ // [Q]
+
+ // Now add D to A if needed.
+ if (D != NULL) {
+ // A = [P]
+ // [Q]
+ // [D]
+ scoped_ptr<BlockSparseMatrix> regularizer(
+ BlockSparseMatrix::CreateDiagonalMatrix(D, bs->cols));
+ m->AppendRows(*regularizer);
+ }
+
+ if (inner_product_computer_.get() == NULL) {
+ inner_product_computer_.reset(InnerProductComputer::Create(
+ *m,
+ options_.subset_preconditioner_start_row_block,
+ bs->rows.size(),
+ sparse_cholesky_->StorageType()));
+ }
+
+ // Compute inner_product = [Q'*Q + D'*D]
+ inner_product_computer_->Compute();
+
+ // Unappend D if needed.
+ if (D != NULL) {
+ // A = [P]
+ // [Q]
+ m->DeleteRowBlocks(bs->cols.size());
+ }
+
+ std::string message;
+ // Compute L. s.t., LL' = Q'*Q + D'*D
+ const LinearSolverTerminationType termination_type =
+ sparse_cholesky_->Factorize(inner_product_computer_->mutable_result(),
+ &message);
+ if (termination_type != LINEAR_SOLVER_SUCCESS) {
+ LOG(ERROR) << "Preconditioner factorization failed: " << message;
+ return false;
+ }
+
+ return true;
+}
+
+} // namespace internal
+} // namespace ceres
diff --git a/internal/ceres/subset_preconditioner.h b/internal/ceres/subset_preconditioner.h
new file mode 100644
index 0000000..062253b
--- /dev/null
+++ b/internal/ceres/subset_preconditioner.h
@@ -0,0 +1,91 @@
+// Ceres Solver - A fast non-linear least squares minimizer
+// Copyright 2017 Google Inc. All rights reserved.
+// http://ceres-solver.org/
+//
+// 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: sameeragarwal@google.com (Sameer Agarwal)
+
+#ifndef CERES_INTERNAL_SUBSET_PRECONDITIONER_H_
+#define CERES_INTERNAL_SUBSET_PRECONDITIONER_H_
+
+#include "ceres/internal/scoped_ptr.h"
+#include "ceres/preconditioner.h"
+
+namespace ceres {
+namespace internal {
+
+class BlockSparseMatrix;
+class SparseCholesky;
+class InnerProductComputer;
+
+// Subset preconditioning, uses a subset of the rows of the Jacobian
+// to construct a preconditioner for the normal equations.
+//
+// To keep the interface simple, we assume that the matrix A has
+// already been re-ordered that the user wishes to some subset of the
+// bottom row blocks of the matrix as the preconditioner. This is
+// controlled by
+// Preconditioner::Options::subset_preconditioner_start_row_block.
+//
+// When using the subset preconditioner, all row blocks starting
+// from this row block are used to construct the preconditioner.
+//
+// More precisely the matrix A is horizontally partitioned as
+//
+// A = [P]
+// [Q]
+//
+// where P as subset_preconditioner_start_row_block row blocks,
+// and the preconditioner is the inverse of the matrix Q'Q.
+//
+// Obviously, the smaller this number, the more accurate and
+// computationally expensive this preconditioner will be.
+//
+// See the tests for example usage.
+class SubsetPreconditioner : public BlockSparseMatrixPreconditioner {
+ public:
+ SubsetPreconditioner(const Preconditioner::Options& options,
+ const BlockSparseMatrix& A);
+ virtual ~SubsetPreconditioner();
+
+ // Preconditioner interface
+ virtual void RightMultiply(const double* x, double* y) const;
+ virtual int num_rows() const { return num_cols_; }
+ virtual int num_cols() const { return num_cols_; }
+
+ private:
+ virtual bool UpdateImpl(const BlockSparseMatrix& A, const double* D);
+
+ const Preconditioner::Options options_;
+ const int num_cols_;
+ scoped_ptr<SparseCholesky> sparse_cholesky_;
+ scoped_ptr<InnerProductComputer> inner_product_computer_;
+};
+
+} // namespace internal
+} // namespace ceres
+
+#endif // CERES_INTERNAL_SUBSET_PRECONDITIONER_H_
diff --git a/internal/ceres/subset_preconditioner_test.cc b/internal/ceres/subset_preconditioner_test.cc
new file mode 100644
index 0000000..63dbbaf
--- /dev/null
+++ b/internal/ceres/subset_preconditioner_test.cc
@@ -0,0 +1,194 @@
+// Ceres Solver - A fast non-linear least squares minimizer
+// Copyright 2017 Google Inc. All rights reserved.
+// http://ceres-solver.org/
+//
+// 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: sameeragarwal@google.com (Sameer Agarwal)
+
+#include "ceres/subset_preconditioner.h"
+#include "Eigen/Dense"
+#include "Eigen/SparseCore"
+#include "ceres/block_sparse_matrix.h"
+#include "ceres/compressed_row_sparse_matrix.h"
+#include "ceres/inner_product_computer.h"
+#include "ceres/internal/eigen.h"
+#include "ceres/internal/scoped_ptr.h"
+#include "glog/logging.h"
+#include "gtest/gtest.h"
+
+namespace ceres {
+namespace internal {
+
+// TODO(sameeragarwal): Refactor the following two functions out of
+// here and sparse_cholesky_test.cc into a more suitable place.
+template <Eigen::UpLoType UpLoType>
+bool SolveLinearSystemUsingEigen(const Matrix& lhs,
+ const Vector rhs,
+ Vector* solution) {
+ Eigen::LLT<Matrix, UpLoType> llt = lhs.selfadjointView<UpLoType>().llt();
+ if (llt.info() != Eigen::Success) {
+ return false;
+ }
+ *solution = llt.solve(rhs);
+ return (llt.info() == Eigen::Success);
+}
+
+// Use Eigen's Dense Cholesky solver to compute the solution to a
+// sparse linear system.
+bool ComputeExpectedSolution(const CompressedRowSparseMatrix& lhs,
+ const Vector& rhs,
+ Vector* solution) {
+ Matrix dense_triangular_lhs;
+ lhs.ToDenseMatrix(&dense_triangular_lhs);
+ if (lhs.storage_type() == CompressedRowSparseMatrix::UPPER_TRIANGULAR) {
+ Matrix full_lhs = dense_triangular_lhs.selfadjointView<Eigen::Upper>();
+ return SolveLinearSystemUsingEigen<Eigen::Upper>(full_lhs, rhs, solution);
+ }
+ return SolveLinearSystemUsingEigen<Eigen::Lower>(
+ dense_triangular_lhs, rhs, solution);
+}
+
+typedef ::testing::tuple<SparseLinearAlgebraLibraryType, bool> Param;
+
+std::string ParamInfoToString(testing::TestParamInfo<Param> info) {
+ Param param = info.param;
+ std::stringstream ss;
+ ss << SparseLinearAlgebraLibraryTypeToString(::testing::get<0>(param)) << "_"
+ << (::testing::get<1>(param) ? "Diagonal" : "NoDiagonal");
+ return ss.str();
+}
+
+class SubsetPreconditionerTest : public ::testing::TestWithParam<Param> {
+ protected:
+ virtual void SetUp() {
+ BlockSparseMatrix::RandomMatrixOptions options;
+ options.num_col_blocks = 4;
+ options.min_col_block_size = 1;
+ options.max_col_block_size = 4;
+ options.num_row_blocks = 8;
+ options.min_row_block_size = 1;
+ options.max_row_block_size = 4;
+ options.block_density = 0.9;
+
+ m_.reset(BlockSparseMatrix::CreateRandomMatrix(options));
+ start_row_block_ = m_->block_structure()->rows.size();
+
+ // Ensure that the bottom part of the matrix has the same column
+ // block structure.
+ options.col_blocks = m_->block_structure()->cols;
+ b_.reset(BlockSparseMatrix::CreateRandomMatrix(options));
+ m_->AppendRows(*b_);
+
+ // Create a Identity block diagonal matrix with the same column
+ // block structure.
+ diagonal_ = Vector::Ones(m_->num_cols());
+ block_diagonal_.reset(BlockSparseMatrix::CreateDiagonalMatrix(
+ diagonal_.data(), b_->block_structure()->cols));
+
+ // Unconditionally add the block diagonal to the matrix b_,
+ // because either it is either part of b_ to make it full rank, or
+ // we pass the same diagonal matrix later as the parameter D. In
+ // either case the preconditioner matrix is b_' b + D'D.
+ b_->AppendRows(*block_diagonal_);
+ inner_product_computer_.reset(InnerProductComputer::Create(
+ *b_, CompressedRowSparseMatrix::UPPER_TRIANGULAR));
+ inner_product_computer_->Compute();
+ }
+
+ scoped_ptr<BlockSparseMatrix> m_;
+ scoped_ptr<BlockSparseMatrix> b_;
+ scoped_ptr<BlockSparseMatrix> block_diagonal_;
+ scoped_ptr<InnerProductComputer> inner_product_computer_;
+ scoped_ptr<Preconditioner> preconditioner_;
+ Vector diagonal_;
+ int start_row_block_;
+};
+
+TEST_P(SubsetPreconditionerTest, foo) {
+ Param param = GetParam();
+ Preconditioner::Options options;
+ options.subset_preconditioner_start_row_block = start_row_block_;
+ options.sparse_linear_algebra_library_type = ::testing::get<0>(param);
+ preconditioner_.reset(new SubsetPreconditioner(options, *m_));
+
+ const bool with_diagonal = ::testing::get<1>(param);
+ if (!with_diagonal) {
+ m_->AppendRows(*block_diagonal_);
+ }
+
+ EXPECT_TRUE(
+ preconditioner_->Update(*m_, with_diagonal ? diagonal_.data() : NULL));
+
+ // Repeatedly apply the preconditioner to random vectors and check
+ // that the preconditioned value is the same as one obtained by
+ // solving the linear system directly.
+ for (int i = 0; i < 5; ++i) {
+ CompressedRowSparseMatrix* lhs = inner_product_computer_->mutable_result();
+ Vector rhs = Vector::Random(lhs->num_rows());
+ Vector expected(lhs->num_rows());
+ EXPECT_TRUE(ComputeExpectedSolution(*lhs, rhs, &expected));
+
+ Vector actual(lhs->num_rows());
+ preconditioner_->RightMultiply(rhs.data(), actual.data());
+
+ Matrix eigen_lhs;
+ lhs->ToDenseMatrix(&eigen_lhs);
+ EXPECT_NEAR((actual - expected).norm() / actual.norm(),
+ 0.0,
+ std::numeric_limits<double>::epsilon() * 10)
+ << "\n"
+ << eigen_lhs << "\n"
+ << expected.transpose() << "\n"
+ << actual.transpose();
+ }
+}
+
+#ifndef CERES_NO_SUITESPARSE
+INSTANTIATE_TEST_CASE_P(SubsetPreconditionerWithSuiteSparse,
+ SubsetPreconditionerTest,
+ ::testing::Combine(::testing::Values(SUITE_SPARSE),
+ ::testing::Values(true, false)),
+ ParamInfoToString);
+#endif
+
+#ifndef CERES_NO_CXSPARSE
+INSTANTIATE_TEST_CASE_P(SubsetPreconditionerWithCXSparse,
+ SubsetPreconditionerTest,
+ ::testing::Combine(::testing::Values(CX_SPARSE),
+ ::testing::Values(true, false)),
+ ParamInfoToString);
+#endif
+
+#ifndef CERES_NO_EIGEN_SPARSE
+INSTANTIATE_TEST_CASE_P(SubsetPreconditionerWithEigenSparse,
+ SubsetPreconditionerTest,
+ ::testing::Combine(::testing::Values(EIGEN_SPARSE),
+ ::testing::Values(true, false)),
+ ParamInfoToString);
+#endif
+
+} // namespace internal
+} // namespace ceres
diff --git a/jni/Android.mk b/jni/Android.mk
index 8dbb593..4ddcdf8 100644
--- a/jni/Android.mk
+++ b/jni/Android.mk
@@ -195,6 +195,7 @@
$(CERES_SRC_PATH)/sparse_normal_cholesky_solver.cc \
$(CERES_SRC_PATH)/split.cc \
$(CERES_SRC_PATH)/stringprintf.cc \
+ $(CERES_SRC_PATH)/subset_preconditioner.cc \
$(CERES_SRC_PATH)/suitesparse.cc \
$(CERES_SRC_PATH)/triplet_sparse_matrix.cc \
$(CERES_SRC_PATH)/trust_region_minimizer.cc \
