diff --git a/GPy/core/gp.py b/GPy/core/gp.py
index e0f5755c..ba7c38f2 100644
--- a/GPy/core/gp.py
+++ b/GPy/core/gp.py
@@ -148,6 +148,28 @@ class GP(Model):
         m, v = self._raw_predict(X,  full_cov=False)
         return self.likelihood.predictive_quantiles(m, v, quantiles, Y_metadata)
 
+    def predictive_gradients(self, Xnew):
+        """
+        Compute the derivatives of the latent function with respect to X*
+
+        Given a set of points at which to predict X* (size [N*,Q]), compute the
+        derivatives of the mean and variance. Resulting arrays are sized:
+         dmu_dX* -- [N*, Q ,D], where D is the number of output in this GP (usually one).
+         dv_dX*  -- [N*, Q],    (since all outputs have the same variance)
+
+        """
+        dmu_dX = np.empty((Xnew.shape[0],Xnew.shape[1],self.output_dim))
+        for i in range(self.output_dim):
+            dmu_dX[:,:,i] = self.kern.gradients_X(self.posterior.woodbury_vector[:,i:i+1].T, Xnew, self.X)
+
+        # gradients wrt the diagonal part k_{xx}
+        dv_dX = self.kern.gradients_X(np.eye(Xnew.shape[0]), Xnew)
+        #grads wrt 'Schur' part K_{xf}K_{ff}^{-1}K_{fx}
+        alpha = -2.*np.dot(self.kern.K(Xnew, self.X),self.posterior.woodbury_inv)
+        dv_dX += self.kern.gradients_X(alpha, Xnew, self.X)
+        return dmu_dX, dv_dX
+
+
     def posterior_samples_f(self,X,size=10, full_cov=True):
         """
         Samples the posterior GP at the points X.
diff --git a/GPy/kern/_src/stationary.py b/GPy/kern/_src/stationary.py
index 74ac525e..c0553238 100644
--- a/GPy/kern/_src/stationary.py
+++ b/GPy/kern/_src/stationary.py
@@ -160,23 +160,20 @@ class Stationary(Kern):
         """
         invdist = self._inv_dist(X, X2)
         dL_dr = self.dK_dr_via_X(X, X2) * dL_dK
+        tmp = invdist*dL_dr
+        if X2 is None:
+            tmp = tmp + tmp.T
+            X2 = X
+
         #The high-memory numpy way:
         #d =  X[:, None, :] - X2[None, :, :]
-        #ret = np.sum((invdist*dL_dr)[:,:,None]*d,1)/self.lengthscale**2
-        #if X2 is None:
-            #ret *= 2.
+        #ret = np.sum(tmp[:,:,None]*d,1)/self.lengthscale**2
 
         #the lower memory way with a loop
-        tmp = invdist*dL_dr
         ret = np.empty(X.shape, dtype=np.float64)
-        if X2 is None:
-            [np.einsum('ij,ij->i', tmp, X[:,q][:,None]-X[:,q][None,:], out=ret[:,q]) for q in xrange(self.input_dim)]
-            ret2 = np.empty(X.shape, dtype=np.float64)
-            [np.einsum('ij,ji->j', tmp, X[:,q][:,None]-X[:,q][None,:], out=ret2[:,q]) for q in xrange(self.input_dim)]
-            ret += ret2
-        else:
-            [np.einsum('ij,ij->i', tmp, X[:,q][:,None]-X2[:,q][None,:], out=ret[:,q]) for q in xrange(self.input_dim)]
+        [np.sum(tmp*(X[:,q][:,None]-X2[:,q][None,:]), axis=1, out=ret[:,q]) for q in xrange(self.input_dim)]
         ret /= self.lengthscale**2
+
         return ret
 
     def gradients_X_diag(self, dL_dKdiag, X):
diff --git a/GPy/models/gplvm.py b/GPy/models/gplvm.py
index 542dcd31..8f5432ba 100644
--- a/GPy/models/gplvm.py
+++ b/GPy/models/gplvm.py
@@ -45,10 +45,10 @@ class GPLVM(GP):
         self.X.gradient = self.kern.gradients_X(self.grad_dict['dL_dK'], self.X, None)
 
     def jacobian(self,X):
-        target = np.zeros((X.shape[0],X.shape[1],self.output_dim))
+        J = np.zeros((X.shape[0],X.shape[1],self.output_dim))
         for i in range(self.output_dim):
-            target[:,:,i]=self.kern.gradients_X(np.dot(self.Ki,self.likelihood.Y[:,i])[None, :],X,self.X)
-        return target
+            J[:,:,i] = self.kern.gradients_X(self.posterior.woodbury_vector[:,i:i+1], X, self.X)
+        return J
 
     def magnification(self,X):
         target=np.zeros(X.shape[0])
diff --git a/GPy/util/multioutput.py b/GPy/util/multioutput.py
index b769d6a5..cc9af29e 100644
--- a/GPy/util/multioutput.py
+++ b/GPy/util/multioutput.py
@@ -40,7 +40,7 @@ def build_likelihood(Y_list,noise_index,likelihoods_list=None):
     return likelihood
 
 
-def ICM(input_dim, num_outputs, kernel, W_rank=1,W=None,kappa=None,name='X'):
+def ICM(input_dim, num_outputs, kernel, W_rank=1,W=None,kappa=None,name='ICM'):
     """
     Builds a kernel for an Intrinsic Coregionalization Model
 
@@ -56,12 +56,10 @@ def ICM(input_dim, num_outputs, kernel, W_rank=1,W=None,kappa=None,name='X'):
         warnings.warn("kernel's input dimension overwritten to fit input_dim parameter.")
 
     K = kernel.prod(GPy.kern.Coregionalize(1, num_outputs, active_dims=[input_dim], rank=W_rank,W=W,kappa=kappa,name='B'),name=name)
-    K['.*variance'] = 1.
-    K['.*variance'].fix()
     return K
 
 
-def LCM(input_dim, num_outputs, kernels_list, W_rank=1,name='X'):
+def LCM(input_dim, num_outputs, kernels_list, W_rank=1,name='ICM'):
     """
     Builds a kernel for an Linear Coregionalization Model
 
@@ -77,6 +75,7 @@ def LCM(input_dim, num_outputs, kernels_list, W_rank=1,name='X'):
     j = 1
     for kernel in kernels_list[1:]:
         K += ICM(input_dim,num_outputs,kernel,W_rank,name='%s%s' %(name,j))
+        j += 1
     return K