Adds support for computing the covariance using Eigen's sparse QR module.
For smaller problems Eigen is faster than SuiteSparseQR. This has been
tested with Eigen 3.2.1. Below are detailed timings. Problem 1 is the
smallest and problem 3 is the largest. The timings below are:
mean +- standard deviation.
Problem 1:
Eigen 0.0009218 +- 0.0002755
SuiteSparse 0.001406 +- 0.001610
Problem 2:
Eigen 0.002338 +- 0.001005
SuiteSparse 0.001910 +- 0.0004513
Problem 3:
Eigen 0.005455 +- 0.001759
SuiteSparse 0.002411 +- 0.0004974
Detailed problem descriptions:
Problem 1 size:
Original Reduced
Parameter blocks 533 54
Parameters 368 104
Effective parameters 1201 94
Residual blocks 233 77
Residual 1194 258
Problem 2 size:
Original Reduced
Parameter blocks 573 84
Parameters 1458 184
Effective parameters 1281 164
Residual blocks 263 107
Residual 1314 378
Problem 3 size:
Original Reduced
Parameter blocks 613 114
Parameters 1548 264
Effective parameters 1361 234
Residual blocks 293 137
Residual 1434 498
Change-Id: I884a67e2f728fe2992812148d82ccf5f27864fd7
diff --git a/internal/ceres/covariance_impl.cc b/internal/ceres/covariance_impl.cc
index 75c80bf..437a03d 100644
--- a/internal/ceres/covariance_impl.cc
+++ b/internal/ceres/covariance_impl.cc
@@ -38,6 +38,8 @@
#include <cstdlib>
#include <utility>
#include <vector>
+#include "Eigen/SparseCore"
+#include "Eigen/SparseQR"
#include "Eigen/SVD"
#include "ceres/compressed_col_sparse_matrix_utils.h"
#include "ceres/compressed_row_sparse_matrix.h"
@@ -401,6 +403,8 @@
case SPARSE_QR:
return ComputeCovarianceValuesUsingSparseQR();
#endif
+ case EIGEN_SPARSE_QR:
+ return ComputeCovarianceValuesUsingEigenSparseQR();
default:
LOG(ERROR) << "Unsupported covariance estimation algorithm type: "
<< CovarianceAlgorithmTypeToString(options_.algorithm_type);
@@ -597,7 +601,7 @@
return false;
#endif // CERES_NO_SUITESPARSE
-};
+}
bool CovarianceImpl::ComputeCovarianceValuesUsingSparseQR() {
EventLogger event_logger(
@@ -851,7 +855,102 @@
}
event_logger.AddEvent("CopyToCovarianceMatrix");
return true;
-};
+}
+
+bool CovarianceImpl::ComputeCovarianceValuesUsingEigenSparseQR() {
+ EventLogger event_logger(
+ "CovarianceImpl::ComputeCovarianceValuesUsingEigenSparseQR");
+ if (covariance_matrix_.get() == NULL) {
+ // Nothing to do, all zeros covariance matrix.
+ return true;
+ }
+
+ CRSMatrix jacobian;
+ problem_->Evaluate(evaluate_options_, NULL, NULL, NULL, &jacobian);
+ event_logger.AddEvent("Evaluate");
+
+ typedef Eigen::SparseMatrix<double, Eigen::ColMajor> EigenSparseMatrix;
+
+ // Convert the matrix to column major order as required by SparseQR.
+ EigenSparseMatrix sparse_jacobian =
+ Eigen::MappedSparseMatrix<double, Eigen::RowMajor>(
+ jacobian.num_rows, jacobian.num_cols,
+ static_cast<int>(jacobian.values.size()),
+ jacobian.rows.data(), jacobian.cols.data(), jacobian.values.data());
+ event_logger.AddEvent("ConvertToSparseMatrix");
+
+ Eigen::SparseQR<EigenSparseMatrix, Eigen::COLAMDOrdering<int> >
+ qr_solver(sparse_jacobian);
+ event_logger.AddEvent("QRDecomposition");
+
+ if(qr_solver.info() != Eigen::Success) {
+ LOG(ERROR) << "Eigen::SparseQR decomposition failed.";
+ return false;
+ }
+
+ if (qr_solver.rank() < jacobian.num_cols) {
+ LOG(ERROR) << "Jacobian matrix is rank deficient. "
+ << "Number of columns: " << jacobian.num_cols
+ << " rank: " << qr_solver.rank();
+ return false;
+ }
+
+ const int* rows = covariance_matrix_->rows();
+ const int* cols = covariance_matrix_->cols();
+ double* values = covariance_matrix_->mutable_values();
+
+ // Compute the inverse column permutation used by QR factorization.
+ Eigen::PermutationMatrix<Eigen::Dynamic, Eigen::Dynamic> inverse_permutation =
+ qr_solver.colsPermutation().inverse();
+
+ // The following loop exploits the fact that the i^th column of A^{-1}
+ // is given by the solution to the linear system
+ //
+ // A x = e_i
+ //
+ // where e_i is a vector with e(i) = 1 and all other entries zero.
+ //
+ // Since the covariance matrix is symmetric, the i^th row and column
+ // are equal.
+ const int num_cols = jacobian.num_cols;
+ const int num_threads = options_.num_threads;
+ scoped_array<double> workspace(new double[num_threads * num_cols]);
+
+#pragma omp parallel for num_threads(num_threads) schedule(dynamic)
+ for (int r = 0; r < num_cols; ++r) {
+ const int row_begin = rows[r];
+ const int row_end = rows[r + 1];
+ if (row_end == row_begin) {
+ continue;
+ }
+
+# ifdef CERES_USE_OPENMP
+ int thread_id = omp_get_thread_num();
+# else
+ int thread_id = 0;
+# endif
+
+ double* solution = workspace.get() + thread_id * num_cols;
+ SolveRTRWithSparseRHS<int>(
+ num_cols,
+ qr_solver.matrixR().innerIndexPtr(),
+ qr_solver.matrixR().outerIndexPtr(),
+ &qr_solver.matrixR().data().value(0),
+ inverse_permutation.indices().coeff(r),
+ solution);
+
+ // Assign the values of the computed covariance using the
+ // inverse permutation used in the QR factorization.
+ for (int idx = row_begin; idx < row_end; ++idx) {
+ const int c = cols[idx];
+ values[idx] = solution[inverse_permutation.indices().coeff(c)];
+ }
+ }
+
+ event_logger.AddEvent("Inverse");
+
+ return true;
+}
} // namespace internal
} // namespace ceres
diff --git a/internal/ceres/covariance_impl.h b/internal/ceres/covariance_impl.h
index 0e7e217..ff640b9 100644
--- a/internal/ceres/covariance_impl.h
+++ b/internal/ceres/covariance_impl.h
@@ -67,6 +67,7 @@
bool ComputeCovarianceValuesUsingSparseCholesky();
bool ComputeCovarianceValuesUsingSparseQR();
bool ComputeCovarianceValuesUsingDenseSVD();
+ bool ComputeCovarianceValuesUsingEigenSparseQR();
const CompressedRowSparseMatrix* covariance_matrix() const {
return covariance_matrix_.get();
diff --git a/internal/ceres/covariance_test.cc b/internal/ceres/covariance_test.cc
index 471b223..6de3fec 100644
--- a/internal/ceres/covariance_test.cc
+++ b/internal/ceres/covariance_test.cc
@@ -409,6 +409,9 @@
options.algorithm_type = DENSE_SVD;
ComputeAndCompareCovarianceBlocks(options, expected_covariance);
+
+ options.algorithm_type = EIGEN_SPARSE_QR;
+ ComputeAndCompareCovarianceBlocks(options, expected_covariance);
}
#ifdef CERES_USE_OPENMP
@@ -457,6 +460,9 @@
options.algorithm_type = DENSE_SVD;
ComputeAndCompareCovarianceBlocks(options, expected_covariance);
+
+ options.algorithm_type = EIGEN_SPARSE_QR;
+ ComputeAndCompareCovarianceBlocks(options, expected_covariance);
}
#endif // CERES_USE_OPENMP
@@ -506,6 +512,9 @@
options.algorithm_type = DENSE_SVD;
ComputeAndCompareCovarianceBlocks(options, expected_covariance);
+
+ options.algorithm_type = EIGEN_SPARSE_QR;
+ ComputeAndCompareCovarianceBlocks(options, expected_covariance);
}
TEST_F(CovarianceTest, LocalParameterization) {
@@ -562,6 +571,9 @@
options.algorithm_type = DENSE_SVD;
ComputeAndCompareCovarianceBlocks(options, expected_covariance);
+
+ options.algorithm_type = EIGEN_SPARSE_QR;
+ ComputeAndCompareCovarianceBlocks(options, expected_covariance);
}
