docs: def _build_likelihood(self):
          """
          Construct a tensorflow function to compute the bound on the marginal
          likelihood.
          """
  
          # FITC approximation to the log marginal likelihood is
          # log ( normal( y | mean, K_fitc ) )
          # where K_fitc = Qff + diag( \nu )
          # where Qff = Kfu Kuu^{-1} Kuf
          # with \nu_i = Kff_{i,i} - Qff_{i,i} + \sigma^2
  
          # We need to compute the Mahalanobis term -0.5* err^T K_fitc^{-1} err
          # (summed over functions).
  
          # We need to deal with the matrix inverse term.
          # K_fitc^{-1} = ( Qff + \diag( \nu ) )^{-1}
          #            = ( V^T V + \diag( \nu ) )^{-1}
          # Applying the Woodbury identity we obtain
          #            = \diag( \nu^{-1} ) - \diag( \nu^{-1} ) V^T ( I + V \diag( \nu^{-1} ) V^T )^{-1) V \diag(\nu^{-1} )
          # Let \beta =  \diag( \nu^{-1} ) err
          # and let \alpha = V \beta
          # then Mahalanobis term = -0.5* ( \beta^T err - \alpha^T Solve( I + V \diag( \nu^{-1} ) V^T, alpha ) )
  
          err, nu, Luu, L, alpha, beta, gamma = self._build_common_terms()
  
          mahalanobisTerm = -0.5 * tf.reduce_sum(tf.square(err) / tf.expand_dims(nu, 1)) \
                            + 0.5 * tf.reduce_sum(tf.square(gamma))
  
          # We need to compute the log normalizing term -N/2 \log 2 pi - 0.5 \log \det( K_fitc )
  
          # We need to deal with the log determinant term.
          # \log \det( K_fitc ) = \log \det( Qff + \diag( \nu ) )
          #                    = \log \det( V^T V + \diag( \nu ) )
          # Applying the determinant lemma we obtain
          #                    = \log [ \det \diag( \nu ) \det( I + V \diag( \nu^{-1} ) V^T ) ]
          #                    = \log [ \det \diag( \nu ) ] + \log [ \det( I + V \diag( \nu^{-1} ) V^T ) ]
  
          constantTerm = -0.5 * self.num_data * tf.log(tf.constant(2. * np.pi, settings.float_type))
          logDeterminantTerm = -0.5 * tf.reduce_sum(tf.log(nu)) - tf.reduce_sum(tf.log(tf.matrix_diag_part(L)))
          logNormalizingTerm = constantTerm + logDeterminantTerm
  
          return mahalanobisTerm + logNormalizingTerm * self.num_latent