diff --git a/GPy/core/gp.py b/GPy/core/gp.py
index 38019fa7..ba81f096 100644
--- a/GPy/core/gp.py
+++ b/GPy/core/gp.py
@@ -42,10 +42,8 @@ class GP(Model):
         assert Y.shape[0] == self.num_data
         _, self.output_dim = self.Y.shape
 
-        if Y_metadata is None:
-            self.Y_metadata = {}
-        else:
-            self.Y_metadata = Y_metadata
+        #TODO: check the type of this is okay?
+        self.Y_metadata = Y_metadata
 
         assert isinstance(kernel, kern.Kern)
         #assert self.input_dim == kernel.input_dim
diff --git a/GPy/inference/latent_function_inference/exact_gaussian_inference.py b/GPy/inference/latent_function_inference/exact_gaussian_inference.py
index c4ed0a0f..95a15fcc 100644
--- a/GPy/inference/latent_function_inference/exact_gaussian_inference.py
+++ b/GPy/inference/latent_function_inference/exact_gaussian_inference.py
@@ -3,6 +3,7 @@
 
 from posterior import Posterior
 from ...util.linalg import pdinv, dpotrs, tdot
+from ...util import diag
 import numpy as np
 log_2_pi = np.log(2*np.pi)
 
@@ -41,7 +42,9 @@ class ExactGaussianInference(object):
 
         K = kern.K(X)
 
-        Wi, LW, LWi, W_logdet = pdinv(K + likelihood.covariance_matrix(Y, Y_metadata))
+        Ky = K.copy()
+        diag.add(Ky, likelihood.gaussian_variance(Y, Y_metadata))
+        Wi, LW, LWi, W_logdet = pdinv(Ky)
 
         alpha, _ = dpotrs(LW, YYT_factor, lower=1)
 
diff --git a/GPy/inference/latent_function_inference/var_dtc.py b/GPy/inference/latent_function_inference/var_dtc.py
index e2aa95f5..59672449 100644
--- a/GPy/inference/latent_function_inference/var_dtc.py
+++ b/GPy/inference/latent_function_inference/var_dtc.py
@@ -65,7 +65,7 @@ class VarDTC(object):
         _, output_dim = Y.shape
 
         #see whether we've got a different noise variance for each datum
-        beta = 1./np.fmax(likelihood.gaussian_variance(Y_metadata), 1e-6)
+        beta = 1./np.fmax(likelihood.gaussian_variance(Y, Y_metadata), 1e-6)
         # VVT_factor is a matrix such that tdot(VVT_factor) = VVT...this is for efficiency!
         #self.YYTfactor = self.get_YYTfactor(Y)
         #VVT_factor = self.get_VVTfactor(self.YYTfactor, beta)
diff --git a/GPy/likelihoods/gaussian.py b/GPy/likelihoods/gaussian.py
index 0c73e485..89ad4486 100644
--- a/GPy/likelihoods/gaussian.py
+++ b/GPy/likelihoods/gaussian.py
@@ -51,14 +51,12 @@ class Gaussian(Likelihood):
             self.log_concave = True
 
     def betaY(self,Y,Y_metadata=None):
+        #TODO: ~Ricardo this does not live here
         return Y/self.gaussian_variance(Y_metadata)
 
     def gaussian_variance(self, Y_metadata=None):
         return self.variance
 
-    def covariance_matrix(self, Y, Y_metadata=None):
-        return np.eye(Y.shape[0]) * self.variance
-
     def update_gradients(self, grad):
         self.variance.gradient = grad