diff --git a/GPy/likelihoods/EP.py b/GPy/likelihoods/EP.py
index 118b226a..8307b6b4 100644
--- a/GPy/likelihoods/EP.py
+++ b/GPy/likelihoods/EP.py
@@ -196,8 +196,9 @@ class EP(likelihood):
                 self.tau_tilde[i] = self.tau_tilde[i] + Delta_tau
                 self.v_tilde[i] = self.v_tilde[i] + Delta_v
                 #Posterior distribution parameters update
-                LLT = LLT + np.outer(Kmn[:,i],Kmn[:,i])*Delta_tau
-                L = jitchol(LLT)
+                #LLT = LLT + np.outer(Kmn[:,i],Kmn[:,i])*Delta_tau
+                #L = jitchol(LLT)
+                cholupdate(L,Kmn[:,i]*np.sqrt(Delta_tau))
                 V,info = linalg.lapack.flapack.dtrtrs(L,Kmn,lower=1)
                 Sigma_diag = np.sum(V*V,-2)
                 si = np.sum(V.T*V[:,i],-1)
diff --git a/GPy/models/sparse_GP.py b/GPy/models/sparse_GP.py
index 14c789b8..c3b9f793 100644
--- a/GPy/models/sparse_GP.py
+++ b/GPy/models/sparse_GP.py
@@ -9,6 +9,11 @@ from .. import kern
 from GP import GP
 from scipy import linalg
 
+def backsub_both_sides(L,X):
+    """ Return L^-T * X * L^-1, assumuing X is symmetrical and L is lower cholesky"""
+    tmp,_ = linalg.lapack.flapack.dtrtrs(L,np.asfortranarray(X),lower=1,trans=1)
+    return linalg.lapack.flapack.dtrtrs(L,np.asfortranarray(tmp.T),lower=1,trans=1)[0].T
+
 class sparse_GP(GP):
     """
     Variational sparse GP model
diff --git a/GPy/util/linalg.py b/GPy/util/linalg.py
index 42a98fea..a62fccb3 100644
--- a/GPy/util/linalg.py
+++ b/GPy/util/linalg.py
@@ -276,3 +276,30 @@ def symmetrify_murray(A):
     nn = A.shape[0]
     A[[range(nn),range(nn)]] /= 2.0
 
+def cholupdate(L,x):
+    """
+    update the LOWER cholesky factor of a pd matrix IN PLACE
+
+    if L is the lower chol. of K, then this function computes L_
+    where L_ is the lower chol of K + x*x^T
+    """
+    support_code = """
+    #include <math.h>
+    """
+    code="""
+    double r,c,s;
+    int j,i;
+    for(j=0; j<N; j++){
+      r = sqrt(L(j,j)*L(j,j) + x(j)*x(j));
+      c = r / L(j,j);
+      s = x(j) / L(j,j);
+      L(j,j) = r;
+      for (i=j+1; i<N; i++){
+        L(i,j) = (L(i,j) + s*x(i))/c;
+        x(i) = c*x(i) - s*L(i,j);
+      }
+    }
+    """
+    x = x.copy()
+    N = x.size
+    weave.inline(code, support_code=support_code, arg_names=['N','L','x'], type_converters=weave.converters.blitz)