docs: def _expectation(p, kern1, feat1, kern2, feat2, nghp=None):
      r"""
      Compute the expectation:
      expectation[n] = < prodK_{Z, x_n} prodK_{x_n, Z} >_p(x_n)
                     = < (\HadamardProd_i Ki_{Z[:, active_dims_i], x[n, active_dims_i]})  <-- Mx1
                 1xM -->  (\HadamardProd_j Kj_{x[n, active_dims_j], Z[:, active_dims_j]}) >_p(x_n)  (MxM)
  
          - \HadamardProd_i Ki_{.,.}, \HadamardProd_j Kj_{.,.} :: Product kernels
          - p                        :: DiagonalGaussian distribution (p.cov NxD)
  
      :return: NxMxM
      """
      if feat1 != feat2:
          raise NotImplementedError("Different features are not supported.")
      if kern1 != kern2:
          raise NotImplementedError("Calculating the expectation over two "
                                    "different Product kernels is not supported.")
  
      kern = kern1
      feat = feat1
  
      if not kern.on_separate_dimensions:
          raise NotImplementedError(
              "Product currently needs to be defined on separate dimensions.")  # pragma: no cover
  
      return functools.reduce(tf.multiply, [
          expectation(p, (k, feat), (k, feat), nghp=nghp) for k in kern.kernels])