diff --git a/examples/bundle_adjuster.cc b/examples/bundle_adjuster.cc index bafdf25f4..61cd41fe5 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_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, diff --git a/internal/ceres/CMakeLists.txt b/internal/ceres/CMakeLists.txt index 7891138e5..03d650aa6 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 @@ 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 diff --git a/internal/ceres/block_diagonal_preconditioner.cc b/internal/ceres/block_diagonal_preconditioner.cc new file mode 100644 index 000000000..0779a91ab --- /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(&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(&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& 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() + .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 000000000..561b0db1c --- /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 +#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 blocks_; + std::vector 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 000000000..94767fb67 --- /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 +#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 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 000000000..2c88e33ea --- /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 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 000000000..db49b0717 --- /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 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 38c3fdb6d..7ceb7ac95 100644 --- a/internal/ceres/evaluator.cc +++ b/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(options, program); @@ -60,11 +61,6 @@ Evaluator* Evaluator::Create(const Evaluator::Options& options, return new ProgramEvaluator(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 48d8453b5..4379ebdbf 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::Summary IterativeSchurComplementSolver::SolveImpl( LinearSolver::Options cg_options; cg_options.type = CONJUGATE_GRADIENTS; cg_options.max_num_iterations = options_.max_num_iterations; - scoped_ptr 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); diff --git a/internal/ceres/linear_solver.cc b/internal/ceres/linear_solver.cc index e3912eb18..fd4025296 100644 --- a/internal/ceres/linear_solver.cc +++ b/internal/ceres/linear_solver.cc @@ -31,7 +31,7 @@ #include "ceres/linear_solver.h" #include -#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 diff --git a/internal/ceres/linear_solver.h b/internal/ceres/linear_solver.h index d28781321..5860ecc8a 100644 --- a/internal/ceres/linear_solver.h +++ b/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 diff --git a/internal/ceres/solver.cc b/internal/ceres/solver.cc index c54668ea0..c45c55f70 100644 --- a/internal/ceres/solver.cc +++ b/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)); diff --git a/internal/ceres/solver_impl.cc b/internal/ceres/solver_impl.cc index ff2f5ea82..ed07d9dc6 100644 --- a/internal/ceres/solver_impl.cc +++ b/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 "