docs: def _expectation(p, rbf_kern, feat1, lin_kern, feat2, nghp=None):
      """
      Compute the expectation:
      expectation[n] = <Ka_{Z1, x_n} Kb_{x_n, Z2}>_p(x_n)
          - K_lin_{.,.} :: RBF kernel
          - K_rbf_{.,.} :: Linear 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
      """
      if rbf_kern.on_separate_dims(lin_kern) and isinstance(p, DiagonalGaussian):  # no joint expectations required
          eKxz1 = expectation(p, (rbf_kern, feat1))
          eKxz2 = expectation(p, (lin_kern, feat2))
          return eKxz1[:, :, None] * eKxz2[:, None, :]
  
      if feat1 != feat2:
          raise NotImplementedError("Features have to be the same for both kernels.")
  
      if rbf_kern.active_dims != lin_kern.active_dims:
          raise NotImplementedError("active_dims have to be the same for both kernels.")
  
      with params_as_tensors_for(rbf_kern, lin_kern, feat1, feat2):
          # use only active dimensions
          Xcov = rbf_kern._slice_cov(tf.matrix_diag(p.cov) if isinstance(p, DiagonalGaussian) else p.cov)
          Z, Xmu = rbf_kern._slice(feat1.Z, p.mu)
  
          N = tf.shape(Xmu)[0]
          D = tf.shape(Xmu)[1]
  
          lin_kern_variances = lin_kern.variance if lin_kern.ARD \
              else tf.zeros((D,), dtype=settings.float_type) + lin_kern.variance
  
          rbf_kern_lengthscales = rbf_kern.lengthscales if rbf_kern.ARD \
              else tf.zeros((D,), dtype=settings.float_type) + rbf_kern.lengthscales  ## Begin RBF eKxz code:
  
          chol_L_plus_Xcov = tf.cholesky(tf.matrix_diag(rbf_kern_lengthscales ** 2) + Xcov)  # NxDxD
  
          Z_transpose = tf.transpose(Z)
          all_diffs = Z_transpose - 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(rbf_kern_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
          eKxz_rbf = rbf_kern.variance * (determinants[:, None] * exponent_mahalanobis)  ## NxM <- End RBF eKxz code
  
          tiled_Z = tf.tile(tf.expand_dims(Z_transpose, 0), (N, 1, 1))  # NxDxM
          z_L_inv_Xcov = tf.matmul(tiled_Z, Xcov / rbf_kern_lengthscales[:, None] ** 2., transpose_a=True)  # NxMxD
  
          cross_eKzxKxz = tf.cholesky_solve(
              chol_L_plus_Xcov, (lin_kern_variances * rbf_kern_lengthscales ** 2.)[..., None] * tiled_Z)  # NxDxM
  
          cross_eKzxKxz = tf.matmul((z_L_inv_Xcov + Xmu[:, None, :]) * eKxz_rbf[..., None], cross_eKzxKxz)  # NxMxM
          return cross_eKzxKxz