From 080d1d04bdf722c3f602833c4c07ac1c5d26fcc0 Mon Sep 17 00:00:00 2001 From: Sameer Agarwal Date: Mon, 12 Aug 2013 16:28:37 -0700 Subject: [PATCH] Use more performant, less conservative Eigen solvers. colPivHouseholderQR -> householderQR ldlt -> llt. The resulting performance differences are significant enough to justify switching. LAPACK's dgels routine used for solving linear least squares problems does not use pivoting either. Similarly, we are not actually using the fact that the matrix being factorized can be indefinite when using LDLT factorization, so its not clear that the performance hit is worth it. These two changes result in Eigen being able to use blocking algorithms, which for Cholesky factorization, brings the performance closer to hardware optimized LAPACK. Similarly for dense QR factorization, on intel there is a 2x speedup. Change-Id: I4459ee0fc8eb87d58e2b299dfaa9e656d539dc5e --- internal/ceres/block_jacobi_preconditioner.cc | 2 +- internal/ceres/dense_normal_cholesky_solver.cc | 2 +- internal/ceres/dense_qr_solver.cc | 2 +- internal/ceres/implicit_schur_complement.cc | 2 +- internal/ceres/implicit_schur_complement_test.cc | 4 ++-- internal/ceres/schur_complement_solver.cc | 2 +- internal/ceres/schur_eliminator_test.cc | 6 +++--- internal/ceres/schur_jacobi_preconditioner.cc | 2 +- internal/ceres/visibility_based_preconditioner_test.cc | 2 +- 9 files changed, 12 insertions(+), 12 deletions(-) diff --git a/internal/ceres/block_jacobi_preconditioner.cc b/internal/ceres/block_jacobi_preconditioner.cc index 749e0b65d..29974d45b 100644 --- a/internal/ceres/block_jacobi_preconditioner.cc +++ b/internal/ceres/block_jacobi_preconditioner.cc @@ -111,7 +111,7 @@ bool BlockJacobiPreconditioner::UpdateImpl(const BlockSparseMatrix& A, } block = block.selfadjointView() - .ldlt() + .llt() .solve(Matrix::Identity(size, size)); } return true; diff --git a/internal/ceres/dense_normal_cholesky_solver.cc b/internal/ceres/dense_normal_cholesky_solver.cc index 96f55115a..8e05dcc12 100644 --- a/internal/ceres/dense_normal_cholesky_solver.cc +++ b/internal/ceres/dense_normal_cholesky_solver.cc @@ -81,7 +81,7 @@ LinearSolver::Summary DenseNormalCholeskySolver::SolveImpl( summary.num_iterations = 1; summary.termination_type = TOLERANCE; VectorRef(x, num_cols) = - lhs.selfadjointView().ldlt().solve(rhs); + lhs.selfadjointView().llt().solve(rhs); event_logger.AddEvent("Solve"); return summary; diff --git a/internal/ceres/dense_qr_solver.cc b/internal/ceres/dense_qr_solver.cc index 1fb9709b4..4ab75abcc 100644 --- a/internal/ceres/dense_qr_solver.cc +++ b/internal/ceres/dense_qr_solver.cc @@ -73,7 +73,7 @@ LinearSolver::Summary DenseQRSolver::SolveImpl( event_logger.AddEvent("Setup"); // Solve the system. - VectorRef(x, num_cols) = A->matrix().colPivHouseholderQr().solve(rhs_); + VectorRef(x, num_cols) = A->matrix().householderQr().solve(rhs_); event_logger.AddEvent("Solve"); if (per_solve_options.D != NULL) { diff --git a/internal/ceres/implicit_schur_complement.cc b/internal/ceres/implicit_schur_complement.cc index 7c934fb88..32722bb6e 100644 --- a/internal/ceres/implicit_schur_complement.cc +++ b/internal/ceres/implicit_schur_complement.cc @@ -161,7 +161,7 @@ void ImplicitSchurComplement::AddDiagonalAndInvert( m = m .selfadjointView() - .ldlt() + .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 bd36672b6..1694273cc 100644 --- a/internal/ceres/implicit_schur_complement_test.cc +++ b/internal/ceres/implicit_schur_complement_test.cc @@ -109,7 +109,7 @@ class ImplicitSchurComplementTest : public ::testing::Test { solution->setZero(); VectorRef schur_solution(solution->data() + num_cols_ - num_schur_rows, num_schur_rows); - schur_solution = lhs->selfadjointView().ldlt().solve(*rhs); + schur_solution = lhs->selfadjointView().llt().solve(*rhs); eliminator->BackSubstitute(A_.get(), b_.get(), D, schur_solution.data(), solution->data()); } @@ -156,7 +156,7 @@ class ImplicitSchurComplementTest : public ::testing::Test { // Reference solution to the f_block. const Vector reference_f_sol = - lhs.selfadjointView().ldlt().solve(rhs); + 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 09f61d717..0df930419 100644 --- a/internal/ceres/schur_complement_solver.cc +++ b/internal/ceres/schur_complement_solver.cc @@ -135,7 +135,7 @@ bool DenseSchurComplementSolver::SolveReducedLinearSystem(double* solution) { VectorRef(solution, num_rows) = ConstMatrixRef(m->values(), num_rows, num_rows) .selfadjointView() - .ldlt() + .llt() .solve(ConstVectorRef(rhs(), num_rows)); return true; diff --git a/internal/ceres/schur_eliminator_test.cc b/internal/ceres/schur_eliminator_test.cc index a7e96ae9e..bed8f3a26 100644 --- a/internal/ceres/schur_eliminator_test.cc +++ b/internal/ceres/schur_eliminator_test.cc @@ -112,7 +112,7 @@ class SchurEliminatorTest : public ::testing::Test { P.block(row, row, block_size, block_size) = P .block(row, row, block_size, block_size) - .ldlt() + .llt() .solve(Matrix::Identity(block_size, block_size)); row += block_size; } @@ -121,7 +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); - sol_expected = H.ldlt().solve(g); + sol_expected = H.llt().solve(g); } void EliminateSolveAndCompare(const VectorRef& diagonal, @@ -160,7 +160,7 @@ class SchurEliminatorTest : public ::testing::Test { Vector reduced_sol = lhs_ref .selfadjointView() - .ldlt() + .llt() .solve(rhs); // Solution to the linear least squares problem. diff --git a/internal/ceres/schur_jacobi_preconditioner.cc b/internal/ceres/schur_jacobi_preconditioner.cc index aa840c572..338df715c 100644 --- a/internal/ceres/schur_jacobi_preconditioner.cc +++ b/internal/ceres/schur_jacobi_preconditioner.cc @@ -128,7 +128,7 @@ void SchurJacobiPreconditioner::RightMultiply(const double* x, VectorRef(y, block_size) = block .selfadjointView() - .ldlt() + .llt() .solve(ConstVectorRef(x, block_size)); x += block_size; diff --git a/internal/ceres/visibility_based_preconditioner_test.cc b/internal/ceres/visibility_based_preconditioner_test.cc index 53d10e1ff..2edbb187e 100644 --- a/internal/ceres/visibility_based_preconditioner_test.cc +++ b/internal/ceres/visibility_based_preconditioner_test.cc @@ -279,7 +279,7 @@ namespace internal { // preconditioner_->RightMultiply(x.data(), y.data()); // z = full_schur_complement // .selfadjointView() -// .ldlt().solve(x); +// .llt().solve(x); // double max_relative_difference = // ((y - z).array() / z.array()).matrix().lpNorm(); // EXPECT_NEAR(max_relative_difference, 0.0, kTolerance);