diff --git a/GPy/testing/linalg_test.py b/GPy/testing/linalg_test.py
index ad524156..ec3aca5a 100644
--- a/GPy/testing/linalg_test.py
+++ b/GPy/testing/linalg_test.py
@@ -53,7 +53,8 @@ class LinalgTests(np.testing.TestCase):
 
     def test_einsum_ijk_ljk_to_ilk(self):
         A = np.random.randn(150, 20, 5)
-        B = np.random.randn(20, 30, 5)
-        pure = np.einsum('ijk,jlk->il', A, B)
-        quick = GPy.util.linalg.ijk_jlk_to_il(A,B)
+        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 bc047654..634a1e0d 100644
--- a/GPy/util/linalg.py
+++ b/GPy/util/linalg.py
@@ -463,7 +463,7 @@ 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):
@@ -472,7 +472,7 @@ def ijk_ljk_to_ilk(A, B):
 
     I.e A.dot(B.T) for every dimension
     """
-    res = np.empty((A.shape[0], B.shape[0], A.shape[-1]))
-    [np.dot(A[:,:,i], B[:,:,i].T, res[i,:,:]) for i in range(A.shape[0])]
-    res = res.swapaxes(0,2)
+    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