Covariance estimation using SuiteSparseQR.

Change-Id: I70d1686e3288fdde5f9723e832e15ffb857d6d85

CMakeLists.txt

@@ -215,6 +215,23 @@ ELSE (EXISTS ${CHOLMOD_INCLUDE})
   SET(CHOLMOD_FOUND FALSE)
 ENDIF (EXISTS ${CHOLMOD_INCLUDE})
 
+SET(SUITESPARSEQR_FOUND TRUE)
+FIND_LIBRARY(SUITESPARSEQR_LIB NAMES spqr PATHS ${SUITESPARSE_SEARCH_LIBS})
+IF (EXISTS ${SUITESPARSEQR_LIB})
+  MESSAGE("-- Found SUITESPARSEQR library: ${SUITESPARSEQR_LIB}")
+ELSE (EXISTS ${SUITESPARSEQR_LIB})
+  MESSAGE("-- Did not find SUITESPARSEQR library")
+  SET(SUITESPARSEQR_FOUND FALSE)
+ENDIF (EXISTS ${SUITESPARSEQR_LIB})
+
+FIND_PATH(SUITESPARSEQR_INCLUDE NAMES SuiteSparseQR.hpp PATHS ${SUITESPARSE_SEARCH_HEADERS})
+IF (EXISTS ${SUITESPARSEQR_INCLUDE})
+  MESSAGE("-- Found SUITESPARSEQR header in: ${SUITESPARSEQR_INCLUDE}")
+ELSE (EXISTS ${SUITESPARSEQR_INCLUDE})
+  MESSAGE("-- Did not find SUITESPARSEQR header")
+  SET(SUITESPARSEQR_FOUND FALSE)
+ENDIF (EXISTS ${SUITESPARSEQR_INCLUDE})
+
 # If SuiteSparse version is >= 4 then SuiteSparse_config is required.
 # For SuiteSparse 3, UFconfig.h is required.
 SET(SUITESPARSE_CONFIG_FOUND TRUE)
@@ -259,6 +276,24 @@ ELSE (EXISTS ${METIS_LIB})
   MESSAGE("-- Did not find METIS library")
 ENDIF (EXISTS ${METIS_LIB})
 
+# SuiteSparseQR may be compiled with Intel Threading Building Blocks.
+SET(TBB_FOUND TRUE)
+FIND_LIBRARY(TBB_LIB NAMES tbb PATHS ${SEARCH_LIBS})
+IF (EXISTS ${TBB_LIB})
+  MESSAGE("-- Found TBB library: ${TBB_LIB}")
+ELSE (EXISTS ${TBB_LIB})
+  MESSAGE("-- Did not find TBB library")
+  SET(TBB_FOUND FALSE)
+ENDIF (EXISTS ${TBB_LIB})
+
+FIND_LIBRARY(TBB_MALLOC_LIB NAMES tbbmalloc PATHS ${SEARCH_LIBS})
+IF (EXISTS ${TBB_MALLOC_LIB})
+  MESSAGE("-- Found TBB Malloc library: ${TBB_MALLOC_LIB}")
+ELSE (EXISTS ${TBB_MALLOC_LIB})
+  MESSAGE("-- Did not find TBB library")
+  SET(TBB_FOUND FALSE)
+ENDIF (EXISTS ${TBB_MALLOC_LIB})
+
 SET(BLAS_AND_LAPACK_FOUND TRUE)
 IF (APPLE)
   # Mac OS X has LAPACK/BLAS bundled in a framework called
@@ -557,6 +592,7 @@ IF (SUITESPARSE)
   INCLUDE_DIRECTORIES(${COLAMD_INCLUDE})
   INCLUDE_DIRECTORIES(${CCOLAMD_INCLUDE})
   INCLUDE_DIRECTORIES(${CHOLMOD_INCLUDE})
+  INCLUDE_DIRECTORIES(${SUITESPARSEQR_INCLUDE})
   IF (SUITESPARSE_CONFIG_FOUND)
     INCLUDE_DIRECTORIES(${SUITESPARSE_CONFIG_INCLUDE})
   ENDIF (SUITESPARSE_CONFIG_FOUND)

include/ceres/covariance.h

@@ -200,34 +200,76 @@ class Covariance {
  public:
   struct Options {
     Options()
-        : num_threads(1),
 #ifndef CERES_NO_SUITESPARSE
-          use_dense_linear_algebra(false),
+        : algorithm_type(SPARSE_QR),
 #else
-          use_dense_linear_algebra(true),
+        : algorithm_type(DENSE_SVD),
 #endif
           min_reciprocal_condition_number(1e-14),
           null_space_rank(0),
+          num_threads(1),
           apply_loss_function(true) {
     }
 
-    // Number of threads to be used for evaluating the Jacobian and
-    // estimation of covariance.
-    int num_threads;
-
+    // Ceres supports three different algorithms for covariance
+    // estimation, which represent different tradeoffs in speed,
+    // accuracy and reliability.
+    //
+    // 1. DENSE_SVD uses Eigen's JacobiSVD to perform the
+    //    computations. It computes the singular value decomposition
+    //
+    //      U * S * V' = J
+    //
+    //    and then uses it to compute the pseudo inverse of J'J as
+    //
+    //      pseudoinverse[J'J]^ = V * pseudoinverse[S] * V'
+    //
+    //    It is an accurate but slow method and should only be used
+    //    for small to moderate sized problems. It can handle
+    //    full-rank as well as rank deficient Jacobians.
+    //
+    // 2. SPARSE_CHOLESKY uses the CHOLMOD sparse Cholesky
+    //    factorization library to compute the decomposition :
+    //
+    //      R'R = J'J
+    //
+    //    and then
+    //
+    //      [J'J]^-1  = [R'R]^-1
+    //
+    //    It a fast algorithm for sparse matrices that should be used
+    //    when the Jacobian matrix J is well conditioned. For
+    //    ill-conditioned matrices, this algorithm can fail
+    //    unpredictabily. This is because Cholesky factorization is
+    //    not a rank-revealing factorization, i.e., it cannot reliably
+    //    detect when the matrix being factorized is not of full
+    //    rank. SuiteSparse/CHOLMOD supplies a heuristic for checking
+    //    if the matrix is rank deficient (cholmod_rcond), but it is
+    //    only a heuristic and can have both false positive and false
+    //    negatives.
+    //
+    //    Recent versions of SuiteSparse (>= 4.2.0) provide a much
+    //    more efficient method for solving for rows of the covariance
+    //    matrix. Therefore, if you are doing SPARSE_CHOLESKY, we
+    //    strongly recommend using a recent version of SuiteSparse.
+    //
+    // 3. SPARSE_QR uses the SuiteSparseQR sparse QR factorization
+    //    library to compute the decomposition
+    //
+    //      Q * R = J
+    //
+    //    [J'J]^-1 = [R*R']^-1
+    //
+    //    It is a moderately fast algorithm for sparse matrices, which
+    //    at the price of more time and memory than the
+    //    SPARSE_CHOLESKY algorithm is numerically better behaved and
+    //    is rank revealing, i.e., it can reliably detect when the
+    //    Jacobian matrix is rank deficient.
+    //
+    // Neither SPARSE_CHOLESKY or SPARSE_QR are capable of computing
+    // the covariance if the Jacobian is rank deficient.
 
-    // When use_dense_linear_algebra = true, Eigen's JacobiSVD
-    // algorithm is used to perform the computations. It is an
-    // accurate but slow method and should only be used for small to
-    // moderate sized problems.
-    //
-    // When use_dense_linear_algebra = false, SuiteSparse/CHOLMOD is
-    // used to perform the computation. Recent versions of SuiteSparse
-    // (>= 4.2.0) provide a much more efficient method for solving for
-    // rows of the covariance matrix. Therefore, if you are doing
-    // large scale covariance estimation, we strongly recommend using
-    // a recent version of SuiteSparse.
-    bool use_dense_linear_algebra;
+    CovarianceAlgorithmType algorithm_type;
 
     // If the Jacobian matrix is near singular, then inverting J'J
     // will result in unreliable results, e.g, if
@@ -240,48 +282,46 @@ class Covariance {
     //   inv(J'J) = [ 2.0471e+14  -2.0471e+14]
     //              [-2.0471e+14   2.0471e+14]
     //
-    // This is not a useful result.
-    //
-    // The reciprocal condition number of a matrix is a measure of
-    // ill-conditioning or how close the matrix is to being
-    // singular/rank deficient. It is defined as the ratio of the
-    // smallest eigenvalue of the matrix to the largest eigenvalue. In
-    // the above case the reciprocal condition number is about
-    // 1e-16. Which is close to machine precision and even though the
-    // inverse exists, it is meaningless, and care should be taken to
-    // interpet the results of such an inversion.
-    //
-    // Matrices with condition number lower than
-    // min_reciprocal_condition_number are considered rank deficient
-    // and by default Covariance::Compute will return false if it
-    // encounters such a matrix.
-    //
-    // use_dense_linear_algebra = false
-    // --------------------------------
-    //
-    // When performing large scale sparse covariance estimation,
-    // computing the exact value of the reciprocal condition number is
-    // not possible as it would require computing the eigenvalues of
-    // J'J.
-    //
-    // In this case we use cholmod_rcond, which uses the ratio of the
-    // smallest to the largest diagonal entries of the Cholesky
-    // factorization as an approximation to the reciprocal condition
-    // number.
-    //
-    // However, care must be taken as this is a heuristic and can
-    // sometimes be a very crude estimate. The default value of
-    // min_reciprocal_condition_number has been set to a conservative
-    // value, and sometimes the Covariance::Compute may return false
-    // even if it is possible to estimate the covariance reliably. In
-    // such cases, the user should exercise their judgement before
-    // lowering the value of min_reciprocal_condition_number.
-    //
-    // use_dense_linear_algebra = true
-    // -------------------------------
-    //
-    // When using dense linear algebra, the user has more control in
-    // dealing with singular and near singular covariance matrices.
+    // This is not a useful result. Therefore, by default
+    // Covariance::Compute will return false if a rank deficient
+    // Jacobian is encountered. How rank deficiency is detected
+    // depends on the algorithm being used.
+    //
+    // 1. DENSE_SVD
+    //
+    //      min_sigma / max_sigma < sqrt(min_reciprocal_condition_number)
+    //
+    //    where min_sigma and max_sigma are the minimum and maxiumum
+    //    singular values of J respectively.
+    //
+    // 2. SPARSE_CHOLESKY
+    //
+    //      cholmod_rcond < min_reciprocal_conditioner_number
+    //
+    //    Here cholmod_rcond is a crude estimate of the reciprocal
+    //    condition number of J'J by using the maximum and minimum
+    //    diagonal entries of the Cholesky factor R. There are no
+    //    theoretical guarantees associated with this test. It can
+    //    give false positives and negatives. Use at your own
+    //    risk. The default value of min_reciprocal_condition_number
+    //    has been set to a conservative value, and sometimes the
+    //    Covariance::Compute may return false even if it is possible
+    //    to estimate the covariance reliably. In such cases, the user
+    //    should exercise their judgement before lowering the value of
+    //    min_reciprocal_condition_number.
+    //
+    // 3. SPARSE_QR
+    //
+    //      rank(J) < num_col(J)
+    //
+    //   Here rank(J) is the estimate of the rank of J returned by the
+    //   SuiteSparseQR algorithm. It is a fairly reliable indication
+    //   of rank deficiency.
+    //
+    double min_reciprocal_condition_number;
+
+    // When using DENSE_SVD, the user has more control in dealing with
+    // singular and near singular covariance matrices.
     //
     // As mentioned above, when the covariance matrix is near
     // singular, instead of computing the inverse of J'J, the
@@ -311,19 +351,12 @@ class Covariance {
     //
     //   lambda_i / lambda_max < min_reciprocal_condition_number.
     //
-    double min_reciprocal_condition_number;
-
-    // Truncate the smallest "null_space_rank" eigenvectors when
-    // computing the pseudo inverse of J'J.
-    //
-    // If null_space_rank = -1, then all eigenvectors with eigenvalues s.t.
-    //
-    //   lambda_i / lambda_max < min_reciprocal_condition_number.
-    //
-    // are dropped. See the documentation for
-    // min_reciprocal_condition_number for more details.
+    // This option has no effect on the SPARSE_CHOLESKY or SPARSE_QR
+    // algorithms.
     int null_space_rank;
 
+    int num_threads;
+
     // Even though the residual blocks in the problem may contain loss
     // functions, setting apply_loss_function to false will turn off
     // the application of the loss function to the output of the cost

include/ceres/types.h

@@ -383,6 +383,12 @@ enum LineSearchInterpolationType {
   CUBIC
 };
 
+enum CovarianceAlgorithmType {
+  DENSE_SVD,
+  SPARSE_CHOLESKY,
+  SPARSE_QR
+};
+
 const char* LinearSolverTypeToString(LinearSolverType type);
 bool StringToLinearSolverType(string value, LinearSolverType* type);
 
@@ -424,6 +430,12 @@ bool StringToLineSearchInterpolationType(
     string value,
     LineSearchInterpolationType* type);
 
+const char* CovarianceAlgorithmTypeToString(
+    CovarianceAlgorithmType type);
+bool StringToCovarianceAlgorithmType(
+    string value,
+    CovarianceAlgorithmType* type);
+
 const char* LinearSolverTerminationTypeToString(
     LinearSolverTerminationType type);
 

internal/ceres/CMakeLists.txt

@@ -152,6 +152,7 @@ IF (GFLAGS)
 ENDIF (GFLAGS)
 
 IF (SUITESPARSE_FOUND)
+  LIST(APPEND CERES_LIBRARY_DEPENDENCIES ${SUITESPARSEQR_LIB})
   LIST(APPEND CERES_LIBRARY_DEPENDENCIES ${CHOLMOD_LIB})
   LIST(APPEND CERES_LIBRARY_DEPENDENCIES ${CCOLAMD_LIB})
   LIST(APPEND CERES_LIBRARY_DEPENDENCIES ${CAMD_LIB})
@@ -170,8 +171,14 @@ IF (SUITESPARSE_FOUND)
   IF (EXISTS ${BLAS_LIB})
     LIST(APPEND CERES_LIBRARY_DEPENDENCIES ${BLAS_LIB})
   ENDIF (EXISTS ${BLAS_LIB})
+
+  IF (TBB_FOUND)
+    LIST(APPEND CERES_LIBRARY_DEPENDENCIES ${TBB_LIB})
+    LIST(APPEND CERES_LIBRARY_DEPENDENCIES ${TBB_MALLOC_LIB})
+  ENDIF (TBB_FOUND)
 ENDIF (SUITESPARSE_FOUND)
 
+
 IF (CXSPARSE_FOUND)
   LIST(APPEND CERES_LIBRARY_DEPENDENCIES ${CXSPARSE_LIB})
 ENDIF (CXSPARSE_FOUND)

internal/ceres/covariance_impl.cc

@@ -48,6 +48,7 @@
 #include "ceres/suitesparse.h"
 #include "ceres/wall_time.h"
 #include "glog/logging.h"
+#include "SuiteSparseQR.hpp"
 
 namespace ceres {
 namespace internal {
@@ -388,28 +389,34 @@ bool CovarianceImpl::ComputeCovarianceSparsity(
 }
 
 bool CovarianceImpl::ComputeCovarianceValues() {
-  if (options_.use_dense_linear_algebra) {
-    return ComputeCovarianceValuesUsingEigen();
-  }
-
+  switch (options_.algorithm_type) {
+    case (DENSE_SVD):
+      return ComputeCovarianceValuesUsingDenseSVD();
 #ifndef CERES_NO_SUITESPARSE
-  return ComputeCovarianceValuesUsingSuiteSparse();
-#else
-  LOG(ERROR) << "Ceres compiled without SuiteSparse. "
-             << "Large scale covariance computation is not possible.";
-  return false;
+    case (SPARSE_CHOLESKY):
+      return ComputeCovarianceValuesUsingSparseCholesky();
+    case (SPARSE_QR):
+      return ComputeCovarianceValuesUsingSparseQR();
 #endif
+    default:
+      LOG(ERROR) << "Unsupported covariance estimation algorithm type: "
+                 << CovarianceAlgorithmTypeToString(options_.algorithm_type);
+      return false;
+  }
+  return false;
 }
 
-bool CovarianceImpl::ComputeCovarianceValuesUsingSuiteSparse() {
+bool CovarianceImpl::ComputeCovarianceValuesUsingSparseCholesky() {
   EventLogger event_logger(
-      "CovarianceImpl::ComputeCovarianceValuesUsingSuiteSparse");
+      "CovarianceImpl::ComputeCovarianceValuesUsingSparseCholesky");
 #ifndef CERES_NO_SUITESPARSE
   if (covariance_matrix_.get() == NULL) {
     // Nothing to do, all zeros covariance matrix.
     return true;
   }
 
+  SuiteSparse ss;
+
   CRSMatrix jacobian;
   problem_->Evaluate(evaluate_options_, NULL, NULL, NULL, &jacobian);
 
@@ -431,12 +438,12 @@ bool CovarianceImpl::ComputeCovarianceValuesUsingSuiteSparse() {
   cholmod_jacobian_view.sorted = 1;
   cholmod_jacobian_view.packed = 1;
 
-  cholmod_factor* factor = ss_.AnalyzeCholesky(&cholmod_jacobian_view);
+  cholmod_factor* factor = ss.AnalyzeCholesky(&cholmod_jacobian_view);
   event_logger.AddEvent("Symbolic Factorization");
-  bool factorization_succeeded = ss_.Cholesky(&cholmod_jacobian_view, factor);
+  bool factorization_succeeded = ss.Cholesky(&cholmod_jacobian_view, factor);
   if (factorization_succeeded) {
     const double reciprocal_condition_number =
-        cholmod_rcond(factor, ss_.mutable_cc());
+        cholmod_rcond(factor, ss.mutable_cc());
     if (reciprocal_condition_number <
         options_.min_reciprocal_condition_number) {
       LOG(WARNING) << "Cholesky factorization of J'J is not reliable. "
@@ -450,7 +457,7 @@ bool CovarianceImpl::ComputeCovarianceValuesUsingSuiteSparse() {
 
   event_logger.AddEvent("Numeric Factorization");
   if (!factorization_succeeded) {
-    ss_.Free(factor);
+    ss.Free(factor);
     LOG(WARNING) << "Cholesky factorization failed.";
     return false;
   }
@@ -474,7 +481,7 @@ bool CovarianceImpl::ComputeCovarianceValuesUsingSuiteSparse() {
   // versions of SuiteSparse have the cholmod_solve2 function which
   // re-uses memory across calls.
 #if (SUITESPARSE_VERSION < 4002)
-  cholmod_dense* rhs = ss_.CreateDenseVector(NULL, num_rows, num_rows);
+  cholmod_dense* rhs = ss.CreateDenseVector(NULL, num_rows, num_rows);
   double* rhs_x = reinterpret_cast<double*>(rhs->x);
 
   for (int r = 0; r < num_rows; ++r) {
@@ -485,17 +492,17 @@ bool CovarianceImpl::ComputeCovarianceValuesUsingSuiteSparse() {
     }
 
     rhs_x[r] = 1.0;
-    cholmod_dense* solution = ss_.Solve(factor, rhs);
+    cholmod_dense* solution = ss.Solve(factor, rhs);
     double* solution_x = reinterpret_cast<double*>(solution->x);
     for (int idx = row_begin; idx < row_end; ++idx) {
       const int c = cols[idx];
       values[idx] = solution_x[c];
     }
-    ss_.Free(solution);
+    ss.Free(solution);
     rhs_x[r] = 0.0;
   }
 
-  ss_.Free(rhs);
+  ss.Free(rhs);
 #else  // SUITESPARSE_VERSION < 4002
 
   const int num_threads = options_.num_threads;
@@ -567,7 +574,7 @@ bool CovarianceImpl::ComputeCovarianceValuesUsingSuiteSparse() {
 
 #endif  // SUITESPARSE_VERSION < 4002
 
-  ss_.Free(factor);
+  ss.Free(factor);
   event_logger.AddEvent("Inversion");
   return true;
 
@@ -578,9 +585,144 @@ bool CovarianceImpl::ComputeCovarianceValuesUsingSuiteSparse() {
 #endif  // CERES_NO_SUITESPARSE
 };
 
-bool CovarianceImpl::ComputeCovarianceValuesUsingEigen() {
+bool CovarianceImpl::ComputeCovarianceValuesUsingSparseQR() {
   EventLogger event_logger(
-      "CovarianceImpl::ComputeCovarianceValuesUsingEigen");
+      "CovarianceImpl::ComputeCovarianceValuesUsingSparseQR");
+
+#ifndef CERES_NO_SUITESPARSE
+  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");
+
+  // Construct a compressed column form of the Jacobian.
+  const int num_rows = jacobian.num_rows;
+  const int num_cols = jacobian.num_cols;
+  const int num_nonzeros = jacobian.values.size();
+
+  vector<SuiteSparse_long> transpose_rows(num_cols + 1, 0);
+  vector<SuiteSparse_long> transpose_cols(num_nonzeros, 0);
+  vector<double> transpose_values(num_nonzeros, 0);
+
+  for (int idx = 0; idx < num_nonzeros; ++idx) {
+    transpose_rows[jacobian.cols[idx] + 1] += 1;
+  }
+
+  for (int i = 1; i < transpose_rows.size(); ++i) {
+    transpose_rows[i] += transpose_rows[i - 1];
+  }
+
+  for (int r = 0; r < num_rows; ++r) {
+    for (int idx = jacobian.rows[r]; idx < jacobian.rows[r + 1]; ++idx) {
+      const int c = jacobian.cols[idx];
+      const int transpose_idx = transpose_rows[c];
+      transpose_cols[transpose_idx] = r;
+      transpose_values[transpose_idx] = jacobian.values[idx];
+      ++transpose_rows[c];
+    }
+  }
+
+  for (int i = transpose_rows.size() - 1; i > 0 ; --i) {
+    transpose_rows[i] = transpose_rows[i - 1];
+  }
+  transpose_rows[0] = 0;
+
+  cholmod_sparse cholmod_jacobian;
+  cholmod_jacobian.nrow = num_rows;
+  cholmod_jacobian.ncol = num_cols;
+  cholmod_jacobian.nzmax = num_nonzeros;
+  cholmod_jacobian.nz = NULL;
+  cholmod_jacobian.p = reinterpret_cast<void*>(&transpose_rows[0]);
+  cholmod_jacobian.i = reinterpret_cast<void*>(&transpose_cols[0]);
+  cholmod_jacobian.x = reinterpret_cast<void*>(&transpose_values[0]);
+  cholmod_jacobian.z = NULL;
+  cholmod_jacobian.stype = 0;  // Matrix is not symmetric.
+  cholmod_jacobian.itype = CHOLMOD_LONG;
+  cholmod_jacobian.xtype = CHOLMOD_REAL;
+  cholmod_jacobian.dtype = CHOLMOD_DOUBLE;
+  cholmod_jacobian.sorted = 1;
+  cholmod_jacobian.packed = 1;
+
+  cholmod_common cc;
+  cholmod_l_start(&cc);
+
+  SuiteSparseQR_factorization<double>* factor =
+      SuiteSparseQR_factorize<double>(SPQR_ORDERING_BESTAMD,
+                                      SPQR_DEFAULT_TOL,
+                                      &cholmod_jacobian,
+                                      &cc);
+  event_logger.AddEvent("Numeric Factorization");
+
+  const int rank = cc.SPQR_istat[4];
+  if (rank < cholmod_jacobian.ncol) {
+    LOG(WARNING) << "Jacobian matrix is rank deficient."
+                 << "Number of columns: " << cholmod_jacobian.ncol
+                 << " rank: " << rank;
+    SuiteSparseQR_free(&factor, &cc);
+    cholmod_l_finish(&cc);
+    return false;
+  }
+
+  const int* rows = covariance_matrix_->rows();
+  const int* cols = covariance_matrix_->cols();
+  double* values = covariance_matrix_->mutable_values();
+
+  // 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.
+
+  cholmod_dense* rhs = cholmod_l_zeros(num_cols, 1, CHOLMOD_REAL, &cc);
+  double* rhs_x = reinterpret_cast<double*>(rhs->x);
+
+  for (int r = 0; r < num_cols; ++r) {
+    int row_begin = rows[r];
+    int row_end = rows[r + 1];
+    if (row_end == row_begin) {
+      continue;
+    }
+
+    rhs_x[r] = 1.0;
+
+    cholmod_dense* y1 = SuiteSparseQR_solve<double>(SPQR_RTX_EQUALS_ETB, factor, rhs, &cc);
+    cholmod_dense* solution = SuiteSparseQR_solve<double>(SPQR_RETX_EQUALS_B, factor, y1, &cc);
+
+    double* solution_x = reinterpret_cast<double*>(solution->x);
+    for (int idx = row_begin; idx < row_end; ++idx) {
+      const int c = cols[idx];
+      values[idx] = solution_x[c];
+    }
+
+    cholmod_l_free_dense(&y1, &cc);
+    cholmod_l_free_dense(&solution, &cc);
+    rhs_x[r] = 0.0;
+  }
+
+  cholmod_l_free_dense(&rhs, &cc);
+  SuiteSparseQR_free(&factor, &cc);
+  cholmod_l_finish(&cc);
+  event_logger.AddEvent("Inversion");
+  return true;
+
+#else  // CERES_NO_SUITESPARSE
+
+  return false;
+
+#endif  // CERES_NO_SUITESPARSE
+};
+
+bool CovarianceImpl::ComputeCovarianceValuesUsingDenseSVD() {
+  EventLogger event_logger(
+      "CovarianceImpl::ComputeCovarianceValuesUsingDenseSVD");
   if (covariance_matrix_.get() == NULL) {
     // Nothing to do, all zeros covariance matrix.
     return true;
@@ -602,6 +744,7 @@ bool CovarianceImpl::ComputeCovarianceValuesUsingEigen() {
 
   Eigen::JacobiSVD<Matrix> svd(dense_jacobian,
                                Eigen::ComputeThinU | Eigen::ComputeThinV);
+
   event_logger.AddEvent("SingularValueDecomposition");
 
   const Vector singular_values = svd.singularValues();

internal/ceres/covariance_impl.h

@@ -64,8 +64,9 @@ class CovarianceImpl {
       ProblemImpl* problem);
 
   bool ComputeCovarianceValues();
-  bool ComputeCovarianceValuesUsingSuiteSparse();
-  bool ComputeCovarianceValuesUsingEigen();
+  bool ComputeCovarianceValuesUsingSparseCholesky();
+  bool ComputeCovarianceValuesUsingSparseQR();
+  bool ComputeCovarianceValuesUsingDenseSVD();
 
   const CompressedRowSparseMatrix* covariance_matrix() const {
     return covariance_matrix_.get();
@@ -80,9 +81,6 @@ class CovarianceImpl {
   map<const double*, int> parameter_block_to_row_index_;
   set<const double*> constant_parameter_blocks_;
   scoped_ptr<CompressedRowSparseMatrix> covariance_matrix_;
-#ifndef CERES_NO_SUITESPARSE
-  SuiteSparse ss_;
-#endif
 };
 
 }  // namespace internal

internal/ceres/covariance_test.cc

@@ -400,11 +400,14 @@ TEST_F(CovarianceTest, NormalBehavior) {
   Covariance::Options options;
 
 #ifndef CERES_NO_SUITESPARSE
-  options.use_dense_linear_algebra = false;
+  options.algorithm_type = SPARSE_CHOLESKY;
+  ComputeAndCompareCovarianceBlocks(options, expected_covariance);
+
+  options.algorithm_type = SPARSE_QR;
   ComputeAndCompareCovarianceBlocks(options, expected_covariance);
 #endif
 
-  options.use_dense_linear_algebra = true;
+  options.algorithm_type = DENSE_SVD;
   ComputeAndCompareCovarianceBlocks(options, expected_covariance);
 }
 
@@ -445,11 +448,14 @@ TEST_F(CovarianceTest, ThreadedNormalBehavior) {
   options.num_threads = 4;
 
 #ifndef CERES_NO_SUITESPARSE
-  options.use_dense_linear_algebra = false;
+  options.algorithm_type = SPARSE_CHOLESKY;
+  ComputeAndCompareCovarianceBlocks(options, expected_covariance);
+
+  options.algorithm_type = SPARSE_QR;
   ComputeAndCompareCovarianceBlocks(options, expected_covariance);
 #endif
 
-  options.use_dense_linear_algebra = true;
+  options.algorithm_type = DENSE_SVD;
   ComputeAndCompareCovarianceBlocks(options, expected_covariance);
 }
 
@@ -491,11 +497,11 @@ TEST_F(CovarianceTest, ConstantParameterBlock) {
   Covariance::Options options;
 
 #ifndef CERES_NO_SUITESPARSE
-  options.use_dense_linear_algebra = false;
+  options.algorithm_type = SPARSE_CHOLESKY;
   ComputeAndCompareCovarianceBlocks(options, expected_covariance);
 #endif
 
-  options.use_dense_linear_algebra = true;
+  options.algorithm_type = DENSE_SVD;
   ComputeAndCompareCovarianceBlocks(options, expected_covariance);
 }
 
@@ -544,11 +550,11 @@ TEST_F(CovarianceTest, LocalParameterization) {
   Covariance::Options options;
 
 #ifndef CERES_NO_SUITESPARSE
-  options.use_dense_linear_algebra = false;
+  options.algorithm_type = SPARSE_CHOLESKY;
   ComputeAndCompareCovarianceBlocks(options, expected_covariance);
 #endif
 
-  options.use_dense_linear_algebra = true;
+  options.algorithm_type = DENSE_SVD;
   ComputeAndCompareCovarianceBlocks(options, expected_covariance);
 }
 
@@ -589,7 +595,7 @@ TEST_F(CovarianceTest, TruncatedRank) {
 
   {
     Covariance::Options options;
-    options.use_dense_linear_algebra = true;
+    options.algorithm_type = DENSE_SVD;
     // Force dropping of the smallest eigenvector.
     options.null_space_rank = 1;
     ComputeAndCompareCovarianceBlocks(options, expected_covariance);
@@ -597,7 +603,7 @@ TEST_F(CovarianceTest, TruncatedRank) {
 
   {
     Covariance::Options options;
-    options.use_dense_linear_algebra = true;
+    options.algorithm_type = DENSE_SVD;
     // Force dropping of the smallest eigenvector via the ratio but
     // automatic truncation.
     options.min_reciprocal_condition_number = 0.044494;
@@ -693,7 +699,7 @@ TEST_F(RankDeficientCovarianceTest, AutomaticTruncation) {
   };
 
   Covariance::Options options;
-  options.use_dense_linear_algebra = true;
+  options.algorithm_type = DENSE_SVD;
   options.null_space_rank = -1;
   ComputeAndCompareCovarianceBlocks(options, expected_covariance);
 }
@@ -723,9 +729,10 @@ class LargeScaleCovarianceTest : public ::testing::Test {
     }
   }
 
-  void ComputeAndCompare(int num_threads) {
+  void ComputeAndCompare(CovarianceAlgorithmType algorithm_type,
+                         int num_threads) {
     Covariance::Options options;
-    options.use_dense_linear_algebra = false;
+    options.algorithm_type = algorithm_type;
     options.num_threads = num_threads;
     Covariance covariance(options);
     EXPECT_TRUE(covariance.Compute(all_covariance_blocks_, &problem_));
@@ -768,7 +775,7 @@ class LargeScaleCovarianceTest : public ::testing::Test {
 #if !defined(CERES_NO_SUITESPARSE) && defined(CERES_USE_OPENMP)
 
 TEST_F(LargeScaleCovarianceTest, Parallel) {
-  ComputeAndCompare(4);
+  ComputeAndCompare(SPARSE_CHOLESKY, 4);
 }
 
 #endif  // !defined(CERES_NO_SUITESPARSE) && defined(CERES_USE_OPENMP)

internal/ceres/types.cc

@@ -239,6 +239,27 @@ bool StringToNonlinearConjugateGradientType(
   return false;
 }
 
+const char* CovarianceAlgorithmTypeToString(
+    CovarianceAlgorithmType type) {
+  switch (type) {
+    CASESTR(DENSE_SVD);
+    CASESTR(SPARSE_CHOLESKY);
+    CASESTR(SPARSE_QR);
+    default:
+      return "UNKNOWN";
+  }
+}
+
+bool StringToCovarianceAlgorithmType(
+    string value,
+    CovarianceAlgorithmType* type) {
+  UpperCase(&value);
+  STRENUM(DENSE_SVD);
+  STRENUM(SPARSE_CHOLESKY);
+  STRENUM(SPARSE_QR);
+  return false;
+}
+
 const char* SolverTerminationTypeToString(SolverTerminationType type) {
   switch (type) {
     CASESTR(NO_CONVERGENCE);