docs: def compute_qu(self):
          """
          Computes the mean and variance of q(u), the variational distribution on
          inducing outputs. SVGP with this q(u) should predict identically to
          SGPR.
          :return: mu, A
          """
          Kuf = features.Kuf(self.feature, self.kern, self.X)
          Kuu = features.Kuu(self.feature, self.kern, jitter=settings.jitter)
  
          Sig = Kuu + (self.likelihood.variance ** -1) * tf.matmul(Kuf, Kuf, transpose_b=True)
          Sig_sqrt = tf.cholesky(Sig)
  
          Sig_sqrt_Kuu = tf.matrix_triangular_solve(Sig_sqrt, Kuu)
  
          A = tf.matmul(Sig_sqrt_Kuu, Sig_sqrt_Kuu, transpose_a=True)
          mu = tf.matmul(Sig_sqrt_Kuu,
                         tf.matrix_triangular_solve(Sig_sqrt, tf.matmul(Kuf, self.Y - self.mean_function(self.X))),
                         transpose_a=True) * self.likelihood.variance ** -1.0
  
          return mu, A