diff --git a/GPy/kern/_src/kern.py b/GPy/kern/_src/kern.py
index e2647901..db440cab 100644
--- a/GPy/kern/_src/kern.py
+++ b/GPy/kern/_src/kern.py
@@ -6,6 +6,7 @@ import numpy as np
 from ...core.parameterization.parameterized import Parameterized
 from kernel_slice_operations import KernCallsViaSlicerMeta
 from ...util.caching import Cache_this
+from GPy.core.parameterization.observable_array import ObsAr
 
 
 
@@ -54,6 +55,20 @@ class Kern(Parameterized):
 
         self._sliced_X = 0
         self.useGPU = self._support_GPU and useGPU
+        self._return_psi2_n_flag = ObsAr(np.zeros(1)).astype(bool)
+
+    @property
+    def return_psi2_n(self):
+        """
+        Flag whether to pass back psi2 as NxMxM or MxM, by summing out N.
+        """
+        return self._return_psi2_n_flag[0]
+    @return_psi2_n.setter
+    def return_psi2_n(self, val):
+        def visit(self):
+            if isinstance(self, Kern):
+                self._return_psi2_n_flag[0]=val
+        self.traverse(visit)
 
     @Cache_this(limit=20)
     def _slice_X(self, X):
@@ -162,7 +177,7 @@ class Kern(Parameterized):
     def __mul__(self, other):
         """ Here we overload the '*' operator. See self.prod for more information"""
         return self.prod(other)
-    
+
     def __imul__(self, other):
         """ Here we overload the '*' operator. See self.prod for more information"""
         return self.prod(other)
diff --git a/GPy/kern/_src/psi_comp/__init__.py b/GPy/kern/_src/psi_comp/__init__.py
index eced8392..a277ff02 100644
--- a/GPy/kern/_src/psi_comp/__init__.py
+++ b/GPy/kern/_src/psi_comp/__init__.py
@@ -10,7 +10,6 @@ import sslinear_psi_comp
 import linear_psi_comp
 
 class PSICOMP_RBF(Pickleable):
-
     @Cache_this(limit=2, ignore_args=(0,))
     def psicomputations(self, variance, lengthscale, Z, variational_posterior):
         if isinstance(variational_posterior, variational.NormalPosterior):
@@ -19,7 +18,7 @@ class PSICOMP_RBF(Pickleable):
             return ssrbf_psi_comp.psicomputations(variance, lengthscale, Z, variational_posterior)
         else:
             raise ValueError, "unknown distriubtion received for psi-statistics"
-                
+
     @Cache_this(limit=2, ignore_args=(0,1,2,3))
     def psiDerivativecomputations(self, dL_dpsi0, dL_dpsi1, dL_dpsi2, variance, lengthscale, Z, variational_posterior):
         if isinstance(variational_posterior, variational.NormalPosterior):
@@ -28,10 +27,10 @@ class PSICOMP_RBF(Pickleable):
             return ssrbf_psi_comp.psiDerivativecomputations(dL_dpsi0, dL_dpsi1, dL_dpsi2, variance, lengthscale, Z, variational_posterior)
         else:
             raise ValueError, "unknown distriubtion received for psi-statistics"
-        
+
     def _setup_observers(self):
         pass
-        
+
 class PSICOMP_Linear(Pickleable):
 
     @Cache_this(limit=2, ignore_args=(0,))
@@ -42,7 +41,7 @@ class PSICOMP_Linear(Pickleable):
             return sslinear_psi_comp.psicomputations(variance, Z, variational_posterior)
         else:
             raise ValueError, "unknown distriubtion received for psi-statistics"
-                
+
     @Cache_this(limit=2, ignore_args=(0,1,2,3))
     def psiDerivativecomputations(self, dL_dpsi0, dL_dpsi1, dL_dpsi2, variance, Z, variational_posterior):
         if isinstance(variational_posterior, variational.NormalPosterior):
@@ -51,6 +50,6 @@ class PSICOMP_Linear(Pickleable):
             return sslinear_psi_comp.psiDerivativecomputations(dL_dpsi0, dL_dpsi1, dL_dpsi2, variance, Z, variational_posterior)
         else:
             raise ValueError, "unknown distriubtion received for psi-statistics"
-        
+
     def _setup_observers(self):
         pass
\ No newline at end of file
diff --git a/GPy/kern/_src/psi_comp/rbf_psi_comp.py b/GPy/kern/_src/psi_comp/rbf_psi_comp.py
index d2d3312d..e7a67d16 100644
--- a/GPy/kern/_src/psi_comp/rbf_psi_comp.py
+++ b/GPy/kern/_src/psi_comp/rbf_psi_comp.py
@@ -139,7 +139,7 @@ def _psi2compDer(dL_dpsi2, variance, lengthscale, Z, mu, S):
     denom2 = np.square(denom)
 
     _psi2 = _psi2computations(variance, lengthscale, Z, mu, S) # NxMxM
-    Lpsi2 = dL_dpsi2[None,:,:]*_psi2
+    Lpsi2 = dL_dpsi2*_psi2 # dL_dpsi2 is MxM, using broadcast to multiply N out
     Lpsi2sum = np.einsum('nmo->n',Lpsi2) #N
     Lpsi2Z = np.einsum('nmo,oq->nq',Lpsi2,Z) #NxQ
     Lpsi2Z2 = np.einsum('nmo,oq,oq->nq',Lpsi2,Z,Z) #NxQ