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
diff --git a/internal/ceres/block_jacobi_preconditioner.cc b/internal/ceres/block_jacobi_preconditioner.cc
index 749e0b6..29974d4 100644
--- a/internal/ceres/block_jacobi_preconditioner.cc
+++ b/internal/ceres/block_jacobi_preconditioner.cc
@@ -111,7 +111,7 @@
}
block = block.selfadjointView<Eigen::Upper>()
- .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 96f5511..8e05dcc 100644
--- a/internal/ceres/dense_normal_cholesky_solver.cc
+++ b/internal/ceres/dense_normal_cholesky_solver.cc
@@ -81,7 +81,7 @@
summary.num_iterations = 1;
summary.termination_type = TOLERANCE;
VectorRef(x, num_cols) =
- lhs.selfadjointView<Eigen::Upper>().ldlt().solve(rhs);
+ lhs.selfadjointView<Eigen::Upper>().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 1fb9709..4ab75ab 100644
--- a/internal/ceres/dense_qr_solver.cc
+++ b/internal/ceres/dense_qr_solver.cc
@@ -73,7 +73,7 @@
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 7c934fb..32722bb 100644
--- a/internal/ceres/implicit_schur_complement.cc
+++ b/internal/ceres/implicit_schur_complement.cc
@@ -161,7 +161,7 @@
m = m
.selfadjointView<Eigen::Upper>()
- .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 bd36672..1694273 100644
--- a/internal/ceres/implicit_schur_complement_test.cc
+++ b/internal/ceres/implicit_schur_complement_test.cc
@@ -109,7 +109,7 @@
solution->setZero();
VectorRef schur_solution(solution->data() + num_cols_ - num_schur_rows,
num_schur_rows);
- schur_solution = lhs->selfadjointView<Eigen::Upper>().ldlt().solve(*rhs);
+ schur_solution = lhs->selfadjointView<Eigen::Upper>().llt().solve(*rhs);
eliminator->BackSubstitute(A_.get(), b_.get(), D,
schur_solution.data(), solution->data());
}
@@ -156,7 +156,7 @@
// Reference solution to the f_block.
const Vector reference_f_sol =
- lhs.selfadjointView<Eigen::Upper>().ldlt().solve(rhs);
+ 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 09f61d7..0df9304 100644
--- a/internal/ceres/schur_complement_solver.cc
+++ b/internal/ceres/schur_complement_solver.cc
@@ -135,7 +135,7 @@
VectorRef(solution, num_rows) =
ConstMatrixRef(m->values(), num_rows, num_rows)
.selfadjointView<Eigen::Upper>()
- .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 a7e96ae..bed8f3a 100644
--- a/internal/ceres/schur_eliminator_test.cc
+++ b/internal/ceres/schur_eliminator_test.cc
@@ -112,7 +112,7 @@
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 @@
.triangularView<Eigen::Upper>() = 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 @@
Vector reduced_sol =
lhs_ref
.selfadjointView<Eigen::Upper>()
- .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 aa840c5..338df71 100644
--- a/internal/ceres/schur_jacobi_preconditioner.cc
+++ b/internal/ceres/schur_jacobi_preconditioner.cc
@@ -128,7 +128,7 @@
VectorRef(y, block_size) =
block
.selfadjointView<Eigen::Upper>()
- .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 53d10e1..2edbb18 100644
--- a/internal/ceres/visibility_based_preconditioner_test.cc
+++ b/internal/ceres/visibility_based_preconditioner_test.cc
@@ -279,7 +279,7 @@
// preconditioner_->RightMultiply(x.data(), y.data());
// z = full_schur_complement
// .selfadjointView<Eigen::Upper>()
-// .ldlt().solve(x);
+// .llt().solve(x);
// double max_relative_difference =
// ((y - z).array() / z.array()).matrix().lpNorm<Eigen::Infinity>();
// EXPECT_NEAR(max_relative_difference, 0.0, kTolerance);
