Add a general sparse iterative solver: CGNR
This adds a new LinearOperator which implements symmetric
products of a matrix, and a new CGNR solver to leverage
CG to directly solve the normal equations. This also
includes a block diagonal preconditioner. In experiments
on problem-16, the non-preconditioned version is about
1/5 the speed of SPARSE_SCHUR, and the preconditioned
version using block cholesky is about 20% slower than
SPARSE_SCHUR.
diff --git a/internal/ceres/CMakeLists.txt b/internal/ceres/CMakeLists.txt
index 7891138..03d650a 100644
--- a/internal/ceres/CMakeLists.txt
+++ b/internal/ceres/CMakeLists.txt
@@ -29,6 +29,7 @@
# Author: keir@google.com (Keir Mierle)
SET(CERES_INTERNAL_SRC
+ block_diagonal_preconditioner.cc
block_evaluate_preparer.cc
block_jacobian_writer.cc
block_random_access_dense_matrix.cc
@@ -37,6 +38,7 @@
block_sparse_matrix.cc
block_structure.cc
canonical_views_clustering.cc
+ cgnr_solver.cc
compressed_row_jacobian_writer.cc
compressed_row_sparse_matrix.cc
conditioned_cost_function.cc
diff --git a/internal/ceres/block_diagonal_preconditioner.cc b/internal/ceres/block_diagonal_preconditioner.cc
new file mode 100644
index 0000000..0779a91
--- /dev/null
+++ b/internal/ceres/block_diagonal_preconditioner.cc
@@ -0,0 +1,118 @@
+// Ceres Solver - A fast non-linear least squares minimizer
+// Copyright 2012 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: keir@google.com (Keir Mierle)
+
+#include "ceres/block_diagonal_preconditioner.h"
+
+#include "Eigen/Cholesky"
+#include "ceres/block_sparse_matrix.h"
+#include "ceres/block_structure.h"
+#include "ceres/casts.h"
+#include "ceres/internal/eigen.h"
+
+namespace ceres {
+namespace internal {
+
+BlockDiagonalPreconditioner::BlockDiagonalPreconditioner(
+ const LinearOperator& A)
+ : block_structure_(
+ *(down_cast<const BlockSparseMatrix*>(&A)->block_structure())) {
+
+ // Calculate the amount of storage needed.
+ int storage_needed = 0;
+ for (int c = 0; c < block_structure_.cols.size(); ++c) {
+ int size = block_structure_.cols[c].size;
+ storage_needed += size * size;
+ }
+
+ // Size the offsets and storage.
+ blocks_.resize(block_structure_.cols.size());
+ block_storage_.resize(storage_needed);
+
+ // Put pointers to the storage in the offsets.
+ double *block_cursor = &block_storage_[0];
+ for (int c = 0; c < block_structure_.cols.size(); ++c) {
+ int size = block_structure_.cols[c].size;
+ blocks_[c] = block_cursor;
+ block_cursor += size * size;
+ }
+}
+
+BlockDiagonalPreconditioner::~BlockDiagonalPreconditioner() {
+}
+
+void BlockDiagonalPreconditioner::Update(const LinearOperator& matrix) {
+ const BlockSparseMatrix& A = *(down_cast<const BlockSparseMatrix*>(&matrix));
+ const CompressedRowBlockStructure* bs = A.block_structure();
+
+ // Compute the diagonal blocks by block inner products.
+ std::fill(block_storage_.begin(), block_storage_.end(), 0.0);
+ for (int r = 0; r < bs->rows.size(); ++r) {
+ const int row_block_size = bs->rows[r].block.size;
+ const vector<Cell>& cells = bs->rows[r].cells;
+ const double* row_values = A.RowBlockValues(r);
+ for (int c = 0; c < cells.size(); ++c) {
+ const int col_block_size = bs->cols[cells[c].block_id].size;
+ ConstMatrixRef m(row_values + cells[c].position,
+ row_block_size,
+ col_block_size);
+
+ MatrixRef(blocks_[cells[c].block_id],
+ col_block_size,
+ col_block_size).noalias() += m.transpose() * m;
+ }
+ }
+
+ // Invert each block.
+ for (int c = 0; c < bs->cols.size(); ++c) {
+ const int size = block_structure_.cols[c].size;
+ MatrixRef D(blocks_[c], size, size);
+ D = D.selfadjointView<Eigen::Upper>()
+ .ldlt()
+ .solve(Matrix::Identity(size, size));
+ }
+}
+
+void BlockDiagonalPreconditioner::RightMultiply(const double* x, double* y) const {
+ for (int c = 0; c < block_structure_.cols.size(); ++c) {
+ const int size = block_structure_.cols[c].size;
+ const int position = block_structure_.cols[c].position;
+ ConstMatrixRef D(blocks_[c], size, size);
+ ConstVectorRef x_block(x + position, size);
+ VectorRef y_block(y + position, size);
+ y_block += D * x_block;
+ }
+}
+
+void BlockDiagonalPreconditioner::LeftMultiply(const double* x, double* y) const {
+ RightMultiply(x, y);
+}
+
+} // namespace internal
+} // namespace ceres
diff --git a/internal/ceres/block_diagonal_preconditioner.h b/internal/ceres/block_diagonal_preconditioner.h
new file mode 100644
index 0000000..561b0db
--- /dev/null
+++ b/internal/ceres/block_diagonal_preconditioner.h
@@ -0,0 +1,68 @@
+// Ceres Solver - A fast non-linear least squares minimizer
+// Copyright 2012 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: keir@google.com (Keir Mierle)
+
+#ifndef CERES_INTERNAL_BLOCK_DIAGONAL_PRECONDITIONER_H_
+#define CERES_INTERNAL_BLOCK_DIAGONAL_PRECONDITIONER_H_
+
+#include <vector>
+#include "ceres/linear_operator.h"
+
+namespace ceres {
+namespace internal {
+
+class CompressedRowBlockStructure;
+class LinearOperator;
+class SparseMatrix;
+
+// A block diagonal preconditioner; also known as block-Jacobi.
+class BlockDiagonalPreconditioner : public LinearOperator {
+ public:
+ BlockDiagonalPreconditioner(const LinearOperator& A);
+ virtual ~BlockDiagonalPreconditioner();
+
+ void Update(const LinearOperator& matrix);
+
+ virtual void RightMultiply(const double* x, double* y) const;
+ virtual void LeftMultiply(const double* x, double* y) const;
+
+ virtual int num_rows() const { return size_; }
+ virtual int num_cols() const { return size_; }
+
+ private:
+ std::vector<double*> blocks_;
+ std::vector<double> block_storage_;
+ int size_;
+ const CompressedRowBlockStructure& block_structure_;
+};
+
+} // namespace internal
+} // namespace ceres
+
+#endif // CERES_INTERNAL_LINEAR_OPERATOR_H_
diff --git a/internal/ceres/cgnr_linear_operator.h b/internal/ceres/cgnr_linear_operator.h
new file mode 100644
index 0000000..94767fb
--- /dev/null
+++ b/internal/ceres/cgnr_linear_operator.h
@@ -0,0 +1,122 @@
+// Ceres Solver - A fast non-linear least squares minimizer
+// Copyright 2012 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: keir@google.com (Keir Mierle)
+
+#ifndef CERES_INTERNAL_CGNR_LINEAR_OPERATOR_H_
+#define CERES_INTERNAL_CGNR_LINEAR_OPERATOR_H_
+
+#include <algorithm>
+#include "ceres/linear_operator.h"
+#include "ceres/internal/scoped_ptr.h"
+#include "ceres/internal/eigen.h"
+
+namespace ceres {
+namespace internal {
+
+class SparseMatrix;
+
+// A linear operator which takes a matrix A and a diagonal vector D and
+// performs products of the form
+//
+// (A^T A + D^T D)x
+//
+// This is used to implement iterative general sparse linear solving with
+// conjugate gradients, where A is the Jacobian and D is a regularizing
+// parameter. A brief proof that D^T D is the correct regularizer:
+//
+// Given a regularized least squares problem:
+//
+// min ||Ax - b||^2 + ||Dx||^2
+// x
+//
+// First expand into matrix notation:
+//
+// (Ax - b)^T (Ax - b) + xD^TDx
+//
+// Then multiply out to get:
+//
+// = xA^TAx - 2xA^T b + b^Tb + xD^TDx
+// = xA^TAx - 2b^T Ax + b^Tb + xD^TDx
+//
+// Take the derivative:
+//
+// 0 = 2A^TAx - 2b^T A + b^Tb + 2 D^TDx
+// 0 = A^TAx - A^T b + D^TDx
+// 0 = (A^TA + D^TD)x - A^T b
+//
+// Thus, the symmetric system we need to solve for CGNR is
+//
+// Sx = z
+//
+// with S = A^TA + D^TD
+// and z = A^T b
+//
+// Note: This class is not thread safe, since it uses some temporary storage.
+class CgnrLinearOperator : public LinearOperator {
+ public:
+ CgnrLinearOperator(LinearOperator* A, double *D)
+ : A_(A), D_(D), z_(new double[A->num_rows()]) {
+ }
+ virtual ~CgnrLinearOperator() {}
+
+ virtual void RightMultiply(const double* x, double* y) const {
+ std::fill(z_.get(), z_.get() + A_->num_rows(), 0.0);
+
+ // z = Ax
+ A_->RightMultiply(x, z_.get());
+
+ // y = y + Atz
+ A_->LeftMultiply(z_.get(), y);
+
+ // y = y + DtDx
+ if (D_ != NULL) {
+ int n = A_->num_cols();
+ for (int i = 0; i < n; ++i) {
+ y[i] += D_[i] * D_[i] * x[i];
+ }
+ }
+ }
+
+ virtual void LeftMultiply(const double* x, double* y) const {
+ RightMultiply(x, y);
+ }
+
+ virtual int num_rows() const { return A_->num_cols(); }
+ virtual int num_cols() const { return A_->num_cols(); }
+
+ private:
+ LinearOperator* A_;
+ double* D_;
+ scoped_array<double> z_;
+};
+
+} // namespace internal
+} // namespace ceres
+
+#endif // CERES_INTERNAL_CGNR_LINEAR_OPERATOR_H_
diff --git a/internal/ceres/cgnr_solver.cc b/internal/ceres/cgnr_solver.cc
new file mode 100644
index 0000000..2c88e33
--- /dev/null
+++ b/internal/ceres/cgnr_solver.cc
@@ -0,0 +1,79 @@
+// Ceres Solver - A fast non-linear least squares minimizer
+// Copyright 2012 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: keir@google.com (Keir Mierle)
+
+#include "ceres/cgnr_solver.h"
+
+#include "ceres/linear_solver.h"
+#include "ceres/cgnr_linear_operator.h"
+#include "ceres/conjugate_gradients_solver.h"
+#include "ceres/block_diagonal_preconditioner.h"
+
+namespace ceres {
+namespace internal {
+
+CgnrSolver::CgnrSolver(const LinearSolver::Options& options)
+ : options_(options),
+ jacobi_preconditioner_(NULL) {
+}
+
+LinearSolver::Summary CgnrSolver::Solve(
+ LinearOperator* A,
+ const double* b,
+ const LinearSolver::PerSolveOptions& per_solve_options,
+ double* x) {
+ // Form z = Atb.
+ scoped_array<double> z(new double[A->num_cols()]);
+ std::fill(z.get(), z.get() + A->num_cols(), 0.0);
+ A->LeftMultiply(b, z.get());
+
+ // Precondition if necessary.
+ LinearSolver::PerSolveOptions cg_per_solve_options = per_solve_options;
+ if (options_.preconditioner_type == JACOBI) {
+ if (jacobi_preconditioner_.get() == NULL) {
+ jacobi_preconditioner_.reset(new BlockDiagonalPreconditioner(*A));
+ }
+ jacobi_preconditioner_->Update(*A);
+ cg_per_solve_options.preconditioner = jacobi_preconditioner_.get();
+ } else if (options_.preconditioner_type != IDENTITY) {
+ // TODO(keir): This should die somehow.
+ }
+
+ // Solve (AtA + DtD)x = z (= Atb).
+ std::fill(x, x + A->num_cols(), 0.0);
+ CgnrLinearOperator AtApDtD(A, per_solve_options.D);
+ ConjugateGradientsSolver conjugate_gradient_solver(options_);
+ return conjugate_gradient_solver.Solve(&AtApDtD,
+ z.get(),
+ cg_per_solve_options,
+ x);
+}
+
+} // namespace internal
+} // namespace ceres
diff --git a/internal/ceres/cgnr_solver.h b/internal/ceres/cgnr_solver.h
new file mode 100644
index 0000000..db49b07
--- /dev/null
+++ b/internal/ceres/cgnr_solver.h
@@ -0,0 +1,66 @@
+// Ceres Solver - A fast non-linear least squares minimizer
+// Copyright 2012 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: keir@google.com (Keir Mierle)
+
+#ifndef CERES_INTERNAL_CGNR_SOLVER_H_
+#define CERES_INTERNAL_CGNR_SOLVER_H_
+
+#include "ceres/internal/scoped_ptr.h"
+#include "ceres/linear_solver.h"
+
+namespace ceres {
+namespace internal {
+
+class BlockDiagonalPreconditioner;
+
+// A conjugate gradients on the normal equations solver. This directly solves
+// for the solution to
+//
+// (A^T A + D^T D)x = A^T b
+//
+// as required for solving for x in the least squares sense. Currently only
+// block diagonal preconditioning is supported.
+class CgnrSolver : public LinearSolver {
+ public:
+ explicit CgnrSolver(const LinearSolver::Options& options);
+ virtual Summary Solve(LinearOperator* A,
+ const double* b,
+ const LinearSolver::PerSolveOptions& per_solve_options,
+ double* x);
+
+ private:
+ const LinearSolver::Options options_;
+ scoped_ptr<BlockDiagonalPreconditioner> jacobi_preconditioner_;
+ DISALLOW_COPY_AND_ASSIGN(CgnrSolver);
+};
+
+} // namespace internal
+} // namespace ceres
+
+#endif // CERES_INTERNAL_CGNR_SOLVER_H_
diff --git a/internal/ceres/evaluator.cc b/internal/ceres/evaluator.cc
index 38c3fdb..7ceb7ac 100644
--- a/internal/ceres/evaluator.cc
+++ b/internal/ceres/evaluator.cc
@@ -53,6 +53,7 @@
case DENSE_SCHUR:
case SPARSE_SCHUR:
case ITERATIVE_SCHUR:
+ case CONJUGATE_GRADIENTS:
return new ProgramEvaluator<BlockEvaluatePreparer,
BlockJacobianWriter>(options,
program);
@@ -60,11 +61,6 @@
return new ProgramEvaluator<ScratchEvaluatePreparer,
CompressedRowJacobianWriter>(options,
program);
- case CONJUGATE_GRADIENTS:
- *error = "CONJUGATE_GRADIENTS is not supported as a linear "
- "solver. If you want to use an iterative linear solver, "
- "use ITERATIVE_SCHUR instead.";
- return NULL;
default:
*error = "Invalid Linear Solver Type. Unable to create evaluator.";
return NULL;
diff --git a/internal/ceres/iterative_schur_complement_solver.cc b/internal/ceres/iterative_schur_complement_solver.cc
index 48d8453..4379ebd 100644
--- a/internal/ceres/iterative_schur_complement_solver.cc
+++ b/internal/ceres/iterative_schur_complement_solver.cc
@@ -38,13 +38,14 @@
#include "Eigen/Dense"
#include "ceres/block_sparse_matrix.h"
#include "ceres/block_structure.h"
+#include "ceres/conjugate_gradients_solver.h"
#include "ceres/implicit_schur_complement.h"
-#include "ceres/linear_solver.h"
-#include "ceres/triplet_sparse_matrix.h"
-#include "ceres/visibility_based_preconditioner.h"
#include "ceres/internal/eigen.h"
#include "ceres/internal/scoped_ptr.h"
+#include "ceres/linear_solver.h"
+#include "ceres/triplet_sparse_matrix.h"
#include "ceres/types.h"
+#include "ceres/visibility_based_preconditioner.h"
namespace ceres {
namespace internal {
@@ -85,7 +86,7 @@
LinearSolver::Options cg_options;
cg_options.type = CONJUGATE_GRADIENTS;
cg_options.max_num_iterations = options_.max_num_iterations;
- scoped_ptr<LinearSolver> cg_solver(LinearSolver::Create(cg_options));
+ ConjugateGradientsSolver cg_solver(cg_options);
LinearSolver::PerSolveOptions cg_per_solve_options;
cg_per_solve_options.r_tolerance = per_solve_options.r_tolerance;
@@ -129,10 +130,10 @@
cg_summary.termination_type = FAILURE;
if (is_preconditioner_good) {
- cg_summary = cg_solver->Solve(schur_complement_.get(),
- schur_complement_->rhs().data(),
- cg_per_solve_options,
- reduced_linear_system_solution_.data());
+ cg_summary = cg_solver.Solve(schur_complement_.get(),
+ schur_complement_->rhs().data(),
+ cg_per_solve_options,
+ reduced_linear_system_solution_.data());
if (cg_summary.termination_type != FAILURE) {
schur_complement_->BackSubstitute(
reduced_linear_system_solution_.data(), x);
diff --git a/internal/ceres/linear_solver.cc b/internal/ceres/linear_solver.cc
index e3912eb..fd40252 100644
--- a/internal/ceres/linear_solver.cc
+++ b/internal/ceres/linear_solver.cc
@@ -31,7 +31,7 @@
#include "ceres/linear_solver.h"
#include <glog/logging.h>
-#include "ceres/conjugate_gradients_solver.h"
+#include "ceres/cgnr_solver.h"
#include "ceres/dense_qr_solver.h"
#include "ceres/iterative_schur_complement_solver.h"
#include "ceres/schur_complement_solver.h"
@@ -47,7 +47,7 @@
LinearSolver* LinearSolver::Create(const LinearSolver::Options& options) {
switch (options.type) {
case CONJUGATE_GRADIENTS:
- return new ConjugateGradientsSolver(options);
+ return new CgnrSolver(options);
case SPARSE_NORMAL_CHOLESKY:
#ifndef CERES_NO_SUITESPARSE
diff --git a/internal/ceres/linear_solver.h b/internal/ceres/linear_solver.h
index d287813..5860ecc 100644
--- a/internal/ceres/linear_solver.h
+++ b/internal/ceres/linear_solver.h
@@ -154,7 +154,7 @@
// In either case, x is the vector that solves the following
// optimization problem.
//
- // arg min_x ||Ax -b||^2 + ||Dx||^2
+ // arg min_x ||Ax - b||^2 + ||Dx||^2
//
// Here A is a matrix of size m x n, with full column rank. If A
// does not have full column rank, the results returned by the
diff --git a/internal/ceres/solver.cc b/internal/ceres/solver.cc
index c54668e..c45c55f 100644
--- a/internal/ceres/solver.cc
+++ b/internal/ceres/solver.cc
@@ -155,7 +155,8 @@
LinearSolverTypeToString(linear_solver_type_given),
LinearSolverTypeToString(linear_solver_type_used));
- if (linear_solver_type_given == ITERATIVE_SCHUR) {
+ if (linear_solver_type_given == CONJUGATE_GRADIENTS ||
+ linear_solver_type_given == ITERATIVE_SCHUR) {
internal::StringAppendF(&report, "Preconditioner %25s%25s\n",
PreconditionerTypeToString(preconditioner_type),
PreconditionerTypeToString(preconditioner_type));
diff --git a/internal/ceres/solver_impl.cc b/internal/ceres/solver_impl.cc
index ff2f5ea..ed07d9d 100644
--- a/internal/ceres/solver_impl.cc
+++ b/internal/ceres/solver_impl.cc
@@ -431,12 +431,6 @@
LinearSolver* SolverImpl::CreateLinearSolver(Solver::Options* options,
string* error) {
- if (options->linear_solver_type == CONJUGATE_GRADIENTS) {
- *error = "CONJUGATE_GRADIENTS is not a valid solver for "
- "linear least squares problems.";
- return NULL;
- }
-
#ifdef CERES_NO_SUITESPARSE
if (options->linear_solver_type == SPARSE_NORMAL_CHOLESKY) {
*error = "Can't use SPARSE_NORMAL_CHOLESKY because SuiteSparse was not "
