docs: def _build_likelihood(self):
          r"""
          q_alpha, q_lambda are variational parameters, size N x R
          This method computes the variational lower bound on the likelihood,
          which is:
              E_{q(F)} [ \log p(Y|F) ] - KL[ q(F) || p(F)]
          with
              q(f) = N(f | K alpha + mean, [K^-1 + diag(square(lambda))]^-1) .
          """
          K = self.kern.K(self.X)
          K_alpha = tf.matmul(K, self.q_alpha)
          f_mean = K_alpha + self.mean_function(self.X)
  
          # compute the variance for each of the outputs
          I = tf.tile(tf.expand_dims(tf.eye(self.num_data, dtype=settings.float_type), 0),
                      [self.num_latent, 1, 1])
          A = I + tf.expand_dims(tf.transpose(self.q_lambda), 1) * \
              tf.expand_dims(tf.transpose(self.q_lambda), 2) * K
          L = tf.cholesky(A)
          Li = tf.matrix_triangular_solve(L, I)
          tmp = Li / tf.expand_dims(tf.transpose(self.q_lambda), 1)
          f_var = 1. / tf.square(self.q_lambda) - tf.transpose(tf.reduce_sum(tf.square(tmp), 1))
  
          # some statistics about A are used in the KL
          A_logdet = 2.0 * tf.reduce_sum(tf.log(tf.matrix_diag_part(L)))
          trAi = tf.reduce_sum(tf.square(Li))
  
          KL = 0.5 * (A_logdet + trAi - self.num_data * self.num_latent +
                      tf.reduce_sum(K_alpha * self.q_alpha))
  
          v_exp = self.likelihood.variational_expectations(f_mean, f_var, self.Y)
          return tf.reduce_sum(v_exp) - KL