diff --git a/include/ceres/covariance.h b/include/ceres/covariance.h
index 6533487..26f94f7 100644
--- a/include/ceres/covariance.h
+++ b/include/ceres/covariance.h
@@ -159,8 +159,12 @@
 // reconstruction is ambiguous upto a similarity transform. This is
 // known as a Gauge Ambiguity. Handling Gauges correctly requires the
 // use of SVD or custom inversion algorithms. For small problems the
-// user can use the dense algorithm. For more details see Morris,
-// Kanatani & Kanade's work the subject.
+// user can use the dense algorithm. For more details see
+//
+// Ken-ichi Kanatani, Daniel D. Morris: Gauges and gauge
+// transformations for uncertainty description of geometric structure
+// with indeterminacy. IEEE Transactions on Information Theory 47(5):
+// 2017-2028 (2001)
 //
 // Example Usage
 // =============
@@ -201,7 +205,7 @@
 #else
           use_dense_linear_algebra(true),
 #endif
-          min_singular_value_threshold(1e-8),
+          min_reciprocal_condition_number(1e-14),
           null_space_rank(0),
           apply_loss_function(true) {
     }
@@ -210,18 +214,87 @@
     // estimation of covariance.
     int num_threads;
 
-    // Use Eigen's JacobiSVD algorithm to compute the covariance. This
-    // is a very accurate but slow algorithm. The up side is that it
-    // can handle numerically rank deficient jacobians. This option
-    // only makes sense for small to moderate sized problems.
+    // Use Eigen's JacobiSVD algorithm to compute the covariance
+    // instead of SuiteSparse. This is a very accurate but slow
+    // algorithm. The up side is that it can handle numerically rank
+    // deficient jacobians. This option only makes sense for small to
+    // moderate sized problems.
     bool use_dense_linear_algebra;
 
-    // When use_dense_linear_algebra is true, singular values less
-    // than min_singular_value_threshold are set to zero.
-    double min_singular_value_threshold;
+    // If the Jacobian matrix is near singular, then inverting J'J
+    // will result in unreliable results, e.g,
+    //
+    //   J = [1.0 1.0         ]
+    //       [1.0 1.0000001   ]
+    //
+    // Which is essentially a rank deficient matrix
+    //
+    //   inv(J'J) = [ 2.0471e+14  -2.0471e+14]
+    //              [-2.0471e+14   2.0471e+14]
+    //
+    // 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.
+    //
+    // Depending on the value of use_dense_linear_algebra this may
+    // have further consequences on the covariance estimation process.
+    //
+    // 1. use_dense_linear_algebra = false
+    //
+    //    If the reciprocal_condition_number of J'J is less than
+    //    min_reciprocal_condition_number, Covariance::Compute() will
+    //    fail and return false.
+    //
+    // 2. use_dense_linear_algebra = true
+    //
+    //    When dense covariance estimation is being used, then rank
+    //    deficiency/singularity of the Jacobian can be handled in a
+    //    more sophisticated manner.
+    //
+    //    If null_space_rank = -1, then instead of computing the
+    //    inverse of J'J, the Moore-Penrose Pseudoinverse is computed. If
+    //    (lambda_i, e_i) are eigenvalue and eigenvector pairs of J'J.
+    //
+    //      pseudoinverse[J'J] = sum_i e_i e_i' / lambda_i
+    //
+    //    if lambda_i / lambda_max >= min_reciprocal_condition_number.
+    //
+    //    If null_space_rank is non-negative, then the smallest
+    //    null_space_rank eigenvalue/eigenvectors are dropped
+    //    irrespective of the magnitude of lambda_i. If the ratio of
+    //    the smallest non-zero eigenvalue to the largest eigenvalue
+    //    in the truncated matrix is still below
+    //    min_reciprocal_condition_number, then the
+    //    Covariance::Compute() will fail and return false.
+    double min_reciprocal_condition_number;
 
-    // When use_dense_linear_algebra is true, the bottom
-    // null_space_rank singular values are set to zero.
+    // When use_dense_linear_algebra is true, null_space_rank
+    // determines how many of the smallest eigenvectors of J'J are
+    // dropped when computing the pseudoinverse.
+    //
+    // If null_space_rank = -1, then instead of computing the inverse
+    // of J'J, the Moore-Penrose Pseudoinverse is computed. If
+    // (lambda_i, e_i) are eigenvalue and eigenvector pairs of J'J.
+    //
+    //   pseudoinverse[J'J] = sum_i e_i e_i' / lambda_i
+    //
+    //   if lambda_i / lambda_max >= min_reciprocal_condition_number.
+    //
+    // If null_space_rank is non-negative, then the smallest
+    // null_space_rank eigenvalue/eigenvectors are dropped
+    // irrespective of the magnitude of lambda_i. If the ratio of the
+    // smallest non-zero eigenvalue to the largest eigenvalue in the
+    // truncated matrix is still below
+    // min_reciprocal_condition_number, then the Covariance::Compute()
+    // will fail and return false.
     int null_space_rank;
 
     // Even though the residual blocks in the problem may contain loss
@@ -254,20 +327,28 @@
   // what parts of the covariance matrix are computed. The full
   // Jacobian is used to do the computation, i.e. they do not have an
   // impact on what part of the Jacobian is used for computation.
+  //
+  // The return value indicates the success or failure of the
+  // covariance computation. Please see the documentation for
+  // Covariance::Options for more on the conditions under which this
+  // function returns false.
   bool Compute(
       const vector<pair<const double*, const double*> >& covariance_blocks,
       Problem* problem);
 
-  // Compute must be called before the first call to
-  // GetCovarianceBlock. covariance_block must point to a memory
-  // location that can store a parameter_block1_size x
-  // parameter_block2_size matrix. The returned covariance will be a
-  // row-major matrix.
+  // Return the block of the covariance matrix corresponding to
+  // parameter_block1 and parameter_block2.
   //
-  // The pair <parameter_block1, parameter_block2> OR the pair
-  // <parameter_block2, parameter_block1> must have been present in
-  // the vector covariance_blocks when Compute was called. Otherwise
+  // Compute must be called before the first call to
+  // GetCovarianceBlock and the pair <parameter_block1,
+  // parameter_block2> OR the pair <parameter_block2,
+  // parameter_block1> must have been present in the vector
+  // covariance_blocks when Compute was called. Otherwise
   // GetCovarianceBlock will return false.
+  //
+  // covariance_block must point to a memory location that can store a
+  // parameter_block1_size x parameter_block2_size matrix. The
+  // returned covariance will be a row-major matrix.
   bool GetCovarianceBlock(const double* parameter_block1,
                           const double* parameter_block2,
                           double* covariance_block) const;
diff --git a/internal/ceres/covariance_impl.cc b/internal/ceres/covariance_impl.cc
index 1a96193..966ba35 100644
--- a/internal/ceres/covariance_impl.cc
+++ b/internal/ceres/covariance_impl.cc
@@ -397,9 +397,23 @@
 
   cholmod_factor* factor = ss_.AnalyzeCholesky(&cholmod_jacobian_view);
   event_logger.AddEvent("Symbolic Factorization");
-  bool status = 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());
+    if (reciprocal_condition_number <
+        options_.min_reciprocal_condition_number) {
+      LOG(WARNING) << "Cholesky factorization of J'J is not reliable. "
+                   << "Reciprocal condition number: "
+                   << reciprocal_condition_number << " "
+                   << "min_reciprocal_condition_number : "
+                   << options_.min_reciprocal_condition_number;
+      factorization_succeeded = false;
+    }
+  }
+
   event_logger.AddEvent("Numeric Factorization");
-  if (!status) {
+  if (!factorization_succeeded) {
     ss_.Free(factor);
     LOG(WARNING) << "Cholesky factorization failed.";
     return false;
@@ -520,22 +534,39 @@
 
   Eigen::JacobiSVD<Matrix> svd(dense_jacobian,
                                Eigen::ComputeThinU | Eigen::ComputeThinV);
-  Vector inverse_singular_values = svd.singularValues();
-  event_logger.AddEvent("SVD");
+  event_logger.AddEvent("SingularValueDecomposition");
 
-  for (int i = 0; i < inverse_singular_values.rows(); ++i) {
-    if (inverse_singular_values[i] > options_.min_singular_value_threshold &&
-        i < (inverse_singular_values.rows() - options_.null_space_rank)) {
-      inverse_singular_values[i] =
-          1.0 / (inverse_singular_values[i] * inverse_singular_values[i]);
-    } else {
-      inverse_singular_values[i] = 0.0;
+  const Vector singular_values = svd.singularValues();
+  const int num_singular_values = singular_values.rows();
+  Vector inverse_squared_singular_values(num_singular_values);
+  inverse_squared_singular_values.setZero();
+
+  const double max_singular_value = singular_values[0];
+  const double min_singular_value_ratio =
+      sqrt(options_.min_reciprocal_condition_number);
+
+  const bool automatic_truncation = (options_.null_space_rank < 0);
+  const int max_rank = min(num_singular_values,
+                           num_singular_values - options_.null_space_rank);
+
+  for (int i = 0; i < max_rank; ++i) {
+    const double singular_value_ratio = singular_values[i] / max_singular_value;
+    if (singular_value_ratio >= min_singular_value_ratio) {
+      inverse_squared_singular_values[i] =
+          1.0 / (singular_values[i] * singular_values[i]);
+    } else if (!automatic_truncation) {
+      LOG(WARNING) << "Cholesky factorization of J'J is not reliable. "
+                   << "Reciprocal condition number: "
+                   << singular_value_ratio * singular_value_ratio << " "
+                   << "min_reciprocal_condition_number : "
+                   << options_.min_reciprocal_condition_number;
+      return false;
     }
   }
 
   Matrix dense_covariance =
       svd.matrixV() *
-      inverse_singular_values.asDiagonal() *
+      inverse_squared_singular_values.asDiagonal() *
       svd.matrixV().transpose();
   event_logger.AddEvent("PseudoInverse");
 
diff --git a/internal/ceres/covariance_test.cc b/internal/ceres/covariance_test.cc
index b74efda..05fb234 100644
--- a/internal/ceres/covariance_test.cc
+++ b/internal/ceres/covariance_test.cc
@@ -532,15 +532,24 @@
     -1.6076e-02,   1.2549e-02,  -3.3329e-04,  -6.6659e-04,  -9.9988e-04,   3.9539e-02
   };
 
-  Covariance::Options options;
-  options.use_dense_linear_algebra = true;
-  options.null_space_rank = 1;
-  ComputeAndCompareCovarianceBlocks(options, expected_covariance);
 
-  options.use_dense_linear_algebra = true;
-  options.null_space_rank = 0;
-  options.min_singular_value_threshold = sqrt(3.5);
-  ComputeAndCompareCovarianceBlocks(options, expected_covariance);
+  {
+    Covariance::Options options;
+    options.use_dense_linear_algebra = true;
+    // Force dropping of the smallest eigenvector.
+    options.null_space_rank = 1;
+    ComputeAndCompareCovarianceBlocks(options, expected_covariance);
+  }
+
+  {
+    Covariance::Options options;
+    options.use_dense_linear_algebra = true;
+    // Force dropping of the smallest eigenvector via the ratio but
+    // automatic truncation.
+    options.min_reciprocal_condition_number = 0.044494;
+    options.null_space_rank = -1;
+    ComputeAndCompareCovarianceBlocks(options, expected_covariance);
+  }
 }
 
 class RankDeficientCovarianceTest : public CovarianceTest {
@@ -598,7 +607,7 @@
   }
 };
 
-TEST_F(RankDeficientCovarianceTest, MinSingularValueTolerance) {
+TEST_F(RankDeficientCovarianceTest, AutomaticTruncation) {
   // J
   //
   //   1  0  0  0  0  0
@@ -631,15 +640,7 @@
 
   Covariance::Options options;
   options.use_dense_linear_algebra = true;
-  ComputeAndCompareCovarianceBlocks(options, expected_covariance);
-
-  options.null_space_rank = 1;
-  ComputeAndCompareCovarianceBlocks(options, expected_covariance);
-
-  options.null_space_rank = 2;
-  ComputeAndCompareCovarianceBlocks(options, expected_covariance);
-
-  options.null_space_rank = 3;
+  options.null_space_rank = -1;
   ComputeAndCompareCovarianceBlocks(options, expected_covariance);
 }