docs: def _build_likelihood(self):
          """
          Construct a tensorflow function to compute the bound on the marginal
          likelihood. For a derivation of the terms in here, see the associated
          SGPR notebook.
          """
  
          num_inducing = len(self.feature)
          num_data = tf.cast(tf.shape(self.Y)[0], settings.float_type)
          output_dim = tf.cast(tf.shape(self.Y)[1], settings.float_type)
  
          err = self.Y - self.mean_function(self.X)
          Kdiag = self.kern.Kdiag(self.X)
          Kuf = features.Kuf(self.feature, self.kern, self.X)
          Kuu = features.Kuu(self.feature, self.kern, jitter=settings.numerics.jitter_level)
          L = tf.cholesky(Kuu)
          sigma = tf.sqrt(self.likelihood.variance)
  
          # Compute intermediate matrices
          A = tf.matrix_triangular_solve(L, Kuf, lower=True) / sigma
          AAT = tf.matmul(A, A, transpose_b=True)
          B = AAT + 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
  
          # compute log marginal bound
          bound = -0.5 * num_data * output_dim * np.log(2 * np.pi)
          bound += tf.negative(output_dim) * tf.reduce_sum(tf.log(tf.matrix_diag_part(LB)))
          bound -= 0.5 * num_data * output_dim * tf.log(self.likelihood.variance)
          bound += -0.5 * tf.reduce_sum(tf.square(err)) / self.likelihood.variance
          bound += 0.5 * tf.reduce_sum(tf.square(c))
          bound += -0.5 * output_dim * tf.reduce_sum(Kdiag) / self.likelihood.variance
          bound += 0.5 * output_dim * tf.reduce_sum(tf.matrix_diag_part(AAT))
  
          return bound