docs: def _build_predict(self, Xnew, full_cov=False):
          """
          Compute the mean and variance of the latent function at some new points.
          Note that this is very similar to the SGPR prediction, for which
          there are notes in the SGPR notebook.
          :param Xnew: Point to predict at.
          """
          pX = DiagonalGaussian(self.X_mean, self.X_var)
  
          num_inducing = len(self.feature)
          psi1 = expectation(pX, (self.kern, self.feature))
          psi2 = tf.reduce_sum(expectation(pX, (self.kern, self.feature), (self.kern, self.feature)), axis=0)
          Kuu = features.Kuu(self.feature, self.kern, jitter=settings.numerics.jitter_level)
          Kus = features.Kuf(self.feature, self.kern, Xnew)
          sigma2 = self.likelihood.variance
          sigma = tf.sqrt(sigma2)
          L = tf.cholesky(Kuu)
  
          A = tf.matrix_triangular_solve(L, tf.transpose(psi1), lower=True) / sigma
          tmp = tf.matrix_triangular_solve(L, psi2, lower=True)
          AAT = tf.matrix_triangular_solve(L, tf.transpose(tmp), lower=True) / sigma2
          B = AAT + tf.eye(num_inducing, dtype=settings.float_type)
          LB = tf.cholesky(B)
          c = tf.matrix_triangular_solve(LB, tf.matmul(A, self.Y), lower=True) / sigma
          tmp1 = tf.matrix_triangular_solve(L, Kus, lower=True)
          tmp2 = tf.matrix_triangular_solve(LB, tmp1, lower=True)
          mean = tf.matmul(tmp2, c, transpose_a=True)
          if full_cov:
              var = self.kern.K(Xnew) + tf.matmul(tmp2, tmp2, transpose_a=True) \
                    - tf.matmul(tmp1, tmp1, transpose_a=True)
              shape = tf.stack([1, 1, tf.shape(self.Y)[1]])
              var = tf.tile(tf.expand_dims(var, 2), shape)
          else:
              var = self.kern.Kdiag(Xnew) + tf.reduce_sum(tf.square(tmp2), 0) \
                    - tf.reduce_sum(tf.square(tmp1), 0)
              shape = tf.stack([1, tf.shape(self.Y)[1]])
              var = tf.tile(tf.expand_dims(var, 1), shape)
          return mean + self.mean_function(Xnew), var