diff --git a/GPy/inference/latent_function_inference/var_gauss.py b/GPy/inference/latent_function_inference/var_gauss.py
index e71416b4..bef1cce6 100644
--- a/GPy/inference/latent_function_inference/var_gauss.py
+++ b/GPy/inference/latent_function_inference/var_gauss.py
@@ -47,7 +47,10 @@ class VarGauss(LatentFunctionInference):
             dL_dthetaL = np.array([])
         dF_da = np.dot(K, dF_dm)
         SigmaB = Sigma*self.beta
-        dF_db = -np.diag(Sigma.dot(np.diag(dF_dv.flatten())).dot(SigmaB))*2
+        #dF_db_ = -np.diag(Sigma.dot(np.diag(dF_dv.flatten())).dot(SigmaB))*2
+        dF_db = -2*np.sum(Sigma**2 * (dF_dv * self.beta), 0)
+        #assert np.allclose(dF_db, dF_db_)
+
         KL = 0.5*(Alogdet + np.trace(Ai) - num_data + np.sum(m*self.alpha))
         dKL_da = m
         A_A2 = Ai - Ai.dot(Ai)
diff --git a/GPy/kern/_src/psi_comp/linear_psi_comp.py b/GPy/kern/_src/psi_comp/linear_psi_comp.py
index f3262d3b..87ceaf66 100644
--- a/GPy/kern/_src/psi_comp/linear_psi_comp.py
+++ b/GPy/kern/_src/psi_comp/linear_psi_comp.py
@@ -21,11 +21,12 @@ def psicomputations(variance, Z, variational_posterior, return_psi2_n=False):
     S = variational_posterior.variance
 
     psi0 = (variance*(np.square(mu)+S)).sum(axis=1)
-    psi1 = np.dot(mu,(variance*Z).T)
-    if not return_psi2_n:
-        psi2 = np.dot(S.sum(axis=0)*np.square(variance)*Z,Z.T)+ tdot(psi1.T)
+    Zv = variance * Z
+    psi1 = np.dot(mu,Zv.T)
+    if return_psi2_n:
+        psi2 = psi1[:,:,None] * psi1[:,None,:] + np.dot(S[:,None,:] * Zv[None,:,:], Zv.T)
     else:
-        raise NotImplementedError
+        psi2 = np.dot(S.sum(axis=0) * Zv, Zv.T) + tdot(psi1.T)
 
     return psi0, psi1, psi2