docs: def _expectation(p, kern, feat, none1, none2, nghp=None):
      """
      Compute the expectation:
      <K_{X, Z}>_p(X)
          - K_{.,.} :: RBF kernel
  
      :return: NxM
      """
      with params_as_tensors_for(kern, feat):
          # use only active dimensions
          Xcov = kern._slice_cov(p.cov)
          Z, Xmu = kern._slice(feat.Z, p.mu)
          D = tf.shape(Xmu)[1]
          if kern.ARD:
              lengthscales = kern.lengthscales
          else:
              lengthscales = tf.zeros((D,), dtype=settings.float_type) + kern.lengthscales
  
          chol_L_plus_Xcov = tf.cholesky(tf.matrix_diag(lengthscales ** 2) + Xcov)  # NxDxD
  
          all_diffs = tf.transpose(Z) - tf.expand_dims(Xmu, 2)  # NxDxM
          exponent_mahalanobis = tf.matrix_triangular_solve(chol_L_plus_Xcov, all_diffs, lower=True)  # NxDxM
          exponent_mahalanobis = tf.reduce_sum(tf.square(exponent_mahalanobis), 1)  # NxM
          exponent_mahalanobis = tf.exp(-0.5 * exponent_mahalanobis)  # NxM
  
          sqrt_det_L = tf.reduce_prod(lengthscales)
          sqrt_det_L_plus_Xcov = tf.exp(tf.reduce_sum(tf.log(tf.matrix_diag_part(chol_L_plus_Xcov)), axis=1))
          determinants = sqrt_det_L / sqrt_det_L_plus_Xcov  # N
  
          return kern.variance * (determinants[:, None] * exponent_mahalanobis)