diff --git a/include/ceres/covariance.h b/include/ceres/covariance.h
index b70f98b..83126b5 100644
--- a/include/ceres/covariance.h
+++ b/include/ceres/covariance.h
@@ -200,34 +200,76 @@
  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;
-
-
-    // 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.
+    // Ceres supports three different algorithms for covariance
+    // estimation, which represent different tradeoffs in speed,
+    // accuracy and reliability.
     //
-    // 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;
+    // 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.
+
+    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 @@
     //   inv(J'J) = [ 2.0471e+14  -2.0471e+14]
     //              [-2.0471e+14   2.0471e+14]
     //
-    // This is not a useful result.
+    // 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.
     //
-    // 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.
+    // 1. DENSE_SVD
     //
-    // 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.
+    //      min_sigma / max_sigma < sqrt(min_reciprocal_condition_number)
     //
-    // use_dense_linear_algebra = false
-    // --------------------------------
+    //    where min_sigma and max_sigma are the minimum and maxiumum
+    //    singular values of J respectively.
     //
-    // 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.
+    // 2. SPARSE_CHOLESKY
     //
-    // 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.
+    //      cholmod_rcond < min_reciprocal_conditioner_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.
+    //    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.
     //
-    // use_dense_linear_algebra = true
-    // -------------------------------
+    // 3. SPARSE_QR
     //
-    // When using dense linear algebra, the user has more control in
-    // dealing with singular and near singular covariance matrices.
+    //      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 @@
     //
     //   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
diff --git a/include/ceres/types.h b/include/ceres/types.h
index 46c2fd2..a967541 100644
--- a/include/ceres/types.h
+++ b/include/ceres/types.h
@@ -383,6 +383,12 @@
   CUBIC
 };
 
+enum CovarianceAlgorithmType {
+  DENSE_SVD,
+  SPARSE_CHOLESKY,
+  SPARSE_QR
+};
+
 const char* LinearSolverTypeToString(LinearSolverType type);
 bool StringToLinearSolverType(string value, LinearSolverType* type);
 
@@ -424,6 +430,12 @@
     string value,
     LineSearchInterpolationType* type);
 
+const char* CovarianceAlgorithmTypeToString(
+    CovarianceAlgorithmType type);
+bool StringToCovarianceAlgorithmType(
+    string value,
+    CovarianceAlgorithmType* type);
+
 const char* LinearSolverTerminationTypeToString(
     LinearSolverTerminationType type);