diff --git a/GPy/kern/_src/prod.py b/GPy/kern/_src/prod.py
index 5e4c0d29..a3b49973 100644
--- a/GPy/kern/_src/prod.py
+++ b/GPy/kern/_src/prod.py
@@ -14,8 +14,9 @@ def numpy_invalid_op_as_exception(func):
     """
     def func_wrapper(*args, **kwargs):
         np.seterr(invalid='raise')
-        func(*args, **kwargs)
+        result = func(*args, **kwargs)
         np.seterr(invalid='warn')
+        return result
     return func_wrapper
 
 
@@ -66,17 +67,23 @@ class Prod(CombinationKernel):
                 to_update = list(set(self.parts) - set(combination))[0]
                 to_update.update_gradients_full(dL_dK * prod, X, X2)
 
-
     def update_gradients_diag(self, dL_dKdiag, X):
         k = self.Kdiag(X)*dL_dKdiag
         for p in self.parts:
             p.update_gradients_diag(k/p.Kdiag(X),X)
 
+    @numpy_invalid_op_as_exception            
     def gradients_X(self, dL_dK, X, X2=None):
         target = np.zeros(X.shape)
         k = self.K(X,X2)*dL_dK
-        for p in self.parts:
-            target += p.gradients_X(k/p.K(X,X2),X,X2)
+        try:
+            for p in self.parts:
+                target += p.gradients_X(k/p.K(X,X2),X,X2)
+        except FloatingPointError:
+            for combination in itertools.combinations(self.parts, len(self.parts) - 1):
+                prod = reduce(np.multiply, [p.K(X, X2) for p in combination])
+                to_update = list(set(self.parts) - set(combination))[0]
+                target += to_update.gradients_X(dL_dK * prod, X, X2)
         return target
 
     def gradients_X_diag(self, dL_dKdiag, X):
diff --git a/GPy/testing/kernel_tests.py b/GPy/testing/kernel_tests.py
index f9d90607..415cc7eb 100644
--- a/GPy/testing/kernel_tests.py
+++ b/GPy/testing/kernel_tests.py
@@ -400,25 +400,27 @@ class Coregionalize_weave_test(unittest.TestCase):
     #reset the weave state for any other tests
     GPy.util.config.config.set('weave', 'working', 'False')
 
+
 class KernelTestsProductWithZeroValues(unittest.TestCase):
 
-    def test_zero_valued_kernel(self):
-        X = np.array([[0,1],[1,0]])
-        Y = np.array([[1],[10]])
-        lin = GPy.kern.Linear(2)
-        bias = GPy.kern.Bias(2)
-        k = lin * bias
-        k.update_gradients_full(1, X)
-        self.assertFalse(np.isnan(k['linear.variances'].gradient),
+    def setUp(self):
+        self.X = np.array([[0,1],[1,0]])
+        self.k = GPy.kern.Linear(2) * GPy.kern.Bias(2)
+
+    def test_zero_valued_kernel_full(self):
+        self.k.update_gradients_full(1, self.X)
+        self.assertFalse(np.isnan(self.k['linear.variances'].gradient),
                          "Gradient resulted in NaN")
 
+    def test_zero_valued_kernel_gradients_X(self):
+        target = self.k.gradients_X(1, self.X)
+        self.assertFalse(np.any(np.isnan(target)),
+                         "Gradient resulted in NaN")
 
 
 if __name__ == "__main__":
     print "Running unit tests, please be (very) patient..."
     unittest.main()
-    #suite = unittest.TestLoader().loadTestsFromTestCase(KernelTestsProductWithZeroValues)
-    #unittest.TextTestRunner().run(suite)
 
 #     np.random.seed(0)
 #     N0 = 3