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
diff --git a/cmake/config.h.in b/cmake/config.h.in
index 41d629e..aea7921 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 6cac9c9..a06f3a5 100644
--- a/internal/ceres/CMakeLists.txt
+++ b/internal/ceres/CMakeLists.txt
@@ -61,7 +61,6 @@
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 afa10de..052690d 100644
--- a/internal/ceres/block_random_access_diagonal_matrix.cc
+++ b/internal/ceres/block_random_access_diagonal_matrix.cc
@@ -34,7 +34,7 @@
#include <set>
#include <utility>
#include <vector>
-#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 @@
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<Eigen::Upper>()
+ .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 0bf9b79..8fa3798 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 <vector>
#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 @@
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 8682ae7..b13cf3f 100644
--- a/internal/ceres/dense_normal_cholesky_solver.cc
+++ b/internal/ceres/dense_normal_cholesky_solver.cc
@@ -32,9 +32,9 @@
#include <cstddef>
+#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 summary;
summary.num_iterations = 1;
- if (SolveUpperTriangularUsingCholesky(num_cols, lhs.data(), rhs.data(), x)
- != Eigen::Success) {
+ summary.termination_type = LINEAR_SOLVER_SUCCESS;
+ Eigen::LLT<Matrix, Eigen::Upper> llt =
+ lhs.selfadjointView<Eigen::Upper>().llt();
+
+ if (llt.info() != Eigen::Success) {
summary.termination_type = LINEAR_SOLVER_FAILURE;
summary.message = "Eigen LLT decomposition failed.";
} else {
@@ -104,6 +107,7 @@
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 7b43ad2..0000000
--- 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<Matrix, Upper> cholesky = m.selfadjointView<Upper>().ldlt();
-#else
- Eigen::LLT<Matrix, Upper> cholesky = m.selfadjointView<Upper>().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<Matrix, Upper> cholesky = lhs.selfadjointView<Upper>().ldlt();
-#else
- Eigen::LLT<Matrix, Upper> cholesky = lhs.selfadjointView<Upper>().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 0d6b18a..0000000
--- 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 b943510..d05f038 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 @@
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<Eigen::Upper>()
+ .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 9bec91f..e586ea1 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 <cstddef>
+#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 @@
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<Eigen::Upper>().llt().solve(*rhs);
+ eliminator->BackSubstitute(A_.get(), b_.get(), D,
+ schur_solution.data(), solution->data());
}
AssertionResult TestImplicitSchurComplement(double* D) {
@@ -163,11 +158,8 @@
}
// 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<Eigen::Upper>().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 6b1a458..2491060 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 @@
summary.num_iterations = 1;
if (options().dense_linear_algebra_library_type == EIGEN) {
- if (SolveUpperTriangularUsingCholesky(num_rows, m->values(), rhs(), solution)
- != Eigen::Success) {
+ Eigen::LLT<Matrix, Eigen::Upper> llt =
+ ConstMatrixRef(m->values(), num_rows, num_rows)
+ .selfadjointView<Eigen::Upper>()
+ .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 50f7e66..7e06806 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 @@
// 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<kEBlockSize, kEBlockSize>::Matrix
- inverse_ete(e_block_size, e_block_size);
- InvertUpperTriangularUsingCholesky(e_block_size,
- ete.data(),
- inverse_ete.data());
+ typename EigenTypes<kEBlockSize, kEBlockSize>::Matrix inverse_ete =
+ ete
+ .template selfadjointView<Eigen::Upper>()
+ .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 @@
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 8217e8c..40bc904 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 @@
.triangularView<Eigen::Upper>() = 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 @@
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<Eigen::Upper>()
+ .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 2a9d337..d49ab4a 100644
--- a/jni/Android.mk
+++ b/jni/Android.mk
@@ -105,12 +105,6 @@
-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 @@
$(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 \
