docs: def _expectation(p, lin_kern, feat1, rbf_kern, feat2, nghp=None):
      """
      Compute the expectation:
      expectation[n] = <Ka_{Z1, x_n} Kb_{x_n, Z2}>_p(x_n)
          - K_lin_{.,.} :: Linear kernel
          - K_rbf_{.,.} :: RBF kernel
      Different Z1 and Z2 are handled if p is diagonal and K_lin and K_rbf have disjoint
      active_dims, in which case the joint expectations simplify into a product of expectations
  
      :return: NxM1xM2
      """
      return tf.matrix_transpose(expectation(p, (rbf_kern, feat2), (lin_kern, feat1)))