diff --git a/examples/bundle_adjuster.cc b/examples/bundle_adjuster.cc index cc8c6a28b..5173e7926 100644 --- a/examples/bundle_adjuster.cc +++ b/examples/bundle_adjuster.cc @@ -77,8 +77,8 @@ DEFINE_double(eta, 1e-2, "Default value for eta. Eta determines the " "accuracy of each linear solve of the truncated newton step. " "Changing this parameter can affect solve performance "); DEFINE_string(solver_type, "sparse_schur", "Options are: " - "sparse_schur, dense_schur, iterative_schur, cholesky, " - "dense_qr, and conjugate_gradients"); + "sparse_schur, dense_schur, iterative_schur, sparse_cholesky, " + "dense_qr, dense_cholesky and conjugate_gradients"); DEFINE_string(preconditioner_type, "jacobi", "Options are: " "identity, jacobi, schur_jacobi, cluster_jacobi, " "cluster_tridiagonal"); @@ -116,7 +116,7 @@ void SetLinearSolver(Solver::Options* options) { options->linear_solver_type = ceres::DENSE_SCHUR; } else if (FLAGS_solver_type == "iterative_schur") { options->linear_solver_type = ceres::ITERATIVE_SCHUR; - } else if (FLAGS_solver_type == "cholesky") { + } else if (FLAGS_solver_type == "sparse_cholesky") { options->linear_solver_type = ceres::SPARSE_NORMAL_CHOLESKY; } else if (FLAGS_solver_type == "cgnr") { options->linear_solver_type = ceres::CGNR; @@ -126,6 +126,12 @@ void SetLinearSolver(Solver::Options* options) { // to store even the smallest of the bundle adjustment jacobian // arrays options->linear_solver_type = ceres::DENSE_QR; + } else if (FLAGS_solver_type == "dense_cholesky") { + // DENSE_NORMAL_CHOLESKY is included here for completeness, but + // actually using this option is a bad idea due to the amount of + // memory needed to store even the smallest of the bundle + // adjustment jacobian arrays + options->linear_solver_type = ceres::DENSE_NORMAL_CHOLESKY; } else { LOG(FATAL) << "Unknown ceres solver type: " << FLAGS_solver_type; diff --git a/include/ceres/types.h b/include/ceres/types.h index 0dcb354e5..d6474cc5f 100644 --- a/include/ceres/types.h +++ b/include/ceres/types.h @@ -59,14 +59,18 @@ enum LinearSolverType { // normal equations A'A x = A'b. They are direct solvers and do not // assume any special problem structure. - // Solve the normal equations using a sparse cholesky solver; based - // on CHOLMOD. - SPARSE_NORMAL_CHOLESKY, + // Solve the normal equations using a dense Cholesky solver; based + // on Eigen. + DENSE_NORMAL_CHOLESKY, // Solve the normal equations using a dense QR solver; based on // Eigen. DENSE_QR, + // Solve the normal equations using a sparse cholesky solver; requires + // SuiteSparse or CXSparse. + SPARSE_NORMAL_CHOLESKY, + // Specialized solvers, specific to problems with a generalized // bi-partitite structure. diff --git a/internal/ceres/CMakeLists.txt b/internal/ceres/CMakeLists.txt index ef0bd6847..a9140051f 100644 --- a/internal/ceres/CMakeLists.txt +++ b/internal/ceres/CMakeLists.txt @@ -45,6 +45,7 @@ SET(CERES_INTERNAL_SRC conditioned_cost_function.cc conjugate_gradients_solver.cc corrector.cc + dense_normal_cholesky_solver.cc dense_qr_solver.cc dense_sparse_matrix.cc detect_structure.cc diff --git a/internal/ceres/dense_normal_cholesky_solver.cc b/internal/ceres/dense_normal_cholesky_solver.cc new file mode 100644 index 000000000..f6bb99abf --- /dev/null +++ b/internal/ceres/dense_normal_cholesky_solver.cc @@ -0,0 +1,86 @@ +// 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: sameeragarwal@google.com (Sameer Agarwal) + +#include "ceres/dense_normal_cholesky_solver.h" + +#include + +#include "Eigen/Dense" +#include "ceres/dense_sparse_matrix.h" +#include "ceres/linear_solver.h" +#include "ceres/internal/eigen.h" +#include "ceres/internal/scoped_ptr.h" +#include "ceres/types.h" + +namespace ceres { +namespace internal { + +DenseNormalCholeskySolver::DenseNormalCholeskySolver( + const LinearSolver::Options& options) + : options_(options) {} + +LinearSolver::Summary DenseNormalCholeskySolver::SolveImpl( + DenseSparseMatrix* A, + const double* b, + const LinearSolver::PerSolveOptions& per_solve_options, + double* x) { + const int num_rows = A->num_rows(); + const int num_cols = A->num_cols(); + + ConstAlignedMatrixRef Aref = A->matrix(); + Matrix lhs(num_cols, num_cols); + lhs.setZero(); + + // lhs += A'A + // + // Using rankUpdate instead of GEMM, exposes the fact that its the + // same matrix being multiplied with itself and that the product is + // symmetric. + lhs.selfadjointView().rankUpdate(Aref.transpose()); + + // rhs = A'b + Vector rhs = Aref.transpose() * ConstVectorRef(b, num_rows); + + if (per_solve_options.D != NULL) { + ConstVectorRef D(per_solve_options.D, num_cols); + lhs += D.array().square().matrix().asDiagonal(); + } + + VectorRef(x, num_cols) = + lhs.selfadjointView().ldlt().solve(rhs); + + LinearSolver::Summary summary; + summary.num_iterations = 1; + summary.termination_type = TOLERANCE; + return summary; +} + +} // namespace internal +} // namespace ceres diff --git a/internal/ceres/dense_normal_cholesky_solver.h b/internal/ceres/dense_normal_cholesky_solver.h new file mode 100644 index 000000000..de4774058 --- /dev/null +++ b/internal/ceres/dense_normal_cholesky_solver.h @@ -0,0 +1,95 @@ +// 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: sameeragarwal@google.com (Sameer Agarwal) +// +// Solve dense rectangular systems Ax = b by forming the normal +// equations and solving them using the Cholesky factorization. + +#ifndef CERES_INTERNAL_DENSE_NORMAL_CHOLESKY_SOLVER_H_ +#define CERES_INTERNAL_DENSE_NORMAL_CHOLESKY_SOLVER_H_ + +#include "ceres/linear_solver.h" +#include "ceres/internal/macros.h" + +namespace ceres { +namespace internal { + +class DenseSparseMatrix; + +// This class implements the LinearSolver interface for solving +// rectangular/unsymmetric (well constrained) linear systems of the +// form +// +// Ax = b +// +// Since there does not usually exist a solution that satisfies these +// equations, the solver instead solves the linear least squares +// problem +// +// min_x |Ax - b|^2 +// +// Setting the gradient of the above optimization problem to zero +// gives us the normal equations +// +// A'Ax = A'b +// +// A'A is a positive definite matrix (hopefully), and the resulting +// linear system can be solved using Cholesky factorization. +// +// If the PerSolveOptions struct has a non-null array D, then the +// augmented/regularized linear system +// +// [ A ]x = [b] +// [ diag(D) ] [0] +// +// is solved. +// +// This class uses the LDLT factorization routines from the Eigen +// library. This solver always returns a solution, it is the user's +// responsibility to judge if the solution is good enough for their +// purposes. +class DenseNormalCholeskySolver: public DenseSparseMatrixSolver { + public: + explicit DenseNormalCholeskySolver(const LinearSolver::Options& options); + + private: + virtual LinearSolver::Summary SolveImpl( + DenseSparseMatrix* A, + const double* b, + const LinearSolver::PerSolveOptions& per_solve_options, + double* x); + + const LinearSolver::Options options_; + CERES_DISALLOW_COPY_AND_ASSIGN(DenseNormalCholeskySolver); +}; + +} // namespace internal +} // namespace ceres + +#endif // CERES_INTERNAL_DENSE_NORMAL_CHOLESKY_SOLVER_H_ diff --git a/internal/ceres/dense_qr_solver.cc b/internal/ceres/dense_qr_solver.cc index 21c5fb595..2b329ee0e 100644 --- a/internal/ceres/dense_qr_solver.cc +++ b/internal/ceres/dense_qr_solver.cc @@ -33,8 +33,8 @@ #include #include "Eigen/Dense" +#include "ceres/dense_sparse_matrix.h" #include "ceres/linear_solver.h" -#include "ceres/triplet_sparse_matrix.h" #include "ceres/internal/eigen.h" #include "ceres/internal/scoped_ptr.h" #include "ceres/types.h" diff --git a/internal/ceres/dense_qr_solver.h b/internal/ceres/dense_qr_solver.h index 2e0b2837e..dd683a8c4 100644 --- a/internal/ceres/dense_qr_solver.h +++ b/internal/ceres/dense_qr_solver.h @@ -28,7 +28,7 @@ // // Author: sameeragarwal@google.com (Sameer Agarwal) // -// Solve dense rectangular systems Ax = b using the QR factoriztion. +// Solve dense rectangular systems Ax = b using the QR factorization. #ifndef CERES_INTERNAL_DENSE_QR_SOLVER_H_ #define CERES_INTERNAL_DENSE_QR_SOLVER_H_ diff --git a/internal/ceres/evaluator.cc b/internal/ceres/evaluator.cc index 204198c31..a3ce6f04b 100644 --- a/internal/ceres/evaluator.cc +++ b/internal/ceres/evaluator.cc @@ -51,6 +51,7 @@ Evaluator* Evaluator::Create(const Evaluator::Options& options, string* error) { switch (options.linear_solver_type) { case DENSE_QR: + case DENSE_NORMAL_CHOLESKY: return new ProgramEvaluator(options, program); diff --git a/internal/ceres/linear_solver.cc b/internal/ceres/linear_solver.cc index 97dbba450..08c3ba110 100644 --- a/internal/ceres/linear_solver.cc +++ b/internal/ceres/linear_solver.cc @@ -31,6 +31,7 @@ #include "ceres/linear_solver.h" #include "ceres/cgnr_solver.h" +#include "ceres/dense_normal_cholesky_solver.h" #include "ceres/dense_qr_solver.h" #include "ceres/iterative_schur_complement_solver.h" #include "ceres/schur_complement_solver.h" @@ -78,6 +79,9 @@ LinearSolver* LinearSolver::Create(const LinearSolver::Options& options) { case DENSE_QR: return new DenseQRSolver(options); + case DENSE_NORMAL_CHOLESKY: + return new DenseNormalCholeskySolver(options); + default: LOG(FATAL) << "Unknown linear solver type :" << options.type; diff --git a/internal/ceres/solver_impl_test.cc b/internal/ceres/solver_impl_test.cc index a6d6aacfb..76ece8dc7 100644 --- a/internal/ceres/solver_impl_test.cc +++ b/internal/ceres/solver_impl_test.cc @@ -577,6 +577,11 @@ TEST(SolverImpl, CreateLinearSolverNormalOperation) { EXPECT_EQ(options.linear_solver_type, DENSE_QR); EXPECT_TRUE(solver.get() != NULL); + options.linear_solver_type = DENSE_NORMAL_CHOLESKY; + solver.reset(SolverImpl::CreateLinearSolver(&options, &error)); + EXPECT_EQ(options.linear_solver_type, DENSE_NORMAL_CHOLESKY); + EXPECT_TRUE(solver.get() != NULL); + #ifndef CERES_NO_SUITESPARSE options.linear_solver_type = SPARSE_NORMAL_CHOLESKY; options.sparse_linear_algebra_library = SUITE_SPARSE; diff --git a/internal/ceres/system_test.cc b/internal/ceres/system_test.cc index 7cdff2184..3dfbc0074 100644 --- a/internal/ceres/system_test.cc +++ b/internal/ceres/system_test.cc @@ -279,7 +279,8 @@ TEST(SystemTest, PowellsFunction) { sparse_linear_algebra_library, \ ordering)) - CONFIGURE(DENSE_QR, SUITE_SPARSE, NATURAL); + CONFIGURE(DENSE_QR, SUITE_SPARSE, NATURAL); + CONFIGURE(DENSE_NORMAL_CHOLESKY, SUITE_SPARSE, NATURAL); CONFIGURE(DENSE_SCHUR, SUITE_SPARSE, SCHUR); #ifndef CERES_NO_SUITESPARSE diff --git a/internal/ceres/types.cc b/internal/ceres/types.cc index 18fefad2a..2e950c51d 100644 --- a/internal/ceres/types.cc +++ b/internal/ceres/types.cc @@ -37,8 +37,9 @@ namespace ceres { const char* LinearSolverTypeToString(LinearSolverType solver_type) { switch (solver_type) { - CASESTR(SPARSE_NORMAL_CHOLESKY); + CASESTR(DENSE_NORMAL_CHOLESKY); CASESTR(DENSE_QR); + CASESTR(SPARSE_NORMAL_CHOLESKY); CASESTR(DENSE_SCHUR); CASESTR(SPARSE_SCHUR); CASESTR(ITERATIVE_SCHUR); diff --git a/internal/ceres/unsymmetric_linear_solver_test.cc b/internal/ceres/unsymmetric_linear_solver_test.cc index f13c5a47f..0b0d593cf 100644 --- a/internal/ceres/unsymmetric_linear_solver_test.cc +++ b/internal/ceres/unsymmetric_linear_solver_test.cc @@ -73,7 +73,8 @@ class UnsymmetricLinearSolverTest : public ::testing::Test { scoped_ptr transformed_A; - if (linear_solver_type == DENSE_QR) { + if (linear_solver_type == DENSE_QR || + linear_solver_type == DENSE_NORMAL_CHOLESKY) { transformed_A.reset(new DenseSparseMatrix(*A_)); } else if (linear_solver_type == SPARSE_NORMAL_CHOLESKY) { transformed_A.reset(new CompressedRowSparseMatrix(*A_)); @@ -118,6 +119,10 @@ TEST_F(UnsymmetricLinearSolverTest, DenseQR) { TestSolver(DENSE_QR, SUITE_SPARSE); } +TEST_F(UnsymmetricLinearSolverTest, DenseNormalCholesky) { + TestSolver(DENSE_NORMAL_CHOLESKY, SUITE_SPARSE); +} + #ifndef CERES_NO_SUITESPARSE TEST_F(UnsymmetricLinearSolverTest, SparseNormalCholeskyUsingSuiteSparse) { TestSolver(SPARSE_NORMAL_CHOLESKY, SUITE_SPARSE); diff --git a/jni/Android.mk b/jni/Android.mk index d107ec4dc..dcd247e66 100644 --- a/jni/Android.mk +++ b/jni/Android.mk @@ -105,6 +105,7 @@ LOCAL_SRC_FILES := $(CERES_SRC_PATH)/array_utils.cc \ $(CERES_SRC_PATH)/conditioned_cost_function.cc \ $(CERES_SRC_PATH)/conjugate_gradients_solver.cc \ $(CERES_SRC_PATH)/corrector.cc \ + $(CERES_SRC_PATH)/dense_normal_cholesky_solver.cc \ $(CERES_SRC_PATH)/dense_qr_solver.cc \ $(CERES_SRC_PATH)/dense_sparse_matrix.cc \ $(CERES_SRC_PATH)/detect_structure.cc \