diff --git a/GPy/kern/_src/linear.py b/GPy/kern/_src/linear.py
index 7d9eeac2..609903aa 100644
--- a/GPy/kern/_src/linear.py
+++ b/GPy/kern/_src/linear.py
@@ -121,7 +121,7 @@ class Linear(Kern):
             gamma = variational_posterior.binary_prob
             mu = variational_posterior.mean
             return np.einsum('nq,q,mq,nq->nm',gamma,self.variances,Z,mu)
-#            return (self.variances*gamma*mu).sum(axis=1)       
+#            return (self.variances*gamma*mu).sum(axis=1)
         else:
             return self.K(variational_posterior.mean, Z) #the variance, it does nothing
 
@@ -177,7 +177,7 @@ class Linear(Kern):
 
             grad =  np.einsum('nm,nq,q,nq->mq',dL_dpsi1,gamma, self.variances,mu) +\
                  np.einsum('nmo,noq->mq',dL_dpsi2,_dpsi2_dZ)
-            
+
             return grad
         else:
             #psi1
@@ -191,15 +191,15 @@ class Linear(Kern):
             gamma = variational_posterior.binary_prob
             mu = variational_posterior.mean
             S = variational_posterior.variance
-            mu2S = np.square(mu)+S            
+            mu2S = np.square(mu)+S
             _, _dpsi2_dgamma, _dpsi2_dmu, _dpsi2_dS, _ = linear_psi_comp._psi2computations(self.variances, Z, mu, S, gamma)
-            
+
             grad_gamma = np.einsum('n,q,nq->nq',dL_dpsi0,self.variances,mu2S) + np.einsum('nm,q,mq,nq->nq',dL_dpsi1,self.variances,Z,mu) +\
                  np.einsum('nmo,nmoq->nq',dL_dpsi2,_dpsi2_dgamma)
             grad_mu = np.einsum('n,nq,q,nq->nq',dL_dpsi0,gamma,2.*self.variances,mu) + np.einsum('nm,nq,q,mq->nq',dL_dpsi1,gamma,self.variances,Z) +\
                  np.einsum('nmo,nmoq->nq',dL_dpsi2,_dpsi2_dmu)
             grad_S = np.einsum('n,nq,q->nq',dL_dpsi0,gamma,self.variances) + np.einsum('nmo,nmoq->nq',dL_dpsi2,_dpsi2_dS)
-            
+
             return grad_mu, grad_S, grad_gamma
         else:
             grad_mu, grad_S = np.zeros(variational_posterior.mean.shape), np.zeros(variational_posterior.mean.shape)
@@ -210,7 +210,7 @@ class Linear(Kern):
             grad_mu += (dL_dpsi1[:, :, None] * (Z * self.variances)).sum(1)
             # psi2
             self._weave_dpsi2_dmuS(dL_dpsi2, Z, variational_posterior, grad_mu, grad_S)
-    
+
             return grad_mu, grad_S
 
     #--------------------------------------------------#
diff --git a/GPy/kern/_src/prod.py b/GPy/kern/_src/prod.py
index e00f38c3..98b60366 100644
--- a/GPy/kern/_src/prod.py
+++ b/GPy/kern/_src/prod.py
@@ -58,8 +58,6 @@ class Prod(CombinationKernel):
     def gradients_X_diag(self, dL_dKdiag, X):
         target = np.zeros(X.shape)
         for k1,k2 in itertools.combinations(self.parts, 2):
-            target += k1.gradients_X(dL_dKdiag*k2.Kdiag(X), X)
-            target += k2.gradients_X(dL_dKdiag*k1.Kdiag(X), X)
+            target += k1.gradients_X_diag(dL_dKdiag*k2.Kdiag(X), X)
+            target += k2.gradients_X_diag(dL_dKdiag*k1.Kdiag(X), X)
         return target
-
-
diff --git a/GPy/testing/kernel_tests.py b/GPy/testing/kernel_tests.py
index 9ed218d8..36d55645 100644
--- a/GPy/testing/kernel_tests.py
+++ b/GPy/testing/kernel_tests.py
@@ -8,6 +8,7 @@ import sys
 
 verbose = 0
 
+np.random.seed(50)
 
 
 class Kern_check_model(GPy.core.Model):
@@ -243,6 +244,17 @@ class KernelGradientTestsContinuous(unittest.TestCase):
         k.randomize()
         self.assertTrue(check_kernel_gradient_functions(k, X=self.X, X2=self.X2, verbose=verbose))
 
+    def test_Prod2(self):
+        k = (GPy.kern.RBF(2, active_dims=[0,4]) * GPy.kern.Linear(self.D))
+        k.randomize()
+        self.assertTrue(check_kernel_gradient_functions(k, X=self.X, X2=self.X2, verbose=verbose))
+
+    def test_Prod3(self):
+        k = GPy.kern.Matern32(2, active_dims=[2,3]) * (GPy.kern.RBF(2, active_dims=[0,4]) + GPy.kern.Linear(self.D))
+        k = (GPy.kern.RBF(2, active_dims=[0,4]) * GPy.kern.Linear(self.D))
+        k.randomize()
+        self.assertTrue(check_kernel_gradient_functions(k, X=self.X, X2=self.X2, verbose=verbose))
+
     def test_Add(self):
         k = GPy.kern.Matern32(2, active_dims=[2,3]) + GPy.kern.RBF(2, active_dims=[0,4]) + GPy.kern.Linear(self.D)
         k += GPy.kern.Matern32(2, active_dims=[2,3]) + GPy.kern.RBF(2, active_dims=[0,4]) + GPy.kern.Linear(self.D)