Gaussian Process Classification Techniques

While Gaussian Processes (GPs) provide an elegant framework for regression tasks with built-in uncertainty quantification, applying them to classification requires navigating a significant hurdle: the likelihood function is no longer Gaussian. In binary classification, we typically deal with Bernoulli likelihoods (outputting 0 or 1), and for multi-class problems, categorical likelihoods. These step away from the convenient conjugate relationship between Gaussian priors and Gaussian likelihoods that makes GP regression analytically tractable.

The Latent Function Approach

The standard way to tackle GP classification is by introducing a latent function, $f(\mathbf{x})$, which follows a Gaussian Process prior:

$$f(\mathbf{x}) \sim \mathcal{GP}(m(\mathbf{x}), k(\mathbf{x}, \mathbf{x}'))$$

Here, $m(\mathbf{x})$ is the mean function (often set to zero) and $k(\mathbf{x}, \mathbf{x}')$ is the covariance (kernel) function, just as in GP regression.

Instead of this latent function directly representing the output, it's passed through a link function (or inverse link function in GLM terminology, often called a "squashing" function here) to produce probabilities for the class labels. For binary classification ($y \in \{0, 1\}$), a common choice is the logistic sigmoid function, $\sigma(z) = 1 / (1 + e^{-z})$:

$$p(y=1 | f(\mathbf{x})) = \sigma(f(\mathbf{x}))$$

Alternatively, the probit link function (the cumulative density function of a standard normal distribution) can be used.

The full model involves the GP prior over the latent function $f$ and the non-Gaussian likelihood for the observed labels $\mathbf{y}$ given the latent function values $\mathbf{f} = [f(\mathbf{x}_1), ..., f(\mathbf{x}_N)]^T$ at the training points $\mathbf{X}$:

$$p(\mathbf{y} | \mathbf{f}) = \prod_{i=1}^N p(y_i | f(\mathbf{x}_i))$$

For the logistic sigmoid link, this likelihood is a product of Bernoulli probabilities: $p(y_i | f(\mathbf{x}_i)) = \sigma(f(\mathbf{x}_i))^{y_i} (1 - \sigma(f(\mathbf{x}_i)))^{1-y_i}$.

The Intractability Challenge

The core difficulty arises when we try to compute the posterior distribution over the latent function, $p(\mathbf{f} | \mathbf{y}, \mathbf{X})$. According to Bayes' theorem:

$$p(\mathbf{f} | \mathbf{y}, \mathbf{X}) = \frac{p(\mathbf{y} | \mathbf{f}) p(\mathbf{f} | \mathbf{X})}{p(\mathbf{y} | \mathbf{X})}$$

Here, $p(\mathbf{f} | \mathbf{X})$ is the GP prior (multivariate Gaussian), but $p(\mathbf{y} | \mathbf{f})$ is non-Gaussian (e.g., product of Bernoullis). The product of a Gaussian and a non-Gaussian term generally does not result in a Gaussian distribution. This means the posterior $p(\mathbf{f} | \mathbf{y}, \mathbf{X})$ is non-Gaussian and intractable to compute analytically.

Furthermore, the marginal likelihood (or evidence), $p(\mathbf{y} | \mathbf{X})$, which is essential for model comparison and hyperparameter optimization, also becomes intractable:

$$p(\mathbf{y} | \mathbf{X}) = \int p(\mathbf{y} | \mathbf{f}) p(\mathbf{f} | \mathbf{X}) d\mathbf{f}$$

This integral involves integrating a non-Gaussian likelihood against a high-dimensional Gaussian prior, which usually lacks a closed-form solution.

Approximation Strategies

To perform inference and make predictions, we need to approximate the intractable posterior $p(\mathbf{f} | \mathbf{y}, \mathbf{X})$. The most common techniques are the Laplace Approximation and Expectation Propagation (EP).

Laplace Approximation

The Laplace approximation finds the mode of the posterior distribution and then fits a Gaussian distribution centered at that mode. The covariance of this Gaussian is determined by the curvature of the log posterior at the mode.

1. Find the Mode: We maximize the log posterior probability with respect to $\mathbf{f}$:

   $$\log p(\mathbf{f} | \mathbf{y}, \mathbf{X}) = \log p(\mathbf{y} | \mathbf{f}) + \log p(\mathbf{f} | \mathbf{X}) - \log p(\mathbf{y} | \mathbf{X})$$

   Since $\log p(\mathbf{y} | \mathbf{X})$ is constant with respect to $\mathbf{f}$, we maximize $\Psi(\mathbf{f}) = \log p(\mathbf{y} | \mathbf{f}) + \log p(\mathbf{f} | \mathbf{X})$. This is equivalent to finding the Maximum A Posteriori (MAP) estimate of $\mathbf{f}$, denoted $\mathbf{f}_{MAP}$. Because $\log p(\mathbf{f} | \mathbf{X})$ is quadratic in $\mathbf{f}$ (from the Gaussian prior) and $\log p(\mathbf{y} | \mathbf{f})$ is typically concave (e.g., for logistic/probit likelihoods), this optimization problem is often convex and can be solved efficiently using methods like Newton-Raphson.

2. Fit Gaussian: We approximate the posterior with a Gaussian $q(\mathbf{f} | \mathbf{y}, \mathbf{X}) = \mathcal{N}(\mathbf{f} | \boldsymbol{\mu}_{LA}, \boldsymbol{\Sigma}_{LA})$, where:

   • Mean $\boldsymbol{\mu}_{LA} = \mathbf{f}_{MAP}$.
   • Covariance $\boldsymbol{\Sigma}_{LA} = (-\nabla \nabla \Psi(\mathbf{f}_{MAP}))^{-1}$. The precision matrix is the negative Hessian of the log posterior evaluated at the mode.

The Laplace approximation provides a Gaussian approximation $q(\mathbf{f})$ to the true posterior $p(\mathbf{f} | \mathbf{y}, \mathbf{X})$, making subsequent calculations tractable. However, its accuracy depends on how well the true posterior is represented by a Gaussian centered at its mode. If the posterior is significantly skewed or multi-modal, the approximation might be poor.

Expectation Propagation (EP)

Expectation Propagation (EP) offers an alternative, often more accurate, approach. It seeks to approximate the posterior $p(\mathbf{f} | \mathbf{y}, \mathbf{X})$ with a distribution $q(\mathbf{f})$ from a chosen family (typically Gaussian), by minimizing the KL divergence $KL(p || q)$ in a specific way.

EP approximates the full posterior by iteratively refining an approximation based on incorporating the likelihood terms one by one. It tries to match certain moments (usually the mean and variance) of the approximate marginal posterior distributions with those of hypothetical distributions obtained by combining the prior with individual likelihood terms.

1. Initialize: Start with an approximate posterior $q(\mathbf{f}) \propto p(\mathbf{f})$. Often, $q(\mathbf{f})$ is initialized to the prior $p(\mathbf{f})$. The approximation $q(\mathbf{f})$ is typically maintained in the exponential family form, often as a Gaussian $\mathcal{N}(\boldsymbol{\mu}_{EP}, \boldsymbol{\Sigma}_{EP})$.
2. Iterative Refinement: EP iteratively updates the parameters of $q(\mathbf{f})$ by considering each likelihood term $p(y_i | f_i)$ (called a "site"). For each site $i$:
   • Form a "cavity distribution" by removing the influence of the current approximation of site $i$ from $q(\mathbf{f})$.
   • Form a "tilted distribution" by multiplying the cavity distribution with the exact likelihood term $p(y_i | f_i)$. This tilted distribution is generally non-Gaussian.
   • Find the parameters of a new Gaussian site approximation that best matches the moments (mean and variance) of the tilted distribution. This often involves computing intractable integrals, which might themselves need approximation, especially for the logistic likelihood (though tractable for the probit likelihood).
   • Update the overall posterior approximation $q(\mathbf{f})$ using the new site approximation.
3. Convergence: Repeat the iterative refinement process until the parameters of $q(\mathbf{f})$ converge.

EP generally provides a more flexible and often more accurate approximation than Laplace, particularly when the true posterior deviates significantly from a simple Gaussian shape around the mode. However, EP is computationally more intensive, its convergence is not always guaranteed, and the moment-matching integrals can be challenging.

Making Predictions

Once we have a Gaussian approximation to the posterior over the latent values at the training points, $q(\mathbf{f} | \mathbf{y}, \mathbf{X}) \approx \mathcal{N}(\boldsymbol{\mu}_{approx}, \boldsymbol{\Sigma}_{approx})$ (where 'approx' is LA or EP), we can derive the approximate predictive distribution for the latent function value $f_*$ at a new test point $\mathbf{x}_*$.

Similar to GP regression, the joint distribution of $\mathbf{f}$ and $f_*$ under the prior is Gaussian. Using the properties of conditioning Gaussian distributions, the approximate predictive distribution $q(f_* | \mathbf{y}, \mathbf{X})$ is also Gaussian:

$$q(f_* | \mathbf{y}, \mathbf{X}) = \mathcal{N}(f_* | \mu_*, \sigma_*^2)$$

where $\mu_*$ and $\sigma_*^2$ depend on the approximate posterior mean $\boldsymbol{\mu}_{approx}$, covariance $\boldsymbol{\Sigma}_{approx}$, the kernel evaluations involving $\mathbf{x}_*$, and the training data $\mathbf{X}$.

To get the final predictive probability for the class label $y_*=1$, we need to average the likelihood $\sigma(f_*)$ over this approximate predictive distribution:

$$p(y_*=1 | \mathbf{x}_*, \mathbf{y}, \mathbf{X}) \approx \int \sigma(f_*) q(f_* | \mathbf{y}, \mathbf{X}) df_*$$

This integral is still often intractable (e.g., integrating the logistic sigmoid over a Gaussian). Common ways to compute it include:

• Numerical quadrature (e.g., Gauss-Hermite quadrature).
• Monte Carlo estimation: Sample $f_*$ values from $q(f_* | \mathbf{y}, \mathbf{X})$, pass them through $\sigma(\cdot)$, and average the results.
• Analytical approximations (possible for the probit link, where the integral involves the CDF of another Gaussian).

The resulting predictive probability incorporates the uncertainty captured by the approximate posterior distribution of $f_*$.

Example visualization of a GP classification result in 1D. The dashed line shows the mean of the approximate posterior latent function $f(x)$. The shaded region represents the uncertainty (e.g., +/- 2 standard deviations) in $f(x)$. Data points from two classes are shown. The probability of a point belonging to Class 1 is estimated by applying the sigmoid function to $f(x)$. The decision boundary occurs where the mean latent function crosses zero (or where the sigmoid output is 0.5). Uncertainty in $f(x)$ translates to uncertainty in the classification, especially near the decision boundary.

Hyperparameter Optimization

Just like in GP regression, the kernel hyperparameters (e.g., lengthscale, variance) and any parameters of the likelihood function need to be determined. This is typically done by maximizing an approximation to the log marginal likelihood $\log p(\mathbf{y} | \mathbf{X}, \boldsymbol{\theta})$, where $\boldsymbol{\theta}$ denotes the hyperparameters. Both the Laplace approximation and EP provide estimates of this quantity as byproducts of their approximation process. Gradient-based optimization methods are commonly used to find the hyperparameters $\boldsymbol{\theta}$ that maximize this approximate marginal likelihood.

In summary, GP classification extends the core ideas of GPs to scenarios with non-Gaussian likelihoods by introducing a latent function and employing approximation techniques like Laplace or EP to handle the resulting analytical intractability. These methods allow us to obtain approximate posterior distributions and predictive probabilities, retaining the Bayesian approach's ability to quantify uncertainty even in classification settings. The choice between Laplace and EP often involves a trade-off between computational cost and the accuracy of the posterior approximation.