diff --git a/GPy/core/gp.py b/GPy/core/gp.py
index 65295521..57e9715a 100644
--- a/GPy/core/gp.py
+++ b/GPy/core/gp.py
@@ -76,14 +76,17 @@ class GP(Model):
 
         """
         Kx = self.kern.K(_Xnew, self.X, which_parts=which_parts).T
-        LiKx, _ = dtrtrs(self.posterior.woodbury_chol, np.asfortranarray(Kx), lower=1)
+        #LiKx, _ = dtrtrs(self.posterior.woodbury_chol, np.asfortranarray(Kx), lower=1)
+        WiKx = np.dot(self.posterior.woodbury_inv, Kx)
         mu = np.dot(Kx.T, self.posterior.woodbury_vector)
         if full_cov:
             Kxx = self.kern.K(_Xnew, which_parts=which_parts)
-            var = Kxx - tdot(LiKx.T)
+            #var = Kxx - tdot(LiKx.T)
+            var = np.dot(Kx.T, WiKx)
         else:
             Kxx = self.kern.Kdiag(_Xnew, which_parts=which_parts)
-            var = Kxx - np.sum(LiKx*LiKx, 0)
+            #var = Kxx - np.sum(LiKx*LiKx, 0)
+            var = Kxx - np.sum(WiKx*Kx, 0)
             var = var.reshape(-1, 1)
         return mu, var
 
diff --git a/GPy/inference/latent_function_inference/laplace.py b/GPy/inference/latent_function_inference/laplace.py
index 9a66d3b6..6e252406 100644
--- a/GPy/inference/latent_function_inference/laplace.py
+++ b/GPy/inference/latent_function_inference/laplace.py
@@ -11,7 +11,7 @@
 #http://gaussianprocess.org/gpml/code.
 
 import numpy as np
-from ...util.linalg import mdot, jitchol, pddet, dpotrs, dtrtrs
+from ...util.linalg import mdot, jitchol, pddet, dpotrs, dtrtrs, dpotri, symmetrify
 from ...util.misc import param_to_array
 from functools import partial as partial_func
 from posterior import Posterior
@@ -216,8 +216,15 @@ class LaplaceInference(object):
         B = np.eye(K.shape[0]) + W_12*K*W_12.T
         L = jitchol(B)
 
-        LiW12, _ = dtrtrs(L, np.diag(W_12[:,0]), lower=1, trans=0)
+        LiW12, _ = dtrtrs(L, np.diagflat(W_12), lower=1, trans=0)
         K_Wi_i = np.dot(LiW12.T, LiW12) # R = W12BiW12, in R&W p 126, eq 5.25
 
+        #here's a better way to compute the required matrix. 
+        # you could do the model finding witha backsub, instead of a dot...
+        #L2 = L/W_12
+        #K_Wi_i_2 , _= dpotri(L2)
+        #symmetrify(K_Wi_i_2)
+
+
         return K_Wi_i, L, LiW12