docs: def _expectation(p, mean1, none1, mean2, none2, nghp=None):
      """
      Compute the expectation:
      expectation[n] = <m1(x_n)^T m2(x_n)>_p(x_n)
          - m1(.) :: Identity mean function
          - m2(.) :: Linear mean function
  
      :return: NxDxQ
      """
      with params_as_tensors_for(mean2):
          N = tf.shape(p.mu)[0]
          e_xxt = p.cov + (p.mu[:, :, None] * p.mu[:, None, :])  # NxDxD
          e_xxt_A = tf.matmul(e_xxt, tf.tile(mean2.A[None, ...], (N, 1, 1)))  # NxDxQ
          e_x_bt = p.mu[:, :, None] * mean2.b[None, None, :]  # NxDxQ
  
          return e_xxt_A + e_x_bt