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.

examples/bundle_adjuster.cc

@@ -66,8 +66,8 @@
 DEFINE_string(input, "", "Input File name");
 
 DEFINE_string(solver_type, "sparse_schur", "Options are: "
-              "sparse_schur, dense_schur, iterative_schur, cholesky "
-              "and dense_qr");
+              "sparse_schur, dense_schur, iterative_schur, cholesky, "
+              "dense_qr, and conjugate_gradients");
 
 DEFINE_string(preconditioner_type, "jacobi", "Options are: "
               "identity, jacobi, schur_jacobi, cluster_jacobi, "
@@ -75,6 +75,7 @@ DEFINE_string(preconditioner_type, "jacobi", "Options are: "
 
 DEFINE_int32(num_iterations, 5, "Number of iterations");
 DEFINE_int32(num_threads, 1, "Number of threads");
+DEFINE_double(eta, 1e-2, "Default value for eta.");
 DEFINE_bool(use_schur_ordering, false, "Use automatic Schur ordering.");
 DEFINE_bool(use_quaternions, false, "If true, uses quaternions to represent "
             "rotations. If false, angle axis is used");
@@ -94,6 +95,8 @@ void SetLinearSolver(Solver::Options* options) {
     options->linear_solver_type = ceres::ITERATIVE_SCHUR;
   } else if (FLAGS_solver_type == "cholesky") {
     options->linear_solver_type = ceres::SPARSE_NORMAL_CHOLESKY;
+  } else if (FLAGS_solver_type == "conjugate_gradients") {
+    options->linear_solver_type = ceres::CONJUGATE_GRADIENTS;
   } else if (FLAGS_solver_type == "dense_qr") {
     // DENSE_QR is included here for completeness, but actually using
     // this opttion is a bad idea due to the amount of memory needed
@@ -105,7 +108,8 @@ void SetLinearSolver(Solver::Options* options) {
                << FLAGS_solver_type;
   }
 
-  if (options->linear_solver_type == ceres::ITERATIVE_SCHUR) {
+  if (options->linear_solver_type == ceres::ITERATIVE_SCHUR ||
+      options->linear_solver_type == ceres::CONJUGATE_GRADIENTS) {
     options->linear_solver_min_num_iterations = 5;
     if (FLAGS_preconditioner_type == "identity") {
       options->preconditioner_type = ceres::IDENTITY;
@@ -178,6 +182,7 @@ void SetMinimizerOptions(Solver::Options* options) {
   options->max_num_iterations = FLAGS_num_iterations;
   options->minimizer_progress_to_stdout = true;
   options->num_threads = FLAGS_num_threads;
+  options->eta = FLAGS_eta;
 }
 
 void SetSolverOptionsFromFlags(BALProblem* bal_problem,

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 @@ SET(CERES_INTERNAL_SRC
     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

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

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_

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_

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

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_

internal/ceres/evaluator.cc

@@ -53,6 +53,7 @@ Evaluator* Evaluator::Create(const Evaluator::Options& options,
     case DENSE_SCHUR:
     case SPARSE_SCHUR:
     case ITERATIVE_SCHUR:
+    case CONJUGATE_GRADIENTS:
       return new ProgramEvaluator<BlockEvaluatePreparer,
                                   BlockJacobianWriter>(options,
                                                        program);
@@ -60,11 +61,6 @@ Evaluator* Evaluator::Create(const Evaluator::Options& options,
       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;

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::Summary IterativeSchurComplementSolver::SolveImpl(
   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 @@ LinearSolver::Summary IterativeSchurComplementSolver::SolveImpl(
   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);

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() {
 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

internal/ceres/linear_solver.h

@@ -154,7 +154,7 @@ class LinearSolver {
     // 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

internal/ceres/solver.cc

@@ -155,7 +155,8 @@ string Solver::Summary::FullReport() const {
                           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));

internal/ceres/solver_impl.cc

@@ -431,12 +431,6 @@ Program* SolverImpl::CreateReducedProgram(Solver::Options* options,
 
 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 "