diff --git a/albatross/covariance_functions/distance_metrics.h b/albatross/covariance_functions/distance_metrics.h
index c9afb8ee..21a13156 100644
--- a/albatross/covariance_functions/distance_metrics.h
+++ b/albatross/covariance_functions/distance_metrics.h
@@ -18,6 +18,8 @@
 
 namespace albatross {
 
+constexpr double EPSILON = 1e-16;
+
 class DistanceMetric : public ParameterHandlingMixin {
 public:
   DistanceMetric(){};
@@ -45,6 +47,12 @@ class EuclideanDistance : public DistanceMetric {
   }
 };
 
+template <typename _Scalar, int _Rows>
+double radial_distance(const Eigen::Matrix<_Scalar, _Rows, 1> &x,
+                       const Eigen::Matrix<_Scalar, _Rows, 1> &y) {
+  return fabs(x.norm() - y.norm());
+}
+
 class RadialDistance : public DistanceMetric {
 public:
   ~RadialDistance(){};
@@ -52,10 +60,25 @@ class RadialDistance : public DistanceMetric {
   std::string get_name() const override { return "radial_distance"; };
 
   double operator()(const Eigen::VectorXd &x, const Eigen::VectorXd &y) const {
-    return fabs(x.norm() - y.norm());
+    return radial_distance(x, y);
   }
 };
 
+template <typename _Scalar, int _Rows>
+double angular_distance(const Eigen::Matrix<_Scalar, _Rows, 1> &x,
+                        const Eigen::Matrix<_Scalar, _Rows, 1> &y) {
+  // The acos operator doesn't behave well near |1|.  acos(1.), for example,
+  // returns NaN, so here we do some special casing,
+  double dot_product = x.dot(y) / (x.norm() * y.norm());
+  if (dot_product > 1. - EPSILON) {
+    return 0.;
+  } else if (dot_product < -1. + EPSILON) {
+    return M_PI;
+  } else {
+    return acos(dot_product);
+  }
+}
+
 class AngularDistance : public DistanceMetric {
 public:
   ~AngularDistance(){};
@@ -65,16 +88,7 @@ class AngularDistance : public DistanceMetric {
   template <typename _Scalar, int _Rows>
   double operator()(const Eigen::Matrix<_Scalar, _Rows, 1> &x,
                     const Eigen::Matrix<_Scalar, _Rows, 1> &y) const {
-    // The acos operator doesn't behave well near |1|.  acos(1.), for example,
-    // returns NaN, so here we do some special casing,
-    double dot_product = x.dot(y) / (x.norm() * y.norm());
-    if (dot_product > 1. - 1e-16) {
-      return 0.;
-    } else if (dot_product < -1. + 1e-16) {
-      return M_PI;
-    } else {
-      return acos(dot_product);
-    }
+    return angular_distance(x, y);
   }
 };
 
diff --git a/albatross/covariance_functions/radial.h b/albatross/covariance_functions/radial.h
index 0f96bc77..69a31e8d 100644
--- a/albatross/covariance_functions/radial.h
+++ b/albatross/covariance_functions/radial.h
@@ -50,6 +50,12 @@ class RadialCovariance : public CovarianceTerm {
   DistanceMetricImpl distance_metric_;
 };
 
+inline double squared_exponential_covariance(double distance,
+                                             double length_scale,
+                                             double sigma = 1.) {
+  return sigma * sigma * exp(-pow(distance / length_scale, 2));
+}
+
 /*
  * SquaredExponential distance
  *  - c(d) = -exp((d/length_scale)^2)
@@ -87,12 +93,16 @@ class SquaredExponential : public RadialCovariance<DistanceMetricImpl> {
   double operator()(const X &x, const X &y) const {
     double distance = this->distance_metric_(x, y);
     double length_scale = this->params_.at("squared_exponential_length_scale");
-    distance /= length_scale;
     double sigma = this->params_.at("sigma_squared_exponential");
-    return sigma * sigma * exp(-distance * distance);
+    return squared_exponential_covariance(distance, length_scale, sigma);
   }
 };
 
+inline double exponential_covariance(double distance, double length_scale,
+                                     double sigma = 1.) {
+  return sigma * sigma * exp(-fabs(distance / length_scale));
+}
+
 /*
  * Exponential distance
  *  - c(d) = -exp(|d|/length_scale)
@@ -121,9 +131,8 @@ class Exponential : public RadialCovariance<DistanceMetricImpl> {
   double operator()(const X &x, const X &y) const {
     double distance = this->distance_metric_(x, y);
     double length_scale = this->params_.at("exponential_length_scale");
-    distance /= length_scale;
     double sigma = this->params_.at("sigma_exponential");
-    return sigma * sigma * exp(-fabs(distance));
+    return exponential_covariance(distance, length_scale, sigma);
   }
 };
 
diff --git a/albatross/models/gp.h b/albatross/models/gp.h
index 77f35759..184cd401 100644
--- a/albatross/models/gp.h
+++ b/albatross/models/gp.h
@@ -147,8 +147,6 @@ class GaussianProcessRegression
         covariance_function_, features, this->model_fit_.train_features);
     // Then we can use the information vector to determine the posterior
     const Eigen::VectorXd pred = cross_cov * this->model_fit_.information;
-    // TODO: right now this is recomputing the LDLT, so is highly inefficient,
-    // Ideally this would get stored inside GaussianProcessFit.
     Eigen::MatrixXd pred_cov =
         symmetric_covariance(covariance_function_, features);
     auto ldlt = this->model_fit_.train_ldlt;