From e712ce1810cf7ec80f58f3998d4091445f5fa906 Mon Sep 17 00:00:00 2001 From: Sameer Agarwal Date: Tue, 7 Apr 2015 14:11:10 -0700 Subject: [PATCH] Revert 81219ff. Eigen upstream was broken a little while ago, and it seemed to be the case that we needed a fix for using the LLT factorization on ARM. This has been fixed and AFAIK there are no stable eigen releases with this bug in it. For full gore, see http://eigen.tuxfamily.org/bz/show_bug.cgi?id=992 In light of the fix, the extra layer of indirection introduced earlier is not needed and we are reverting to normal programming. Change-Id: I16929d2145253b38339b573b27b6b8fabd523704 --- cmake/config.h.in | 4 -- internal/ceres/CMakeLists.txt | 1 - .../block_random_access_diagonal_matrix.cc | 9 ++- ...lock_random_access_diagonal_matrix_test.cc | 9 ++- .../ceres/dense_normal_cholesky_solver.cc | 10 ++- internal/ceres/eigen_dense_cholesky.cc | 65 ------------------ internal/ceres/eigen_dense_cholesky.h | 67 ------------------- internal/ceres/implicit_schur_complement.cc | 14 ++-- .../ceres/implicit_schur_complement_test.cc | 24 +++---- internal/ceres/schur_complement_solver.cc | 11 ++- internal/ceres/schur_eliminator_impl.h | 22 ++---- internal/ceres/schur_eliminator_test.cc | 16 ++--- jni/Android.mk | 7 -- 13 files changed, 53 insertions(+), 206 deletions(-) delete mode 100644 internal/ceres/eigen_dense_cholesky.cc delete mode 100644 internal/ceres/eigen_dense_cholesky.h diff --git a/cmake/config.h.in b/cmake/config.h.in index 41d629ee4..aea79216f 100644 --- a/cmake/config.h.in +++ b/cmake/config.h.in @@ -44,10 +44,6 @@ // If defined, use the LGPL code in Eigen. @CERES_USE_EIGEN_SPARSE@ -// If defined, use Eigen's LDLT factorization instead of LLT -// factorization. -@CERES_USE_LDLT_FOR_EIGEN_CHOLESKY@ - // If defined, Ceres was compiled without LAPACK. @CERES_NO_LAPACK@ diff --git a/internal/ceres/CMakeLists.txt b/internal/ceres/CMakeLists.txt index 6cac9c9d9..a06f3a5b9 100644 --- a/internal/ceres/CMakeLists.txt +++ b/internal/ceres/CMakeLists.txt @@ -61,7 +61,6 @@ SET(CERES_INTERNAL_SRC dogleg_strategy.cc dynamic_compressed_row_jacobian_writer.cc dynamic_compressed_row_sparse_matrix.cc - eigen_dense_cholesky.cc evaluator.cc file.cc gradient_checking_cost_function.cc diff --git a/internal/ceres/block_random_access_diagonal_matrix.cc b/internal/ceres/block_random_access_diagonal_matrix.cc index afa10def4..052690d18 100644 --- a/internal/ceres/block_random_access_diagonal_matrix.cc +++ b/internal/ceres/block_random_access_diagonal_matrix.cc @@ -34,7 +34,7 @@ #include #include #include -#include "ceres/eigen_dense_cholesky.h" +#include "Eigen/Dense" #include "ceres/internal/port.h" #include "ceres/internal/scoped_ptr.h" #include "ceres/stl_util.h" @@ -125,7 +125,12 @@ void BlockRandomAccessDiagonalMatrix::Invert() { double* values = tsm_->mutable_values(); for (int i = 0; i < blocks_.size(); ++i) { const int block_size = blocks_[i]; - InvertUpperTriangularUsingCholesky(block_size, values, values); + MatrixRef block(values, block_size, block_size); + block = + block + .selfadjointView() + .llt() + .solve(Matrix::Identity(block_size, block_size)); values += block_size * block_size; } } diff --git a/internal/ceres/block_random_access_diagonal_matrix_test.cc b/internal/ceres/block_random_access_diagonal_matrix_test.cc index 0bf9b7957..8fa3798c3 100644 --- a/internal/ceres/block_random_access_diagonal_matrix_test.cc +++ b/internal/ceres/block_random_access_diagonal_matrix_test.cc @@ -32,10 +32,10 @@ #include #include "ceres/block_random_access_diagonal_matrix.h" -#include "ceres/eigen_dense_cholesky.h" #include "ceres/internal/eigen.h" #include "glog/logging.h" #include "gtest/gtest.h" +#include "Eigen/Cholesky" namespace ceres { namespace internal { @@ -147,10 +147,9 @@ TEST_F(BlockRandomAccessDiagonalMatrixTest, Invert) { const TripletSparseMatrix* tsm = m_->matrix(); Matrix dense; tsm->ToDenseMatrix(&dense); - Matrix expected_inverse(dense.rows(), dense.rows()); - InvertUpperTriangularUsingCholesky(dense.rows(), - dense.data(), - expected_inverse.data()); + Matrix expected_inverse = + dense.llt().solve(Matrix::Identity(dense.rows(), dense.rows())); + m_->Invert(); tsm->ToDenseMatrix(&dense); diff --git a/internal/ceres/dense_normal_cholesky_solver.cc b/internal/ceres/dense_normal_cholesky_solver.cc index 8682ae778..b13cf3fc9 100644 --- a/internal/ceres/dense_normal_cholesky_solver.cc +++ b/internal/ceres/dense_normal_cholesky_solver.cc @@ -32,9 +32,9 @@ #include +#include "Eigen/Dense" #include "ceres/blas.h" #include "ceres/dense_sparse_matrix.h" -#include "ceres/eigen_dense_cholesky.h" #include "ceres/internal/eigen.h" #include "ceres/internal/scoped_ptr.h" #include "ceres/lapack.h" @@ -95,8 +95,11 @@ LinearSolver::Summary DenseNormalCholeskySolver::SolveUsingEigen( LinearSolver::Summary summary; summary.num_iterations = 1; - if (SolveUpperTriangularUsingCholesky(num_cols, lhs.data(), rhs.data(), x) - != Eigen::Success) { + summary.termination_type = LINEAR_SOLVER_SUCCESS; + Eigen::LLT llt = + lhs.selfadjointView().llt(); + + if (llt.info() != Eigen::Success) { summary.termination_type = LINEAR_SOLVER_FAILURE; summary.message = "Eigen LLT decomposition failed."; } else { @@ -104,6 +107,7 @@ LinearSolver::Summary DenseNormalCholeskySolver::SolveUsingEigen( summary.message = "Success."; } + VectorRef(x, num_cols) = llt.solve(rhs); event_logger.AddEvent("Solve"); return summary; } diff --git a/internal/ceres/eigen_dense_cholesky.cc b/internal/ceres/eigen_dense_cholesky.cc deleted file mode 100644 index 7b43ad294..000000000 --- a/internal/ceres/eigen_dense_cholesky.cc +++ /dev/null @@ -1,65 +0,0 @@ -#include "ceres/eigen_dense_cholesky.h" -#include "ceres/internal/eigen.h" -#include "ceres/internal/port.h" -#include "Eigen/Dense" - -namespace ceres { -namespace internal { - -using Eigen::Upper; - -Eigen::ComputationInfo -InvertUpperTriangularUsingCholesky(const int size, - const double* values, - double* inverse_values) { - ConstMatrixRef m(values, size, size); - MatrixRef inverse(inverse_values, size, size); - - // On ARM we have experienced significant numerical problems with - // Eigen's LLT implementation. Defining - // CERES_USE_LDLT_FOR_EIGEN_CHOLESKY switches to using the slightly - // more expensive but much more numerically well behaved LDLT - // factorization algorithm. - -#ifdef CERES_USE_LDLT_FOR_EIGEN_CHOLESKY - Eigen::LDLT cholesky = m.selfadjointView().ldlt(); -#else - Eigen::LLT cholesky = m.selfadjointView().llt(); -#endif - - inverse = cholesky.solve(Matrix::Identity(size, size)); - return cholesky.info(); -} - -Eigen::ComputationInfo -SolveUpperTriangularUsingCholesky(int size, - const double* lhs_values, - const double* rhs_values, - double* solution_values) { - ConstMatrixRef lhs(lhs_values, size, size); - - // On ARM we have experienced significant numerical problems with - // Eigen's LLT implementation. Defining - // CERES_USE_LDLT_FOR_EIGEN_CHOLESKY switches to using the slightly - // more expensive but much more numerically well behaved LDLT - // factorization algorithm. - -#ifdef CERES_USE_LDLT_FOR_EIGEN_CHOLESKY - Eigen::LDLT cholesky = lhs.selfadjointView().ldlt(); -#else - Eigen::LLT cholesky = lhs.selfadjointView().llt(); -#endif - - VectorRef solution(solution_values, size); - if (solution_values == rhs_values) { - cholesky.solveInPlace(solution); - } else { - ConstVectorRef rhs(rhs_values, size); - solution = cholesky.solve(rhs); - } - - return cholesky.info(); -} - -} // namespace internal -} // namespace ceres diff --git a/internal/ceres/eigen_dense_cholesky.h b/internal/ceres/eigen_dense_cholesky.h deleted file mode 100644 index 0d6b18a22..000000000 --- a/internal/ceres/eigen_dense_cholesky.h +++ /dev/null @@ -1,67 +0,0 @@ -// Ceres Solver - A fast non-linear least squares minimizer -// Copyright 2015 Google Inc. All rights reserved. -// http://ceres-solver.org/ -// -// 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) -// -// Wrappers around Eigen's dense Cholesky factorization -// routines. Normally, if CERES_USE_LDLT_FOR_EIGEN_CHOLESKY is not -// defined Eigen's LLT factorization is used, which is faster than the -// LDLT factorization, however on ARM the LLT factorization seems to -// give significantly worse results than LDLT and we are forced to use -// LDLT instead. -// -// These wrapper functions provide a level of indirection to deal with -// this switching and hide it from the rest of the code base. - -#ifndef CERES_INTERNAL_ARRAY_UTILS_H_ -#define CERES_INTERNAL_ARRAY_UTILS_H_ - -#include "ceres/internal/eigen.h" - -namespace ceres { -namespace internal { - -// Invert a matrix using Eigen's dense Cholesky factorization. values -// and inverse_values can point to the same array. -Eigen::ComputationInfo -InvertUpperTriangularUsingCholesky(int size, - const double* values, - double* inverse_values); - -// Solve a linear system using Eigen's dense Cholesky -// factorization. rhs_values and solution can point to the same array. -Eigen::ComputationInfo -SolveUpperTriangularUsingCholesky(int size, - const double* lhs, - const double* rhs, - double* solution); - -} // namespace internal -} // namespace ceres - -#endif // CERES_INTERNAL_ARRAY_UTILS_H_ diff --git a/internal/ceres/implicit_schur_complement.cc b/internal/ceres/implicit_schur_complement.cc index b943510cf..d05f03817 100644 --- a/internal/ceres/implicit_schur_complement.cc +++ b/internal/ceres/implicit_schur_complement.cc @@ -30,9 +30,9 @@ #include "ceres/implicit_schur_complement.h" +#include "Eigen/Dense" #include "ceres/block_sparse_matrix.h" #include "ceres/block_structure.h" -#include "ceres/eigen_dense_cholesky.h" #include "ceres/internal/eigen.h" #include "ceres/internal/scoped_ptr.h" #include "ceres/linear_solver.h" @@ -148,16 +148,18 @@ void ImplicitSchurComplement::AddDiagonalAndInvert( const int row_block_pos = block_diagonal_structure->rows[r].block.position; const int row_block_size = block_diagonal_structure->rows[r].block.size; const Cell& cell = block_diagonal_structure->rows[r].cells[0]; - double* block_values = block_diagonal->mutable_values() + cell.position; - MatrixRef m(block_values, row_block_size, row_block_size); + MatrixRef m(block_diagonal->mutable_values() + cell.position, + row_block_size, row_block_size); + if (D != NULL) { ConstVectorRef d(D + row_block_pos, row_block_size); m += d.array().square().matrix().asDiagonal(); } - InvertUpperTriangularUsingCholesky(row_block_size, - block_values, - block_values); + m = m + .selfadjointView() + .llt() + .solve(Matrix::Identity(row_block_size, row_block_size)); } } diff --git a/internal/ceres/implicit_schur_complement_test.cc b/internal/ceres/implicit_schur_complement_test.cc index 9bec91f3e..e586ea14e 100644 --- a/internal/ceres/implicit_schur_complement_test.cc +++ b/internal/ceres/implicit_schur_complement_test.cc @@ -31,10 +31,10 @@ #include "ceres/implicit_schur_complement.h" #include +#include "Eigen/Dense" #include "ceres/block_random_access_dense_matrix.h" #include "ceres/block_sparse_matrix.h" #include "ceres/casts.h" -#include "ceres/eigen_dense_cholesky.h" #include "ceres/internal/eigen.h" #include "ceres/internal/scoped_ptr.h" #include "ceres/linear_least_squares_problems.h" @@ -107,16 +107,11 @@ class ImplicitSchurComplementTest : public ::testing::Test { solution->resize(num_cols_); solution->setZero(); - double* schur_solution = solution->data() + num_cols_ - num_schur_rows; - SolveUpperTriangularUsingCholesky(num_schur_rows, - lhs->data(), - rhs->data(), - schur_solution); - eliminator->BackSubstitute(A_.get(), - b_.get(), - D, - schur_solution, - solution->data()); + VectorRef schur_solution(solution->data() + num_cols_ - num_schur_rows, + num_schur_rows); + schur_solution = lhs->selfadjointView().llt().solve(*rhs); + eliminator->BackSubstitute(A_.get(), b_.get(), D, + schur_solution.data(), solution->data()); } AssertionResult TestImplicitSchurComplement(double* D) { @@ -163,11 +158,8 @@ class ImplicitSchurComplementTest : public ::testing::Test { } // Reference solution to the f_block. - Vector reference_f_sol(rhs.rows()); - SolveUpperTriangularUsingCholesky(lhs.rows(), - lhs.data(), - rhs.data(), - reference_f_sol.data()); + const Vector reference_f_sol = + lhs.selfadjointView().llt().solve(rhs); // Backsubstituted solution from the implicit schur solver using the // reference solution to the f_block. diff --git a/internal/ceres/schur_complement_solver.cc b/internal/ceres/schur_complement_solver.cc index 6b1a45889..2491060dc 100644 --- a/internal/ceres/schur_complement_solver.cc +++ b/internal/ceres/schur_complement_solver.cc @@ -43,7 +43,6 @@ #include "ceres/conjugate_gradients_solver.h" #include "ceres/cxsparse.h" #include "ceres/detect_structure.h" -#include "ceres/eigen_dense_cholesky.h" #include "ceres/internal/eigen.h" #include "ceres/internal/scoped_ptr.h" #include "ceres/lapack.h" @@ -53,6 +52,7 @@ #include "ceres/triplet_sparse_matrix.h" #include "ceres/types.h" #include "ceres/wall_time.h" +#include "Eigen/Dense" #include "Eigen/SparseCore" namespace ceres { @@ -198,13 +198,18 @@ DenseSchurComplementSolver::SolveReducedLinearSystem( summary.num_iterations = 1; if (options().dense_linear_algebra_library_type == EIGEN) { - if (SolveUpperTriangularUsingCholesky(num_rows, m->values(), rhs(), solution) - != Eigen::Success) { + Eigen::LLT llt = + ConstMatrixRef(m->values(), num_rows, num_rows) + .selfadjointView() + .llt(); + if (llt.info() != Eigen::Success) { summary.termination_type = LINEAR_SOLVER_FAILURE; summary.message = "Eigen failure. Unable to perform dense Cholesky factorization."; return summary; } + + VectorRef(solution, num_rows) = llt.solve(ConstVectorRef(rhs(), num_rows)); } else { VectorRef(solution, num_rows) = ConstVectorRef(rhs(), num_rows); summary.termination_type = diff --git a/internal/ceres/schur_eliminator_impl.h b/internal/ceres/schur_eliminator_impl.h index 50f7e66d7..7e06806a5 100644 --- a/internal/ceres/schur_eliminator_impl.h +++ b/internal/ceres/schur_eliminator_impl.h @@ -57,7 +57,6 @@ #include "ceres/block_random_access_matrix.h" #include "ceres/block_sparse_matrix.h" #include "ceres/block_structure.h" -#include "ceres/eigen_dense_cholesky.h" #include "ceres/internal/eigen.h" #include "ceres/internal/fixed_array.h" #include "ceres/internal/scoped_ptr.h" @@ -269,11 +268,11 @@ Eliminate(const BlockSparseMatrix* A, // which case its much faster to compute the inverse once and // use it to multiply other matrices/vectors instead of doing a // Solve call over and over again. - typename EigenTypes::Matrix - inverse_ete(e_block_size, e_block_size); - InvertUpperTriangularUsingCholesky(e_block_size, - ete.data(), - inverse_ete.data()); + typename EigenTypes::Matrix inverse_ete = + ete + .template selfadjointView() + .llt() + .solve(Matrix::Identity(e_block_size, e_block_size)); // For the current chunk compute and update the rhs of the reduced // linear system. @@ -362,18 +361,7 @@ BackSubstitute(const BlockSparseMatrix* A, ete.data(), 0, 0, e_block_size, e_block_size); } - // On ARM we have experienced significant numerical problems with - // Eigen's LLT implementation. Defining - // CERES_USE_LDLT_FOR_EIGEN_CHOLESKY switches to using the slightly - // more expensive but much more numerically well behaved LDLT - // factorization algorithm. - -#ifdef CERES_USE_LDLT_FOR_EIGEN_CHOLESKY - ete.ldlt().solveInPlace(y_block); -#else ete.llt().solveInPlace(y_block); -#endif - } } diff --git a/internal/ceres/schur_eliminator_test.cc b/internal/ceres/schur_eliminator_test.cc index 8217e8c05..40bc904ba 100644 --- a/internal/ceres/schur_eliminator_test.cc +++ b/internal/ceres/schur_eliminator_test.cc @@ -35,7 +35,6 @@ #include "ceres/block_sparse_matrix.h" #include "ceres/casts.h" #include "ceres/detect_structure.h" -#include "ceres/eigen_dense_cholesky.h" #include "ceres/internal/eigen.h" #include "ceres/internal/scoped_ptr.h" #include "ceres/linear_least_squares_problems.h" @@ -122,10 +121,7 @@ class SchurEliminatorTest : public ::testing::Test { .triangularView() = R - Q.transpose() * P * Q; rhs_expected = g.tail(schur_size) - Q.transpose() * P * g.head(num_eliminate_cols); - SolveUpperTriangularUsingCholesky(H.rows(), - H.data(), - g.data(), - sol_expected.data()); + sol_expected = H.llt().solve(g); } void EliminateSolveAndCompare(const VectorRef& diagonal, @@ -161,11 +157,11 @@ class SchurEliminatorTest : public ::testing::Test { eliminator->Eliminate(A.get(), b.get(), diagonal.data(), &lhs, rhs.data()); MatrixRef lhs_ref(lhs.mutable_values(), lhs.num_rows(), lhs.num_cols()); - Vector reduced_sol(lhs.num_rows()); - SolveUpperTriangularUsingCholesky(lhs.num_cols(), - lhs.values(), - rhs.data(), - reduced_sol.data()); + Vector reduced_sol = + lhs_ref + .selfadjointView() + .llt() + .solve(rhs); // Solution to the linear least squares problem. Vector sol(num_cols); diff --git a/jni/Android.mk b/jni/Android.mk index 2a9d3374b..d49ab4a2c 100644 --- a/jni/Android.mk +++ b/jni/Android.mk @@ -105,12 +105,6 @@ LOCAL_CFLAGS := $(CERES_EXTRA_DEFINES) \ -DCERES_NO_CXSPARSE \ -DCERES_STD_UNORDERED_MAP -# On ARM we have experienced significant numerical problems with -# Eigen's LLT implementation. Defining -# CERES_USE_LDLT_FOR_EIGEN_CHOLESKY switches to using the slightly -# more expensive but much more numerically well behaved LDLT -# factorization algorithm. -LOCAL_CFLAGS += -DCERES_USE_LDLT_FOR_EIGEN_CHOLESKY # If the user did not enable threads in CERES_EXTRA_DEFINES, then add # CERES_NO_THREADS. @@ -150,7 +144,6 @@ LOCAL_SRC_FILES := $(CERES_SRC_PATH)/array_utils.cc \ $(CERES_SRC_PATH)/dogleg_strategy.cc \ $(CERES_SRC_PATH)/dynamic_compressed_row_jacobian_writer.cc \ $(CERES_SRC_PATH)/dynamic_compressed_row_sparse_matrix.cc \ - $(CERES_SRC_PATH)/eigen_dense_cholesky.cc \ $(CERES_SRC_PATH)/evaluator.cc \ $(CERES_SRC_PATH)/file.cc \ $(CERES_SRC_PATH)/gradient_checking_cost_function.cc \