diff --git a/GPy/kern/parts/linear.py b/GPy/kern/parts/linear.py
index c1da3944..04bd78a4 100644
--- a/GPy/kern/parts/linear.py
+++ b/GPy/kern/parts/linear.py
@@ -5,6 +5,7 @@
 from kernpart import Kernpart
 import numpy as np
 from ...util.linalg import tdot
+from ...util.misc import fast_array_equal
 from scipy import weave
 
 class Linear(Kernpart):
@@ -266,7 +267,7 @@ class Linear(Kernpart):
     #---------------------------------------#
 
     def _K_computations(self, X, X2):
-        if not (np.array_equal(X, self._Xcache) and np.array_equal(X2, self._X2cache)):
+        if not (fast_array_equal(X, self._Xcache) and fast_array_equal(X2, self._X2cache)):
             self._Xcache = X.copy()
             if X2 is None:
                 self._dot_product = tdot(X)
@@ -277,8 +278,8 @@ class Linear(Kernpart):
 
     def _psi_computations(self, Z, mu, S):
         # here are the "statistics" for psi1 and psi2
-        Zv_changed = not (np.array_equal(Z, self._Z) and np.array_equal(self.variances, self._variances))
-        muS_changed = not (np.array_equal(mu, self._mu) and np.array_equal(S, self._S))
+        Zv_changed = not (fast_array_equal(Z, self._Z) and fast_array_equal(self.variances, self._variances))
+        muS_changed = not (fast_array_equal(mu, self._mu) and fast_array_equal(S, self._S))
         if Zv_changed:
             # Z has changed, compute Z specific stuff
             # self.ZZ = Z[:,None,:]*Z[None,:,:] # num_inducing,num_inducing,input_dim
diff --git a/GPy/kern/parts/rbf.py b/GPy/kern/parts/rbf.py
index fce1722b..d317e4b7 100644
--- a/GPy/kern/parts/rbf.py
+++ b/GPy/kern/parts/rbf.py
@@ -7,6 +7,7 @@ import numpy as np
 import hashlib
 from scipy import weave
 from ...util.linalg import tdot
+from ...util.misc import fast_array_equal
 
 class RBF(Kernpart):
     """
@@ -222,7 +223,7 @@ class RBF(Kernpart):
     #---------------------------------------#
 
     def _K_computations(self, X, X2):
-        if not (np.array_equal(X, self._X) and np.array_equal(X2, self._X2) and np.array_equal(self._params , self._get_params())):
+        if not (fast_array_equal(X, self._X) and fast_array_equal(X2, self._X2) and fast_array_equal(self._params , self._get_params())):
             self._X = X.copy()
             self._params == self._get_params().copy()
             if X2 is None:
@@ -239,14 +240,14 @@ class RBF(Kernpart):
 
     def _psi_computations(self, Z, mu, S):
         # here are the "statistics" for psi1 and psi2
-        if not np.array_equal(Z, self._Z):
+        if not fast_array_equal(Z, self._Z):
             #Z has changed, compute Z specific stuff
             self._psi2_Zhat = 0.5*(Z[:,None,:] +Z[None,:,:]) # M,M,Q
             self._psi2_Zdist = 0.5*(Z[:,None,:]-Z[None,:,:]) # M,M,Q
             self._psi2_Zdist_sq = np.square(self._psi2_Zdist/self.lengthscale) # M,M,Q
             self._Z = Z
 
-        if not (np.array_equal(Z, self._Z) and np.array_equal(mu, self._mu) and np.array_equal(S, self._S)):
+        if not (fast_array_equal(Z, self._Z) and fast_array_equal(mu, self._mu) and fast_array_equal(S, self._S)):
             #something's changed. recompute EVERYTHING
 
             #psi1
diff --git a/GPy/util/misc.py b/GPy/util/misc.py
index e0f70703..5b3a70da 100644
--- a/GPy/util/misc.py
+++ b/GPy/util/misc.py
@@ -1,8 +1,8 @@
 # Copyright (c) 2012, GPy authors (see AUTHORS.txt).
 # Licensed under the BSD 3-clause license (see LICENSE.txt)
 
-
 import numpy as np
+from scipy import weave
 
 def opt_wrapper(m, **kwargs):
     """
@@ -58,6 +58,47 @@ def kmm_init(X, m = 10):
     inducing = np.array(inducing)
     return X[inducing]
 
+def fast_array_equal(A, B):
+    code="""
+    int i, j;
+    return_val = 1;
+
+    #pragma omp parallel for private(i, j)
+    for(i=0;i<N;i++){
+       for(j=0;j<D;j++){
+          if(A(i, j) != B(i, j)){
+              return_val = 0;
+              break;
+          }
+       }
+    }
+    """
+
+    support_code = """
+    #include <omp.h>
+    #include <math.h>
+    """
+
+    weave_options = {'headers'           : ['<omp.h>'],
+                     'extra_compile_args': ['-fopenmp -O3'],
+                     'extra_link_args'   : ['-lgomp']}
+
+
+    value = False
+
+    if A is not None and B is not None and A.shape == B.shape:
+        if len(A.shape) == 2:
+            N, D = A.shape
+            value = weave.inline(code, support_code=support_code, libraries=['gomp'],
+                                 arg_names=['A', 'B', 'N', 'D'],
+                                 type_converters=weave.converters.blitz,**weave_options)
+        else:
+            value = np.array_equal(A,B)
+
+
+    return value
+
+
 if __name__ == '__main__':
     import pylab as plt
     X = np.linspace(1,10, 100)[:, None]