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

include/ceres/types.h

@@ -400,7 +400,8 @@ enum LineSearchInterpolationType {
 enum CovarianceAlgorithmType {
   DENSE_SVD,
   SPARSE_CHOLESKY,
-  SPARSE_QR
+  SPARSE_QR,
+  EIGEN_SPARSE_QR
 };
 
 CERES_EXPORT const char* LinearSolverTypeToString(

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 @@ bool CovarianceImpl::ComputeCovarianceValues() {
     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 @@ bool CovarianceImpl::ComputeCovarianceValuesUsingSparseCholesky() {
   return false;
 
 #endif  // CERES_NO_SUITESPARSE
-};
+}
 
 bool CovarianceImpl::ComputeCovarianceValuesUsingSparseQR() {
   EventLogger event_logger(
@@ -851,7 +855,102 @@ bool CovarianceImpl::ComputeCovarianceValuesUsingDenseSVD() {
   }
   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

internal/ceres/covariance_impl.h

@@ -67,6 +67,7 @@ class CovarianceImpl {
   bool ComputeCovarianceValuesUsingSparseCholesky();
   bool ComputeCovarianceValuesUsingSparseQR();
   bool ComputeCovarianceValuesUsingDenseSVD();
+  bool ComputeCovarianceValuesUsingEigenSparseQR();
 
   const CompressedRowSparseMatrix* covariance_matrix() const {
     return covariance_matrix_.get();

internal/ceres/covariance_test.cc

@@ -409,6 +409,9 @@ TEST_F(CovarianceTest, NormalBehavior) {
 
   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 @@ TEST_F(CovarianceTest, ThreadedNormalBehavior) {
 
   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 @@ TEST_F(CovarianceTest, ConstantParameterBlock) {
 
   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 @@ TEST_F(CovarianceTest, LocalParameterization) {
 
   options.algorithm_type = DENSE_SVD;
   ComputeAndCompareCovarianceBlocks(options, expected_covariance);
+
+  options.algorithm_type = EIGEN_SPARSE_QR;
+  ComputeAndCompareCovarianceBlocks(options, expected_covariance);
 }