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
G-P-S · Apr 7, 2015 · e712ce1 · e712ce1
1 parent e78a97a
commit e712ce1
Show file tree

Hide file tree

Showing 13 changed files with 53 additions and 206 deletions.
diff --git 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
@@ -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
@@ -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 @@ 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<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
@@ -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 @@ 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
@@ -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,15 +95,19 @@ 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<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 {
     summary.termination_type = LINEAR_SOLVER_SUCCESS;
     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
diff --git a/internal/ceres/eigen_dense_cholesky.h b/internal/ceres/eigen_dense_cholesky.h
diff --git 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<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
@@ -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 @@ 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<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 @@ 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<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
@@ -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<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
@@ -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<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 @@ 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
-
   }
 }