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/examples/bundle_adjuster.cc b/examples/bundle_adjuster.cc index bafdf25..61cd41f 100644 --- a/examples/bundle_adjuster.cc +++ b/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_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 @@ 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 @@ << 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 @@ 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,
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 "