docs: def _build_predict(self, Xnew, full_cov=False):
          """
          Compute the mean and variance of the latent function at some new points
          Xnew. For a derivation of the terms in here, see the associated SGPR
          notebook. 
          """
          num_inducing = len(self.feature)
          err = self.Y - self.mean_function(self.X)
          Kuf = features.Kuf(self.feature, self.kern, self.X)
          Kuu = features.Kuu(self.feature, self.kern, jitter=settings.numerics.jitter_level)
          Kus = features.Kuf(self.feature, self.kern, Xnew)
          sigma = tf.sqrt(self.likelihood.variance)
          L = tf.cholesky(Kuu)
          A = tf.matrix_triangular_solve(L, Kuf, lower=True) / sigma
          B = tf.matmul(A, A, transpose_b=True) + tf.eye(num_inducing, dtype=settings.float_type)
          LB = tf.cholesky(B)
          Aerr = tf.matmul(A, err)
          c = tf.matrix_triangular_solve(LB, Aerr, 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)
              var = tf.tile(var[None, ...], [self.num_latent, 1, 1])  # P x N x N
          else:
              var = self.kern.Kdiag(Xnew) + tf.reduce_sum(tf.square(tmp2), 0) \
                    - tf.reduce_sum(tf.square(tmp1), 0)
              var = tf.tile(var[:, None], [1, self.num_latent])
          return mean + self.mean_function(Xnew), var