diff --git a/GPy/models/sparse_GP_regression.py b/GPy/models/sparse_GP_regression.py
index 561aaca4..6313f4f7 100644
--- a/GPy/models/sparse_GP_regression.py
+++ b/GPy/models/sparse_GP_regression.py
@@ -269,7 +269,9 @@ class sgp_debugE(sparse_GP_regression):
         self.dL_dpsi2 = - 0.5 * self.beta * (self.G)
 
         # Compute dL_dKmm
-        self.dL_dKmm =  np.dot(np.dot(self.G,self.beta*self.psi2) - np.dot(self.LBL_inv, self.psi1VVpsi1), self.Kmmi) + 0.5*self.G # dE
+        tmp = mdot(self.beta*self.psi2, self.Kmmi, self.psi1VVpsi1)
+        self.dL_dKmm = -0.5*mdot(self.Kmmi,tmp + tmp.T + self.psi1VVpsi1,self.Kmmi)
+        #self.dL_dKmm =  np.dot(np.dot(self.G,self.beta*self.psi2) - np.dot(self.LBL_inv, self.psi1VVpsi1), self.Kmmi) + 0.5*self.G # dE
 
     def log_likelihood(self):
         A = -0.5*self.N*self.D*(np.log(2.*np.pi) - np.log(self.beta))
diff --git a/GPy/util/linalg.py b/GPy/util/linalg.py
index 092589a9..7414eb29 100644
--- a/GPy/util/linalg.py
+++ b/GPy/util/linalg.py
@@ -63,8 +63,8 @@ def jitchol(A,maxtries=5):
             raise linalg.LinAlgError, "not pd: negative diagonal elements"
         jitter= diagA.mean()*1e-6
         for i in range(1,maxtries+1):
+            print 'Warning: adding jitter of '+str(jitter)
             try:
-                print 'Warning: adding jitter of '+str(jitter)
                 return linalg.cholesky(A+np.eye(A.shape[0])*jitter, lower = True)
             except:
                 jitter *= 10
diff --git a/grid_parameters.py b/grid_parameters.py
index e15aab6e..011204ee 100644
--- a/grid_parameters.py
+++ b/grid_parameters.py
@@ -8,41 +8,50 @@ pb.close('all')
 
 N = 1000
 M = 10
-resolution=3
+resolution=5
 
 X = np.linspace(0,12,N)[:,None]
 Z = np.linspace(0,12,M)[:,None] # inducing points (fixed for now)
 Y = np.sin(X) + np.random.randn(*X.shape)/np.sqrt(50.)
-k = GPy.kern.rbf(1)
-
+#k = GPy.kern.rbf(1)
+k = GPy.kern.Matern32(1) + GPy.kern.white(1)
 
 models = [GPy.models.sparse_GP_regression(X,Y,Z=Z,kernel=k),
     GPy.models.sgp_debugB(X,Y,Z=Z,kernel=k),
-    GPy.models.sgp_debugC(X,Y,Z=Z,kernel=k),
-    GPy.models.sgp_debugE(X,Y,Z=Z,kernel=k)]
-#[m.constrain_fixed('iip') for m in models]
-#m.constrain_fixed('white',1e-6)
-#[m.constrain_fixed('precision',50) for m in models]
-#[m.ensure_default_constraints() for m in models]
+    GPy.models.sgp_debugC(X,Y,Z=Z,kernel=k)]#,
+    #GPy.models.sgp_debugE(X,Y,Z=Z,kernel=k)]
 
+[m.constrain_fixed('white',0.001) for m in models]
 
-xx,yy = np.mgrid[1.5:4:0+resolution*1j,-2:2:0+resolution*1j]
+#xx,yy = np.mgrid[1.5:4:0+resolution*1j,-2:2:0+resolution*1j]
+xx,yy = np.mgrid[3:16:0+resolution*1j,-2:2:0+resolution*1j]
 
 lls = []
 cgs = []
+grads = []
+count = 0
 for l,v in zip(xx.flatten(),yy.flatten()):
+    count += 1
+    print count, 'of', resolution**2
+    sys.stdout.flush()
+
     [m.set('lengthscale',l) for m in models]
-    [m.set('rbf_variance',10.**v) for m in models]
-    lls.append(models[0].log_likelihood())
+    [m.set('_variance',10.**v) for m in models]
+    lls.append([m.log_likelihood() for m in models])
+    grads.append([m.log_likelihood_gradients() for m in models])
     cgs.append([m.checkgrad(verbose=0,return_ratio=True) for m in models])
 
-lls = np.array(lls).reshape(resolution,resolution)
-cgs = np.array(zip(*cgs),dtype=np.float64).reshape(-1,resolution,resolution)
+lls = np.array(zip(*lls)).reshape(-1,resolution,resolution)
+cgs = np.array(zip(*cgs)).reshape(-1,resolution,resolution)
 
-for cg in cgs:
+for ll,cg in zip(lls,cgs):
     pb.figure()
-    pb.contourf(xx,yy,lls,50,cmap=pb.cm.gray)
+    pb.contourf(xx,yy,ll,100,cmap=pb.cm.gray)
     pb.colorbar()
+    try:
+        pb.contour(xx,yy,np.exp(ll),colors='k')
+    except:
+        pass
     pb.scatter(xx.flatten(),yy.flatten(),20,np.log(np.abs(cg.flatten())),cmap=pb.cm.jet,linewidth=0)
     pb.colorbar()