Add DenseQR Interface
1. Add EigenDenseQR & tests.
This implementation now uses an in place decomposition,
which means that we are not allocating, deallocating
memory every call.
2. Add LAPACKDenseQR and tests.
The LAPACK implementation instead of using dgels which is a
routine which does the factorization and solve in one
call, now uses dgeqrf for factorization and then
dormqr and dtrtrs for solving. This allows us to
have a factorize and solve interface like DenseCholesky.
And opens the door to iterative refinement and mixed
precision solves.
3. The refactor also allows us to simplify the interface to
DenseSparseMatrix considerably. The internals of this
class were complicated because we had the AppendDiagonal
and RemoveDiagonal methods and we did not want to allocate
deallocate memory every call. But since we pay the cost
of the copy anyways, we can just hold that buffer
in DenseQRSolver.
4. Delete lapack.cc/h
5. The net result is that everything seems to be a bit faster.
For LAPACK we are not doing some of the scaling work that
dgels was doing. For Eigen I think it maybe the inplace
decomposition.
Benchmark Time CPU Time Old Time New CPU Old CPU New
----------------------------------------------------------------------------------------------------------------------------------------------------------
BM_DenseSolver<ceres::EIGEN, ceres::DENSE_QR>/1/1 -0.1154 -0.1159 692 612 691 611
BM_DenseSolver<ceres::EIGEN, ceres::DENSE_QR>/2/1 -0.1601 -0.1553 717 603 712 601
BM_DenseSolver<ceres::EIGEN, ceres::DENSE_QR>/3/1 -0.1673 -0.1575 733 610 724 610
BM_DenseSolver<ceres::EIGEN, ceres::DENSE_QR>/6/2 -0.1008 -0.1003 886 797 884 796
BM_DenseSolver<ceres::EIGEN, ceres::DENSE_QR>/10/3 -0.1489 -0.1514 1283 1092 1281 1087
BM_DenseSolver<ceres::EIGEN, ceres::DENSE_QR>/12/4 -0.1040 -0.1104 1556 1394 1553 1381
BM_DenseSolver<ceres::EIGEN, ceres::DENSE_QR>/20/5 -0.0007 -0.0097 1911 1910 1908 1890
BM_DenseSolver<ceres::EIGEN, ceres::DENSE_QR>/40/5 -0.1033 -0.1022 2981 2673 2957 2655
BM_DenseSolver<ceres::EIGEN, ceres::DENSE_QR>/100/10 -0.0147 +0.0015 9275 9138 9026 9040
BM_DenseSolver<ceres::EIGEN, ceres::DENSE_QR>/200/10 -0.1408 -0.1284 15093 12968 14778 12880
BM_DenseSolver<ceres::EIGEN, ceres::DENSE_QR>/200/20 -0.0310 -0.0355 38973 37765 38837 37460
BM_DenseSolver<ceres::LAPACK, ceres::DENSE_QR>/1/1 -0.1228 -0.1256 736 646 731 640
BM_DenseSolver<ceres::LAPACK, ceres::DENSE_QR>/2/1 -0.1401 -0.1396 740 636 735 633
BM_DenseSolver<ceres::LAPACK, ceres::DENSE_QR>/3/1 -0.1731 -0.1695 744 615 738 613
BM_DenseSolver<ceres::LAPACK, ceres::DENSE_QR>/6/2 -0.1399 -0.1408 1121 965 1113 956
BM_DenseSolver<ceres::LAPACK, ceres::DENSE_QR>/10/3 -0.1110 -0.1145 1571 1397 1560 1382
BM_DenseSolver<ceres::LAPACK, ceres::DENSE_QR>/12/4 -0.1411 -0.1417 2006 1722 1993 1710
BM_DenseSolver<ceres::LAPACK, ceres::DENSE_QR>/20/5 -0.1740 -0.1729 2741 2264 2724 2253
BM_DenseSolver<ceres::LAPACK, ceres::DENSE_QR>/40/5 -0.0966 -0.1123 3462 3128 3425 3040
BM_DenseSolver<ceres::LAPACK, ceres::DENSE_QR>/100/10 -0.0387 -0.0998 10365 9964 10339 9307
BM_DenseSolver<ceres::LAPACK, ceres::DENSE_QR>/200/10 -0.2044 -0.2049 16031 12754 15998 12720
BM_DenseSolver<ceres::LAPACK, ceres::DENSE_QR>/200/20 -0.2391 -0.2386 35777 27223 35716 27193
Change-Id: I782f0d7664efe1435eebda92ddf47a0fe66c9c72
diff --git a/internal/ceres/dense_linear_solver_test.cc b/internal/ceres/dense_linear_solver_test.cc
index 419ec3e..4dff5af 100644
--- a/internal/ceres/dense_linear_solver_test.cc
+++ b/internal/ceres/dense_linear_solver_test.cc
@@ -89,24 +89,22 @@
solver->Solve(&lhs, rhs.data(), per_solve_options, solution.data());
EXPECT_EQ(summary.termination_type, LINEAR_SOLVER_SUCCESS);
- // If solving for the regularized solution, add the diagonal to the
- // matrix. This makes subsequent computations simpler.
- if (testing::get<2>(param)) {
- lhs.AppendDiagonal(problem->D.get());
- };
+ Vector normal_rhs = lhs.matrix().transpose() * rhs.head(num_rows);
+ Matrix normal_lhs = lhs.matrix().transpose() * lhs.matrix();
- Vector tmp = Vector::Zero(num_rows + num_cols);
- lhs.RightMultiply(solution.data(), tmp.data());
- Vector actual_normal_rhs = Vector::Zero(num_cols);
- lhs.LeftMultiply(tmp.data(), actual_normal_rhs.data());
+ if (regularized) {
+ ConstVectorRef diagonal(problem->D.get(), num_cols);
+ normal_lhs += diagonal.array().square().matrix().asDiagonal();
+ }
- Vector expected_normal_rhs = Vector::Zero(num_cols);
- lhs.LeftMultiply(rhs.data(), expected_normal_rhs.data());
- const double residual = (expected_normal_rhs - actual_normal_rhs).norm() /
- expected_normal_rhs.norm();
+ Vector actual_normal_rhs = normal_lhs * solution;
- EXPECT_NEAR(residual, 0.0, 10 * std::numeric_limits<double>::epsilon())
- << "\nexpected: " << expected_normal_rhs.transpose()
+ const double normalized_residual =
+ (normal_rhs - actual_normal_rhs).norm() / normal_rhs.norm();
+
+ EXPECT_NEAR(
+ normalized_residual, 0.0, 10 * std::numeric_limits<double>::epsilon())
+ << "\nexpected: " << normal_rhs.transpose()
<< "\nactual: " << actual_normal_rhs.transpose();
}
