diff --git a/internal/ceres/covariance_test.cc b/internal/ceres/covariance_test.cc
index dea0723..229173f 100644
--- a/internal/ceres/covariance_test.cc
+++ b/internal/ceres/covariance_test.cc
@@ -31,12 +31,13 @@
 #include "ceres/covariance.h"
 
 #include <algorithm>
-#include <cstdint>
 #include <cmath>
+#include <cstdint>
 #include <map>
 #include <memory>
 #include <utility>
 
+#include "ceres/autodiff_cost_function.h"
 #include "ceres/compressed_row_sparse_matrix.h"
 #include "ceres/cost_function.h"
 #include "ceres/covariance_impl.h"
@@ -53,7 +54,7 @@
 using std::pair;
 using std::vector;
 
-class UnaryCostFunction: public CostFunction {
+class UnaryCostFunction : public CostFunction {
  public:
   UnaryCostFunction(const int num_residuals,
                     const int32_t parameter_block_size,
@@ -63,9 +64,9 @@
     mutable_parameter_block_sizes()->push_back(parameter_block_size);
   }
 
-  virtual bool Evaluate(double const* const* parameters,
-                        double* residuals,
-                        double** jacobians) const {
+  bool Evaluate(double const* const* parameters,
+                double* residuals,
+                double** jacobians) const final {
     for (int i = 0; i < num_residuals(); ++i) {
       residuals[i] = 1;
     }
@@ -85,8 +86,7 @@
   vector<double> jacobian_;
 };
 
-
-class BinaryCostFunction: public CostFunction {
+class BinaryCostFunction : public CostFunction {
  public:
   BinaryCostFunction(const int num_residuals,
                      const int32_t parameter_block1_size,
@@ -102,9 +102,9 @@
     mutable_parameter_block_sizes()->push_back(parameter_block2_size);
   }
 
-  virtual bool Evaluate(double const* const* parameters,
-                        double* residuals,
-                        double** jacobians) const {
+  bool Evaluate(double const* const* parameters,
+                double* residuals,
+                double** jacobians) const final {
     for (int i = 0; i < num_residuals(); ++i) {
       residuals[i] = 2;
     }
@@ -134,22 +134,22 @@
  public:
   virtual ~PolynomialParameterization() {}
 
-  virtual bool Plus(const double* x,
-                    const double* delta,
-                    double* x_plus_delta) const {
+  bool Plus(const double* x,
+            const double* delta,
+            double* x_plus_delta) const final {
     x_plus_delta[0] = delta[0] * x[0];
     x_plus_delta[1] = delta[0] * x[1];
     return true;
   }
 
-  virtual bool ComputeJacobian(const double* x, double* jacobian) const {
+  bool ComputeJacobian(const double* x, double* jacobian) const final {
     jacobian[0] = x[0];
     jacobian[1] = x[1];
     return true;
   }
 
-  virtual int GlobalSize() const { return 2; }
-  virtual int LocalSize() const { return 1; }
+  int GlobalSize() const final { return 2; }
+  int LocalSize() const final { return 1; }
 };
 
 TEST(CovarianceImpl, ComputeCovarianceSparsity) {
@@ -192,6 +192,7 @@
   //  . . . . . . X X X X
   //  . . . . . . X X X X
 
+  // clang-format off
   int expected_rows[] = {0, 5, 10, 15, 18, 21, 24, 28, 32, 36, 40};
   int expected_cols[] = {0, 6, 7, 8, 9,
                          1, 2, 3, 4, 5,
@@ -203,7 +204,7 @@
                          6, 7, 8, 9,
                          6, 7, 8, 9,
                          6, 7, 8, 9};
-
+  // clang-format on
 
   vector<pair<const double*, const double*>> covariance_blocks;
   covariance_blocks.push_back(make_pair(block1, block1));
@@ -215,8 +216,8 @@
 
   Covariance::Options options;
   CovarianceImpl covariance_impl(options);
-  EXPECT_TRUE(covariance_impl
-              .ComputeCovarianceSparsity(covariance_blocks, &problem));
+  EXPECT_TRUE(
+      covariance_impl.ComputeCovarianceSparsity(covariance_blocks, &problem));
 
   const CompressedRowSparseMatrix* crsm = covariance_impl.covariance_matrix();
 
@@ -227,17 +228,13 @@
   const int* rows = crsm->rows();
   for (int r = 0; r < crsm->num_rows() + 1; ++r) {
     EXPECT_EQ(rows[r], expected_rows[r])
-        << r << " "
-        << rows[r] << " "
-        << expected_rows[r];
+        << r << " " << rows[r] << " " << expected_rows[r];
   }
 
   const int* cols = crsm->cols();
   for (int c = 0; c < crsm->num_nonzeros(); ++c) {
     EXPECT_EQ(cols[c], expected_cols[c])
-        << c << " "
-        << cols[c] << " "
-        << expected_cols[c];
+        << c << " " << cols[c] << " " << expected_cols[c];
   }
 }
 
@@ -279,6 +276,7 @@
   //  . . . X X X X
   //  . . . X X X X
 
+  // clang-format off
   int expected_rows[] = {0, 5, 7, 9, 13, 17, 21, 25};
   int expected_cols[] = {0, 3, 4, 5, 6,
                          1, 2,
@@ -287,6 +285,7 @@
                          3, 4, 5, 6,
                          3, 4, 5, 6,
                          3, 4, 5, 6};
+  // clang-format on
 
   vector<pair<const double*, const double*>> covariance_blocks;
   covariance_blocks.push_back(make_pair(block1, block1));
@@ -298,8 +297,8 @@
 
   Covariance::Options options;
   CovarianceImpl covariance_impl(options);
-  EXPECT_TRUE(covariance_impl
-              .ComputeCovarianceSparsity(covariance_blocks, &problem));
+  EXPECT_TRUE(
+      covariance_impl.ComputeCovarianceSparsity(covariance_blocks, &problem));
 
   const CompressedRowSparseMatrix* crsm = covariance_impl.covariance_matrix();
 
@@ -310,17 +309,13 @@
   const int* rows = crsm->rows();
   for (int r = 0; r < crsm->num_rows() + 1; ++r) {
     EXPECT_EQ(rows[r], expected_rows[r])
-        << r << " "
-        << rows[r] << " "
-        << expected_rows[r];
+        << r << " " << rows[r] << " " << expected_rows[r];
   }
 
   const int* cols = crsm->cols();
   for (int c = 0; c < crsm->num_nonzeros(); ++c) {
     EXPECT_EQ(cols[c], expected_cols[c])
-        << c << " "
-        << cols[c] << " "
-        << expected_cols[c];
+        << c << " " << cols[c] << " " << expected_cols[c];
   }
 }
 
@@ -360,6 +355,7 @@
   //  . . . X X X X
   //  . . . X X X X
 
+  // clang-format off
   int expected_rows[] = {0, 5, 7, 9, 13, 17, 21, 25};
   int expected_cols[] = {0, 3, 4, 5, 6,
                          1, 2,
@@ -368,6 +364,7 @@
                          3, 4, 5, 6,
                          3, 4, 5, 6,
                          3, 4, 5, 6};
+  // clang-format on
 
   vector<pair<const double*, const double*>> covariance_blocks;
   covariance_blocks.push_back(make_pair(block1, block1));
@@ -379,8 +376,8 @@
 
   Covariance::Options options;
   CovarianceImpl covariance_impl(options);
-  EXPECT_TRUE(covariance_impl
-              .ComputeCovarianceSparsity(covariance_blocks, &problem));
+  EXPECT_TRUE(
+      covariance_impl.ComputeCovarianceSparsity(covariance_blocks, &problem));
 
   const CompressedRowSparseMatrix* crsm = covariance_impl.covariance_matrix();
 
@@ -391,17 +388,13 @@
   const int* rows = crsm->rows();
   for (int r = 0; r < crsm->num_rows() + 1; ++r) {
     EXPECT_EQ(rows[r], expected_rows[r])
-        << r << " "
-        << rows[r] << " "
-        << expected_rows[r];
+        << r << " " << rows[r] << " " << expected_rows[r];
   }
 
   const int* cols = crsm->cols();
   for (int c = 0; c < crsm->num_nonzeros(); ++c) {
     EXPECT_EQ(cols[c], expected_cols[c])
-        << c << " "
-        << cols[c] << " "
-        << expected_cols[c];
+        << c << " " << cols[c] << " " << expected_cols[c];
   }
 }
 
@@ -409,7 +402,7 @@
  protected:
   typedef map<const double*, pair<int, int>> BoundsMap;
 
-  virtual void SetUp() {
+  void SetUp() override {
     double* x = parameters_;
     double* y = x + 2;
     double* z = y + 3;
@@ -422,40 +415,33 @@
     z[0] = 3;
 
     {
-      double jacobian[] = { 1.0, 0.0, 0.0, 1.0};
+      double jacobian[] = {1.0, 0.0, 0.0, 1.0};
       problem_.AddResidualBlock(new UnaryCostFunction(2, 2, jacobian), NULL, x);
     }
 
     {
-      double jacobian[] = { 2.0, 0.0, 0.0, 0.0, 2.0, 0.0, 0.0, 0.0, 2.0 };
+      double jacobian[] = {2.0, 0.0, 0.0, 0.0, 2.0, 0.0, 0.0, 0.0, 2.0};
       problem_.AddResidualBlock(new UnaryCostFunction(3, 3, jacobian), NULL, y);
     }
 
     {
       double jacobian = 5.0;
-      problem_.AddResidualBlock(new UnaryCostFunction(1, 1, &jacobian),
-                                NULL,
-                                z);
+      problem_.AddResidualBlock(
+          new UnaryCostFunction(1, 1, &jacobian), NULL, z);
     }
 
     {
-      double jacobian1[] = { 1.0, 2.0, 3.0 };
-      double jacobian2[] = { -5.0, -6.0 };
+      double jacobian1[] = {1.0, 2.0, 3.0};
+      double jacobian2[] = {-5.0, -6.0};
       problem_.AddResidualBlock(
-          new BinaryCostFunction(1, 3, 2, jacobian1, jacobian2),
-          NULL,
-          y,
-          x);
+          new BinaryCostFunction(1, 3, 2, jacobian1, jacobian2), NULL, y, x);
     }
 
     {
-      double jacobian1[] = {2.0 };
-      double jacobian2[] = { 3.0, -2.0 };
+      double jacobian1[] = {2.0};
+      double jacobian2[] = {3.0, -2.0};
       problem_.AddResidualBlock(
-          new BinaryCostFunction(1, 1, 2, jacobian1, jacobian2),
-          NULL,
-          z,
-          x);
+          new BinaryCostFunction(1, 1, 2, jacobian1, jacobian2), NULL, z, x);
     }
 
     all_covariance_blocks_.push_back(make_pair(x, x));
@@ -481,8 +467,7 @@
 
   // Computes covariance in tangent space.
   void ComputeAndCompareCovarianceBlocksInTangentSpace(
-                                         const Covariance::Options& options,
-                                         const double* expected_covariance) {
+      const Covariance::Options& options, const double* expected_covariance) {
     ComputeAndCompareCovarianceBlocksInTangentOrAmbientSpace(
         options,
         false,  // tangent
@@ -548,8 +533,9 @@
                                     bool lift_covariance_to_ambient_space,
                                     const Covariance& covariance,
                                     const double* expected_covariance) {
-    const BoundsMap& column_bounds = lift_covariance_to_ambient_space ?
-        column_bounds_ : local_column_bounds_;
+    const BoundsMap& column_bounds = lift_covariance_to_ambient_space
+                                         ? column_bounds_
+                                         : local_column_bounds_;
     const int row_begin = FindOrDie(column_bounds, block1).first;
     const int row_end = FindOrDie(column_bounds, block1).second;
     const int col_begin = FindOrDie(column_bounds, block2).first;
@@ -557,13 +543,10 @@
 
     Matrix actual(row_end - row_begin, col_end - col_begin);
     if (lift_covariance_to_ambient_space) {
-      EXPECT_TRUE(covariance.GetCovarianceBlock(block1,
-                                                block2,
-                                                actual.data()));
+      EXPECT_TRUE(covariance.GetCovarianceBlock(block1, block2, actual.data()));
     } else {
-      EXPECT_TRUE(covariance.GetCovarianceBlockInTangentSpace(block1,
-                                                              block2,
-                                                              actual.data()));
+      EXPECT_TRUE(covariance.GetCovarianceBlockInTangentSpace(
+          block1, block2, actual.data()));
     }
 
     int dof = 0;  // degrees of freedom = sum of LocalSize()s
@@ -571,22 +554,22 @@
       dof = std::max(dof, bound.second.second);
     }
     ConstMatrixRef expected(expected_covariance, dof, dof);
-    double diff_norm = (expected.block(row_begin,
-                                       col_begin,
-                                       row_end - row_begin,
-                                       col_end - col_begin) - actual).norm();
+    double diff_norm =
+        (expected.block(
+             row_begin, col_begin, row_end - row_begin, col_end - col_begin) -
+         actual)
+            .norm();
     diff_norm /= (row_end - row_begin) * (col_end - col_begin);
 
     const double kTolerance = 1e-5;
     EXPECT_NEAR(diff_norm, 0.0, kTolerance)
         << "rows: " << row_begin << " " << row_end << "  "
         << "cols: " << col_begin << " " << col_end << "  "
-        << "\n\n expected: \n " << expected.block(row_begin,
-                                                  col_begin,
-                                                  row_end - row_begin,
-                                                  col_end - col_begin)
-        << "\n\n actual: \n " << actual
-        << "\n\n full expected: \n" << expected;
+        << "\n\n expected: \n "
+        << expected.block(
+               row_begin, col_begin, row_end - row_begin, col_end - col_begin)
+        << "\n\n actual: \n " << actual << "\n\n full expected: \n"
+        << expected;
   }
 
   double parameters_[6];
@@ -596,7 +579,6 @@
   BoundsMap local_column_bounds_;
 };
 
-
 TEST_F(CovarianceTest, NormalBehavior) {
   // J
   //
@@ -619,6 +601,7 @@
   //    6  -4  0   0   0 29
 
   // inv(J'J) computed using octave.
+  // clang-format off
   double expected_covariance[] = {
      7.0747e-02,  -8.4923e-03,   1.6821e-02,   3.3643e-02,   5.0464e-02,  -1.5809e-02,  // NOLINT
     -8.4923e-03,   8.1352e-02,   2.4758e-02,   4.9517e-02,   7.4275e-02,   1.2978e-02,  // NOLINT
@@ -627,6 +610,7 @@
      5.0464e-02,   7.4275e-02,  -2.8906e-03,  -5.7813e-03,   2.4133e-01,  -1.9598e-04,  // NOLINT
     -1.5809e-02,   1.2978e-02,  -6.5325e-05,  -1.3065e-04,  -1.9598e-04,   3.9544e-02,  // NOLINT
   };
+  // clang-format on
 
   Covariance::Options options;
 
@@ -668,6 +652,7 @@
   //    6  -4  0   0   0 29
 
   // inv(J'J) computed using octave.
+  // clang-format off
   double expected_covariance[] = {
      7.0747e-02,  -8.4923e-03,   1.6821e-02,   3.3643e-02,   5.0464e-02,  -1.5809e-02,  // NOLINT
     -8.4923e-03,   8.1352e-02,   2.4758e-02,   4.9517e-02,   7.4275e-02,   1.2978e-02,  // NOLINT
@@ -676,6 +661,7 @@
      5.0464e-02,   7.4275e-02,  -2.8906e-03,  -5.7813e-03,   2.4133e-01,  -1.9598e-04,  // NOLINT
     -1.5809e-02,   1.2978e-02,  -6.5325e-05,  -1.3065e-04,  -1.9598e-04,   3.9544e-02,  // NOLINT
   };
+  // clang-format on
 
   Covariance::Options options;
   options.num_threads = 4;
@@ -720,6 +706,7 @@
   //  0  0  0  0  0 29
 
   // pinv(J'J) computed using octave.
+  // clang-format off
   double expected_covariance[] = {
               0,            0,            0,            0,            0,            0,  // NOLINT
               0,            0,            0,            0,            0,            0,  // NOLINT
@@ -727,6 +714,7 @@
               0,            0,     -0.02778,      0.19444,     -0.08333,     -0.00000,  // NOLINT
               0,            0,     -0.04167,     -0.08333,      0.12500,     -0.00000,  // NOLINT
               0,            0,     -0.00000,     -0.00000,     -0.00000,      0.03448   // NOLINT
+      // clang-format on
   };
 
   Covariance::Options options;
@@ -777,6 +765,7 @@
 
   // A * inv((J*A)'*(J*A)) * A'
   // Computed using octave.
+  // clang-format off
   double expected_covariance[] = {
     0.01766,   0.01766,   0.02158,   0.04316,   0.00000,  -0.00122,
     0.01766,   0.01766,   0.02158,   0.04316,   0.00000,  -0.00122,
@@ -785,6 +774,7 @@
     0.00000,   0.00000,   0.00000,   0.00000,   0.00000,   0.00000,
    -0.00122,  -0.00122,  -0.00149,  -0.00298,   0.00000,   0.03457
   };
+  // clang-format on
 
   Covariance::Options options;
 
@@ -839,12 +829,14 @@
 
   // inv((J*A)'*(J*A))
   // Computed using octave.
+  // clang-format off
   double expected_covariance[] = {
     0.01766,   0.02158,   0.04316,   -0.00122,
     0.02158,   0.24860,  -0.00281,   -0.00149,
     0.04316,  -0.00281,   0.24439,   -0.00298,
    -0.00122,  -0.00149,  -0.00298,    0.03457  // NOLINT
   };
+  // clang-format on
 
   Covariance::Options options;
 
@@ -902,12 +894,14 @@
 
   // pinv((J*A)'*(J*A))
   // Computed using octave.
+  // clang-format off
   double expected_covariance[] = {
     0.0, 0.0, 0.0, 0.0,
     0.0, 0.0, 0.0, 0.0,
     0.0, 0.0, 0.0, 0.0,
     0.0, 0.0, 0.0, 0.034482 // NOLINT
   };
+  // clang-format on
 
   Covariance::Options options;
 
@@ -950,6 +944,7 @@
   // 3.4142 is the smallest eigen value of J'J. The following matrix
   // was obtained by dropping the eigenvector corresponding to this
   // eigenvalue.
+  // clang-format off
   double expected_covariance[] = {
      5.4135e-02,  -3.5121e-02,   1.7257e-04,   3.4514e-04,   5.1771e-04,  -1.6076e-02,  // NOLINT
     -3.5121e-02,   3.8667e-02,  -1.9288e-03,  -3.8576e-03,  -5.7864e-03,   1.2549e-02,  // NOLINT
@@ -958,7 +953,7 @@
      5.1771e-04,  -5.7864e-03,  -5.2946e-02,  -1.0589e-01,   9.1162e-02,  -9.9988e-04,  // NOLINT
     -1.6076e-02,   1.2549e-02,  -3.3329e-04,  -6.6659e-04,  -9.9988e-04,   3.9539e-02   // NOLINT
   };
-
+  // clang-format on
 
   {
     Covariance::Options options;
@@ -1102,46 +1097,39 @@
 
 class RankDeficientCovarianceTest : public CovarianceTest {
  protected:
-  virtual void SetUp() {
+  void SetUp() final {
     double* x = parameters_;
     double* y = x + 2;
     double* z = y + 3;
 
     {
-      double jacobian[] = { 1.0, 0.0, 0.0, 1.0};
+      double jacobian[] = {1.0, 0.0, 0.0, 1.0};
       problem_.AddResidualBlock(new UnaryCostFunction(2, 2, jacobian), NULL, x);
     }
 
     {
-      double jacobian[] = { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 };
+      double jacobian[] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
       problem_.AddResidualBlock(new UnaryCostFunction(3, 3, jacobian), NULL, y);
     }
 
     {
       double jacobian = 5.0;
-      problem_.AddResidualBlock(new UnaryCostFunction(1, 1, &jacobian),
-                                NULL,
-                                z);
+      problem_.AddResidualBlock(
+          new UnaryCostFunction(1, 1, &jacobian), NULL, z);
     }
 
     {
-      double jacobian1[] = { 0.0, 0.0, 0.0 };
-      double jacobian2[] = { -5.0, -6.0 };
+      double jacobian1[] = {0.0, 0.0, 0.0};
+      double jacobian2[] = {-5.0, -6.0};
       problem_.AddResidualBlock(
-          new BinaryCostFunction(1, 3, 2, jacobian1, jacobian2),
-          NULL,
-          y,
-          x);
+          new BinaryCostFunction(1, 3, 2, jacobian1, jacobian2), NULL, y, x);
     }
 
     {
-      double jacobian1[] = {2.0 };
-      double jacobian2[] = { 3.0, -2.0 };
+      double jacobian1[] = {2.0};
+      double jacobian2[] = {3.0, -2.0};
       problem_.AddResidualBlock(
-          new BinaryCostFunction(1, 1, 2, jacobian1, jacobian2),
-          NULL,
-          z,
-          x);
+          new BinaryCostFunction(1, 1, 2, jacobian1, jacobian2), NULL, z, x);
     }
 
     all_covariance_blocks_.push_back(make_pair(x, x));
@@ -1179,6 +1167,7 @@
   //   6 -4  0  0  0 29
 
   // pinv(J'J) computed using octave.
+  // clang-format off
   double expected_covariance[] = {
      0.053998,  -0.033145,   0.000000,   0.000000,   0.000000,  -0.015744,
     -0.033145,   0.045067,   0.000000,   0.000000,   0.000000,   0.013074,
@@ -1187,6 +1176,7 @@
      0.000000,   0.000000,   0.000000,   0.000000,   0.000000,   0.000000,
     -0.015744,   0.013074,   0.000000,   0.000000,   0.000000,   0.039543
   };
+  // clang-format on
 
   Covariance::Options options;
   options.algorithm_type = DENSE_SVD;
@@ -1194,9 +1184,90 @@
   ComputeAndCompareCovarianceBlocks(options, expected_covariance);
 }
 
+struct LinearCostFunction {
+  template <typename T>
+  bool operator()(const T* x, const T* y, T* residual) const {
+    residual[0] = T(10.0) - *x;
+    residual[1] = T(5.0) - *y;
+    return true;
+  }
+  static CostFunction* Create() {
+    return new AutoDiffCostFunction<LinearCostFunction, 2, 1, 1>(
+        new LinearCostFunction);
+  }
+};
+
+TEST(Covariance, ZeroSizedLocalParameterizationGetCovariance) {
+  double x = 0.0;
+  double y = 1.0;
+  Problem problem;
+  problem.AddResidualBlock(LinearCostFunction::Create(), nullptr, &x, &y);
+  problem.SetParameterization(&y, new SubsetParameterization(1, {0}));
+  // J = [-1 0]
+  //     [ 0 0]
+  Covariance::Options options;
+  options.algorithm_type = DENSE_SVD;
+  Covariance covariance(options);
+  vector<pair<const double*, const double*>> covariance_blocks;
+  covariance_blocks.push_back(std::make_pair(&x, &x));
+  covariance_blocks.push_back(std::make_pair(&x, &y));
+  covariance_blocks.push_back(std::make_pair(&y, &x));
+  covariance_blocks.push_back(std::make_pair(&y, &y));
+  EXPECT_TRUE(covariance.Compute(covariance_blocks, &problem));
+
+  double value = -1;
+  covariance.GetCovarianceBlock(&x, &x, &value);
+  EXPECT_NEAR(value, 1.0, std::numeric_limits<double>::epsilon());
+
+  value = -1;
+  covariance.GetCovarianceBlock(&x, &y, &value);
+  EXPECT_NEAR(value, 0.0, std::numeric_limits<double>::epsilon());
+
+  value = -1;
+  covariance.GetCovarianceBlock(&y, &x, &value);
+  EXPECT_NEAR(value, 0.0, std::numeric_limits<double>::epsilon());
+
+  value = -1;
+  covariance.GetCovarianceBlock(&y, &y, &value);
+  EXPECT_NEAR(value, 0.0, std::numeric_limits<double>::epsilon());
+}
+
+TEST(Covariance, ZeroSizedLocalParameterizationGetCovarianceInTangentSpace) {
+  double x = 0.0;
+  double y = 1.0;
+  Problem problem;
+  problem.AddResidualBlock(LinearCostFunction::Create(), nullptr, &x, &y);
+  problem.SetParameterization(&y, new SubsetParameterization(1, {0}));
+  // J = [-1 0]
+  //     [ 0 0]
+  Covariance::Options options;
+  options.algorithm_type = DENSE_SVD;
+  Covariance covariance(options);
+  vector<pair<const double*, const double*>> covariance_blocks;
+  covariance_blocks.push_back(std::make_pair(&x, &x));
+  covariance_blocks.push_back(std::make_pair(&x, &y));
+  covariance_blocks.push_back(std::make_pair(&y, &x));
+  covariance_blocks.push_back(std::make_pair(&y, &y));
+  EXPECT_TRUE(covariance.Compute(covariance_blocks, &problem));
+
+  double value = -1;
+  covariance.GetCovarianceBlockInTangentSpace(&x, &x, &value);
+  EXPECT_NEAR(value, 1.0, std::numeric_limits<double>::epsilon());
+
+  value = -1;
+  // The following three calls, should not touch this value, since the
+  // tangent space is of size zero
+  covariance.GetCovarianceBlockInTangentSpace(&x, &y, &value);
+  EXPECT_EQ(value, -1);
+  covariance.GetCovarianceBlockInTangentSpace(&y, &x, &value);
+  EXPECT_EQ(value, -1);
+  covariance.GetCovarianceBlockInTangentSpace(&y, &y, &value);
+  EXPECT_EQ(value, -1);
+}
+
 class LargeScaleCovarianceTest : public ::testing::Test {
  protected:
-  virtual void SetUp() {
+  void SetUp() final {
     num_parameter_blocks_ = 2000;
     parameter_block_size_ = 5;
     parameters_.reset(
@@ -1208,11 +1279,11 @@
       jacobian *= (i + 1);
 
       double* block_i = parameters_.get() + i * parameter_block_size_;
-      problem_.AddResidualBlock(new UnaryCostFunction(parameter_block_size_,
-                                                      parameter_block_size_,
-                                                      jacobian.data()),
-                                NULL,
-                                block_i);
+      problem_.AddResidualBlock(
+          new UnaryCostFunction(
+              parameter_block_size_, parameter_block_size_, jacobian.data()),
+          NULL,
+          block_i);
       for (int j = i; j < num_parameter_blocks_; ++j) {
         double* block_j = parameters_.get() + j * parameter_block_size_;
         all_covariance_blocks_.push_back(make_pair(block_i, block_j));
@@ -1244,8 +1315,10 @@
       covariance.GetCovarianceBlock(block_i, block_i, actual.data());
       EXPECT_NEAR((expected - actual).norm(), 0.0, kTolerance)
           << "block: " << i << ", " << i << "\n"
-          << "expected: \n" << expected << "\n"
-          << "actual: \n" << actual;
+          << "expected: \n"
+          << expected << "\n"
+          << "actual: \n"
+          << actual;
 
       expected.setZero();
       for (int j = i + 1; j < num_parameter_blocks_; ++j) {
@@ -1253,8 +1326,10 @@
         covariance.GetCovarianceBlock(block_i, block_j, actual.data());
         EXPECT_NEAR((expected - actual).norm(), 0.0, kTolerance)
             << "block: " << i << ", " << j << "\n"
-            << "expected: \n" << expected << "\n"
-            << "actual: \n" << actual;
+            << "expected: \n"
+            << expected << "\n"
+            << "actual: \n"
+            << actual;
       }
     }
   }