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";
}
