Address some of the comments on CGNR patch - Rename BlockDiagonalPreconditioner to BlockJacobiPreconditioner - Include the diagonal in the block jacobi preconditioner. - Better flag help for eta. - Enable test for CGNR - Rename CONJUGATE_GRADIENTS to CGNR. - etc.
diff --git a/examples/bundle_adjuster.cc b/examples/bundle_adjuster.cc index 885dbe2..0331f32 100644 --- a/examples/bundle_adjuster.cc +++ b/examples/bundle_adjuster.cc
@@ -97,8 +97,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 == "cgnr") { + options->linear_solver_type = ceres::CGNR; } else if (FLAGS_solver_type == "dense_qr") { // DENSE_QR is included here for completeness, but actually using // this option is a bad idea due to the amount of memory needed @@ -110,14 +110,14 @@ << FLAGS_solver_type; } - if (options->linear_solver_type == ceres::CONJUGATE_GRADIENTS) { + if (options->linear_solver_type == ceres::CGNR) { options->linear_solver_min_num_iterations = 5; if (FLAGS_preconditioner_type == "identity") { options->preconditioner_type = ceres::IDENTITY; } else if (FLAGS_preconditioner_type == "jacobi") { options->preconditioner_type = ceres::JACOBI; } else { - LOG(FATAL) << "For CONJUGATE_GRADIENTS, only identity and jacobian " + LOG(FATAL) << "For CGNR, only identity and jacobian " << "preconditioners are supported. Got: " << FLAGS_preconditioner_type; }
diff --git a/include/ceres/types.h b/include/ceres/types.h index ea41b4d..b83a266 100644 --- a/include/ceres/types.h +++ b/include/ceres/types.h
@@ -83,12 +83,8 @@ // problems. ITERATIVE_SCHUR, - // Symmetric positive definite solvers - - // This is not meant for direct client use; it is used under the - // hood while using ITERATIVE_SCHUR. Once there is a decent - // preconditioner, this will make sense for general sparse problems. - CONJUGATE_GRADIENTS + // Conjugate gradients on the normal equations. + CGNR }; enum PreconditionerType {
diff --git a/internal/ceres/CMakeLists.txt b/internal/ceres/CMakeLists.txt index 03d650a..f105fcd 100644 --- a/internal/ceres/CMakeLists.txt +++ b/internal/ceres/CMakeLists.txt
@@ -29,8 +29,8 @@ # Author: keir@google.com (Keir Mierle) SET(CERES_INTERNAL_SRC - block_diagonal_preconditioner.cc block_evaluate_preparer.cc + block_jacobi_preconditioner.cc block_jacobian_writer.cc block_random_access_dense_matrix.cc block_random_access_matrix.cc
diff --git a/internal/ceres/block_diagonal_preconditioner.h b/internal/ceres/block_diagonal_preconditioner.h deleted file mode 100644 index 561b0db..0000000 --- a/internal/ceres/block_diagonal_preconditioner.h +++ /dev/null
@@ -1,68 +0,0 @@ -// 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/block_diagonal_preconditioner.cc b/internal/ceres/block_jacobi_preconditioner.cc similarity index 74% rename from internal/ceres/block_diagonal_preconditioner.cc rename to internal/ceres/block_jacobi_preconditioner.cc index 0779a91..f043200 100644 --- a/internal/ceres/block_diagonal_preconditioner.cc +++ b/internal/ceres/block_jacobi_preconditioner.cc
@@ -28,22 +28,23 @@ // // Author: keir@google.com (Keir Mierle) -#include "ceres/block_diagonal_preconditioner.h" +#include "ceres/block_jacobi_preconditioner.h" #include "Eigen/Cholesky" #include "ceres/block_sparse_matrix.h" #include "ceres/block_structure.h" #include "ceres/casts.h" +#include "ceres/integral_types.h" #include "ceres/internal/eigen.h" namespace ceres { namespace internal { -BlockDiagonalPreconditioner::BlockDiagonalPreconditioner( +BlockJacobiPreconditioner::BlockJacobiPreconditioner( const LinearOperator& A) : block_structure_( - *(down_cast<const BlockSparseMatrix*>(&A)->block_structure())) { - + *(down_cast<const BlockSparseMatrix*>(&A)->block_structure())), + num_rows_(A.num_rows()) { // Calculate the amount of storage needed. int storage_needed = 0; for (int c = 0; c < block_structure_.cols.size(); ++c) { @@ -56,7 +57,7 @@ block_storage_.resize(storage_needed); // Put pointers to the storage in the offsets. - double *block_cursor = &block_storage_[0]; + 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; @@ -64,10 +65,10 @@ } } -BlockDiagonalPreconditioner::~BlockDiagonalPreconditioner() { +BlockJacobiPreconditioner::~BlockJacobiPreconditioner() { } -void BlockDiagonalPreconditioner::Update(const LinearOperator& matrix) { +void BlockJacobiPreconditioner::Update(const LinearOperator& matrix, const double* D) { const BlockSparseMatrix& A = *(down_cast<const BlockSparseMatrix*>(&matrix)); const CompressedRowBlockStructure* bs = A.block_structure(); @@ -86,20 +87,35 @@ MatrixRef(blocks_[cells[c].block_id], col_block_size, col_block_size).noalias() += m.transpose() * m; + + // TODO(keir): Figure out when the below expression is actually faster + // than doing the full rank update. The issue is that for smaller sizes, + // the rankUpdate() function is slower than the full product done above. + // + // On the typical bundling problems, the above product is ~5% faster. + // + // MatrixRef(blocks_[cells[c].block_id], + // col_block_size, + // col_block_size).selfadjointView<Eigen::Upper>().rankUpdate(m); + // } } - // Invert each block. + // Add the diagonal and 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)); + const int position = block_structure_.cols[c].position; + MatrixRef DD(blocks_[c], size, size); + + DD.diagonal() += ConstVectorRef(D + position, size).array().square().matrix(); + + DD = DD.selfadjointView<Eigen::Upper>() + .ldlt() + .solve(Matrix::Identity(size, size)); } } -void BlockDiagonalPreconditioner::RightMultiply(const double* x, double* y) const { +void BlockJacobiPreconditioner::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; @@ -110,7 +126,7 @@ } } -void BlockDiagonalPreconditioner::LeftMultiply(const double* x, double* y) const { +void BlockJacobiPreconditioner::LeftMultiply(const double* x, double* y) const { RightMultiply(x, y); }
diff --git a/internal/ceres/block_jacobi_preconditioner.h b/internal/ceres/block_jacobi_preconditioner.h new file mode 100644 index 0000000..91cfedd --- /dev/null +++ b/internal/ceres/block_jacobi_preconditioner.h
@@ -0,0 +1,84 @@ +// 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_JACOBI_PRECONDITIONER_H_ +#define CERES_INTERNAL_BLOCK_JACOBI_PRECONDITIONER_H_ + +#include <vector> +#include "ceres/linear_operator.h" + +namespace ceres { +namespace internal { + +class CompressedRowBlockStructure; +class LinearOperator; +class SparseMatrix; + +// A block Jacobi preconditioner. This is intended for use with conjugate +// gradients, or other iterative symmetric solvers. To use the preconditioner, +// create one by passing a BlockSparseMatrix as the linear operator "A" to the +// constructor. This fixes the sparsity pattern to the pattern of the matrix +// A^TA. +// +// Before each use of the preconditioner in a solve with conjugate gradients, +// update the matrix by running Update(A, D). The values of the matrix A are +// inspected to construct the preconditioner. The vector D is applied as the +// D^TD diagonal term. +class BlockJacobiPreconditioner : public LinearOperator { + public: + // A must remain valid while the BlockJacobiPreconditioner is. + BlockJacobiPreconditioner(const LinearOperator& A); + virtual ~BlockJacobiPreconditioner(); + + // Update the preconditioner with the values found in A. The sparsity pattern + // must match that of the A passed to the constructor. D is a vector that + // must have the same number of rows as A, and is applied as a diagonal in + // addition to the block diagonals of A. + void Update(const LinearOperator& A, const double* D); + + // LinearOperator interface. + virtual void RightMultiply(const double* x, double* y) const; + virtual void LeftMultiply(const double* x, double* y) const; + virtual int num_rows() const { return num_rows_; } + virtual int num_cols() const { return num_rows_; } + + private: + std::vector<double*> blocks_; + std::vector<double> block_storage_; + int num_rows_; + + // The block structure of the matrix this preconditioner is for (e.g. J). + const CompressedRowBlockStructure& block_structure_; +}; + +} // namespace internal +} // namespace ceres + +#endif // CERES_INTERNAL_BLOCK_JACOBI_PRECONDITIONER_H_
diff --git a/internal/ceres/cgnr_linear_operator.h b/internal/ceres/cgnr_linear_operator.h index 94767fb..f32d8d9 100644 --- a/internal/ceres/cgnr_linear_operator.h +++ b/internal/ceres/cgnr_linear_operator.h
@@ -61,12 +61,11 @@ // // 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 = 2A^TAx - 2A^T b + 2 D^TDx // 0 = A^TAx - A^T b + D^TDx // 0 = (A^TA + D^TD)x - A^T b // @@ -80,26 +79,25 @@ // 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()]) { + CgnrLinearOperator(const LinearOperator& A, const 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); + std::fill(z_.get(), z_.get() + A_.num_rows(), 0.0); // z = Ax - A_->RightMultiply(x, z_.get()); + A_.RightMultiply(x, z_.get()); // y = y + Atz - A_->LeftMultiply(z_.get(), y); + 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]; - } + int n = A_.num_cols(); + VectorRef(y, n).array() += ConstVectorRef(D_, n).array().square() * + ConstVectorRef(x, n).array(); } } @@ -107,12 +105,12 @@ RightMultiply(x, y); } - virtual int num_rows() const { return A_->num_cols(); } - virtual int num_cols() const { return A_->num_cols(); } + virtual int num_rows() const { return A_.num_cols(); } + virtual int num_cols() const { return A_.num_cols(); } private: - LinearOperator* A_; - double* D_; + const LinearOperator& A_; + const double* D_; scoped_array<double> z_; };
diff --git a/internal/ceres/cgnr_solver.cc b/internal/ceres/cgnr_solver.cc index 2c88e33..ccc8026 100644 --- a/internal/ceres/cgnr_solver.cc +++ b/internal/ceres/cgnr_solver.cc
@@ -30,10 +30,11 @@ #include "ceres/cgnr_solver.h" +#include "glog/logging.h" #include "ceres/linear_solver.h" #include "ceres/cgnr_linear_operator.h" #include "ceres/conjugate_gradients_solver.h" -#include "ceres/block_diagonal_preconditioner.h" +#include "ceres/block_jacobi_preconditioner.h" namespace ceres { namespace internal { @@ -57,19 +58,19 @@ 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_.reset(new BlockJacobiPreconditioner(*A)); } - jacobi_preconditioner_->Update(*A); + jacobi_preconditioner_->Update(*A, per_solve_options.D); cg_per_solve_options.preconditioner = jacobi_preconditioner_.get(); } else if (options_.preconditioner_type != IDENTITY) { - // TODO(keir): This should die somehow. + LOG(FATAL) << "CGNR only supports IDENTITY and JACOBI preconditioners."; } // Solve (AtA + DtD)x = z (= Atb). std::fill(x, x + A->num_cols(), 0.0); - CgnrLinearOperator AtApDtD(A, per_solve_options.D); + CgnrLinearOperator lhs(*A, per_solve_options.D); ConjugateGradientsSolver conjugate_gradient_solver(options_); - return conjugate_gradient_solver.Solve(&AtApDtD, + return conjugate_gradient_solver.Solve(&lhs, z.get(), cg_per_solve_options, x);
diff --git a/internal/ceres/cgnr_solver.h b/internal/ceres/cgnr_solver.h index db49b07..dd36f99 100644 --- a/internal/ceres/cgnr_solver.h +++ b/internal/ceres/cgnr_solver.h
@@ -37,7 +37,7 @@ namespace ceres { namespace internal { -class BlockDiagonalPreconditioner; +class BlockJacobiPreconditioner; // A conjugate gradients on the normal equations solver. This directly solves // for the solution to @@ -56,7 +56,7 @@ private: const LinearSolver::Options options_; - scoped_ptr<BlockDiagonalPreconditioner> jacobi_preconditioner_; + scoped_ptr<BlockJacobiPreconditioner> jacobi_preconditioner_; DISALLOW_COPY_AND_ASSIGN(CgnrSolver); };
diff --git a/internal/ceres/evaluator.cc b/internal/ceres/evaluator.cc index 7ceb7ac..ea05aef 100644 --- a/internal/ceres/evaluator.cc +++ b/internal/ceres/evaluator.cc
@@ -53,7 +53,7 @@ case DENSE_SCHUR: case SPARSE_SCHUR: case ITERATIVE_SCHUR: - case CONJUGATE_GRADIENTS: + case CGNR: return new ProgramEvaluator<BlockEvaluatePreparer, BlockJacobianWriter>(options, program);
diff --git a/internal/ceres/iterative_schur_complement_solver.cc b/internal/ceres/iterative_schur_complement_solver.cc index 4379ebd..42d7e86 100644 --- a/internal/ceres/iterative_schur_complement_solver.cc +++ b/internal/ceres/iterative_schur_complement_solver.cc
@@ -84,7 +84,6 @@ // Instantiate a conjugate gradient solver that runs on the Schur complement // matrix with the block diagonal of the matrix F'F as the preconditioner. LinearSolver::Options cg_options; - cg_options.type = CONJUGATE_GRADIENTS; cg_options.max_num_iterations = options_.max_num_iterations; ConjugateGradientsSolver cg_solver(cg_options); LinearSolver::PerSolveOptions cg_per_solve_options;
diff --git a/internal/ceres/linear_solver.cc b/internal/ceres/linear_solver.cc index fd40252..7221a74 100644 --- a/internal/ceres/linear_solver.cc +++ b/internal/ceres/linear_solver.cc
@@ -46,7 +46,7 @@ LinearSolver* LinearSolver::Create(const LinearSolver::Options& options) { switch (options.type) { - case CONJUGATE_GRADIENTS: + case CGNR: return new CgnrSolver(options); case SPARSE_NORMAL_CHOLESKY:
diff --git a/internal/ceres/solver.cc b/internal/ceres/solver.cc index c45c55f..77f04d1 100644 --- a/internal/ceres/solver.cc +++ b/internal/ceres/solver.cc
@@ -155,7 +155,7 @@ LinearSolverTypeToString(linear_solver_type_given), LinearSolverTypeToString(linear_solver_type_used)); - if (linear_solver_type_given == CONJUGATE_GRADIENTS || + if (linear_solver_type_given == CGNR || linear_solver_type_given == ITERATIVE_SCHUR) { internal::StringAppendF(&report, "Preconditioner %25s%25s\n", PreconditionerTypeToString(preconditioner_type),
diff --git a/internal/ceres/system_test.cc b/internal/ceres/system_test.cc index 01573ac..a9a6747 100644 --- a/internal/ceres/system_test.cc +++ b/internal/ceres/system_test.cc
@@ -476,6 +476,8 @@ CONFIGURE(DENSE_SCHUR, USER, IDENTITY, 1); CONFIGURE(DENSE_SCHUR, SCHUR, IDENTITY, 1); + CONFIGURE(CGNR, USER, JACOBI, 1); + CONFIGURE(ITERATIVE_SCHUR, USER, JACOBI, 1); #ifndef CERES_NO_SUITESPARSE
diff --git a/internal/ceres/types.cc b/internal/ceres/types.cc index 55910b0..860f8a4 100644 --- a/internal/ceres/types.cc +++ b/internal/ceres/types.cc
@@ -42,7 +42,7 @@ CASESTR(DENSE_SCHUR); CASESTR(SPARSE_SCHUR); CASESTR(ITERATIVE_SCHUR); - CASESTR(CONJUGATE_GRADIENTS); + CASESTR(CGNR); default: return "UNKNOWN"; }