diff --git a/GPy/kern/_src/basis_funcs.py b/GPy/kern/_src/basis_funcs.py
index e18be3fb..a6c1f36c 100644
--- a/GPy/kern/_src/basis_funcs.py
+++ b/GPy/kern/_src/basis_funcs.py
@@ -121,7 +121,7 @@ class LinearSlopeBasisFuncKernel(BasisFuncKernel):
     
 class ChangePointBasisFuncKernel(BasisFuncKernel):
     def __init__(self, input_dim, changepoint, variance=1., active_dims=None, ARD=False, name='changepoint'):
-        self.changepoint = changepoint
+        self.changepoint = np.array(changepoint)
         super(ChangePointBasisFuncKernel, self).__init__(input_dim, variance, active_dims, ARD, name)
     
     @Cache_this(limit=3, ignore_args=())
diff --git a/GPy/testing/linalg_test.py b/GPy/testing/linalg_test.py
index 81cb0368..ec3aca5a 100644
--- a/GPy/testing/linalg_test.py
+++ b/GPy/testing/linalg_test.py
@@ -50,3 +50,11 @@ class LinalgTests(np.testing.TestCase):
         pure = np.einsum('ijk,jlk->il', A, B)
         quick = GPy.util.linalg.ijk_jlk_to_il(A,B)
         np.testing.assert_allclose(pure, quick)
+
+    def test_einsum_ijk_ljk_to_ilk(self):
+        A = np.random.randn(150, 20, 5)
+        B = np.random.randn(150, 20, 5)
+        #B = A.copy()
+        pure = np.einsum('ijk,ljk->ilk', A, B)
+        quick = GPy.util.linalg.ijk_ljk_to_ilk(A,B)
+        np.testing.assert_allclose(pure, quick)
diff --git a/GPy/util/linalg.py b/GPy/util/linalg.py
index 285701e7..634a1e0d 100644
--- a/GPy/util/linalg.py
+++ b/GPy/util/linalg.py
@@ -463,5 +463,16 @@ def ijk_jlk_to_il(A, B):
     Faster version of einsum einsum('ijk,jlk->il', A,B)
     """
     res = np.zeros((A.shape[0], B.shape[1]))
-    [np.add(np.dot(A[:,:,k], B[:,:,k]), res, res) for k in range(B.shape[-1])]
+    [np.add(np.dot(A[:,:,k], B[:,:,k]), res, out=res) for k in range(B.shape[-1])]
+    return res
+
+def ijk_ljk_to_ilk(A, B):
+    """
+    Faster version of einsum np.einsum('ijk,ljk->ilk', A, B)
+
+    I.e A.dot(B.T) for every dimension
+    """
+    res = np.zeros((A.shape[-1], A.shape[0], B.shape[0]))
+    [np.dot(A[:,:,i], B[:,:,i].T, out=res[i,:,:]) for i in range(A.shape[-1])]
+    res = res.swapaxes(0, 2).swapaxes(0,1)
     return res