diff --git a/GPy/kern/src/linear.py b/GPy/kern/src/linear.py
index 25f6299d..0c897c74 100644
--- a/GPy/kern/src/linear.py
+++ b/GPy/kern/src/linear.py
@@ -73,14 +73,15 @@ class Linear(Kern):
         return np.sum(self.variances * np.square(X), -1)
 
     def update_gradients_full(self, dL_dK, X, X2=None):
+        if X2 is None: dL_dK = (dL_dK+dL_dK.T)/2
         if self.ARD:
             if X2 is None:
                 #self.variances.gradient = np.array([np.sum(dL_dK * tdot(X[:, i:i + 1])) for i in range(self.input_dim)])
-                self.variances.gradient = np.einsum('ij,iq,jq->q', dL_dK, X, X)
+                self.variances.gradient = (dL_dK.dot(X)*X).sum(0) #np.einsum('ij,iq,jq->q', dL_dK, X, X)
             else:
                 #product = X[:, None, :] * X2[None, :, :]
                 #self.variances.gradient = (dL_dK[:, :, None] * product).sum(0).sum(0)
-                self.variances.gradient = np.einsum('ij,iq,jq->q', dL_dK, X, X2)
+                self.variances.gradient = (dL_dK.dot(X2)*X).sum(0)  #np.einsum('ij,iq,jq->q', dL_dK, X, X2)
         else:
             self.variances.gradient = np.sum(self._dot_product(X, X2) * dL_dK)
 
@@ -93,13 +94,15 @@ class Linear(Kern):
 
 
     def gradients_X(self, dL_dK, X, X2=None):
+        if X2 is None: dL_dK = (dL_dK+dL_dK.T)/2
         if X2 is None:
-            return np.einsum('jq,q,ij->iq', X, 2*self.variances, dL_dK)
+            return dL_dK.dot(X)*(2*self.variances) #np.einsum('jq,q,ij->iq', X, 2*self.variances, dL_dK)
         else:
             #return (((X2[None,:, :] * self.variances)) * dL_dK[:, :, None]).sum(1)
-            return np.einsum('jq,q,ij->iq', X2, self.variances, dL_dK)
+            return dL_dK.dot(X2)*self.variances #np.einsum('jq,q,ij->iq', X2, self.variances, dL_dK)
 
     def gradients_XX(self, dL_dK, X, X2=None):
+        if X2 is None: dL_dK = (dL_dK+dL_dK.T)/2
         if X2 is None:
             return 2*np.ones(X.shape)*self.variances
         else:
@@ -162,6 +165,7 @@ class LinearFull(Kern):
         return np.einsum('ij,jk,lk->il', X, P, X if X2 is None else X2)
 
     def update_gradients_full(self, dL_dK, X, X2=None):
+        if X2 is None: dL_dK = (dL_dK+dL_dK.T)/2
         self.kappa.gradient = np.einsum('ij,ik,kj->j', X, dL_dK, X if X2 is None else X2)
         self.W.gradient = np.einsum('ij,kl,ik,lm->jm', X, X if X2 is None else X2, dL_dK, self.W)
         self.W.gradient += np.einsum('ij,kl,ik,jm->lm', X, X if X2 is None else X2, dL_dK, self.W)
@@ -175,6 +179,7 @@ class LinearFull(Kern):
         self.W.gradient = 2.*np.einsum('ij,ik,jl,i->kl', X, X, self.W, dL_dKdiag)
 
     def gradients_X(self, dL_dK, X, X2=None):
+        if X2 is None: dL_dK = (dL_dK+dL_dK.T)/2
         P = np.dot(self.W, self.W.T) + np.diag(self.kappa)
         if X2 is None:
             return 2.*np.einsum('ij,jk,kl->il', dL_dK, X, P)
diff --git a/GPy/testing/kernel_tests.py b/GPy/testing/kernel_tests.py
index e2f4324a..bbfb565b 100644
--- a/GPy/testing/kernel_tests.py
+++ b/GPy/testing/kernel_tests.py
@@ -32,9 +32,9 @@ class Kern_check_model(GPy.core.Model):
             X = np.random.randn(20, kernel.input_dim)
         if dL_dK is None:
             if X2 is None:
-                dL_dK = np.ones((X.shape[0], X.shape[0]))
+                dL_dK = np.random.rand(X.shape[0], X.shape[0])
             else:
-                dL_dK = np.ones((X.shape[0], X2.shape[0]))
+                dL_dK = np.random.rand(X.shape[0], X2.shape[0])
 
         self.kernel = kernel
         self.X = X
@@ -311,7 +311,7 @@ class KernelGradientTestsContinuous(unittest.TestCase):
         self.assertTrue(check_kernel_gradient_functions(k, X=self.X, X2=self.X2, verbose=verbose))
 
     def test_RBF(self):
-        k = GPy.kern.RBF(self.D)
+        k = GPy.kern.RBF(self.D, ARD=True)
         k.randomize()
         self.assertTrue(check_kernel_gradient_functions(k, X=self.X, X2=self.X2, verbose=verbose))