Probabilistic Artificial Intelligence - Bayesian Linear Regression

Bayesian Linear Regression

Linear Regression

Given a set of (x, y) pairs, linear regression aims to find a linear model that fits the data optimally.
Given the linear model $y=w^Tx$, we want to find optimal weights $w$.
There are many ways of estimating w from data, the most common being the least squares estimator:

$$\hat{w}_\mathrm{ls}\doteq\arg\min_{\boldsymbol{w}\in\mathbb{R}^d}\sum_{i=1}^n(y_i-\boldsymbol{w}^\top\boldsymbol{x}_i)^2=\arg\min_{\boldsymbol{w}\in\mathbb{R}^d}\|\boldsymbol{y}-\boldsymbol{X}\boldsymbol{w}\|_2^2,$$

A slightly different estimator is used for ridge regression:

$$\hat{w}_{\mathrm{ridge}}\doteq\arg\min_{w\in\mathbb{R}^d}\|y-Xw\|_2^2+\lambda\|w\|_2^2$$

As the formula shows, the squared $l_2$ regularization term $\lambda|w|_2^2$ penalizes large $w$ and thus reduces the complexity of the resulting model,so Ridge regression is more robust than standard linear regression in the presence of multicollinearity. Multicollinearity occurs when multiple independent inputs are highly correlated. In this case, their individual effects on the predicted variable cannot be estimated well. Classical linear regression is highly volatile to small input changes. The regularization of ridge regression reduces this volatility by introducing a bias on the weights towards 0.

uncertainty

In practice, our data D is merely a sample of the process we are modeling. In these cases, we are looking for models that generalize to unseen data.
Therefore, it is useful to express the uncertainty about our model due to the lack of data. This uncertainty is commonly referred to as the epistemic uncertainty.
Usually, there is another source of uncertainty called aleatoric uncertainty, which originates directly from the process that we are modeling. This uncertainty is the noise in the labels that cannot be explained by the inputs.

Weight-space View

The most immediate and natural Bayesian interpretation of linear regression is to simply impose a prior on the weights $w$.

Assume that prior $w\sim\mathcal{N}(0,\sigma_p^2I)$ and likelihood $y_i\mid x_i,w\sim\mathcal{N}(w^\top x_i,\sigma_n^2).$ are both Gaussian, we will get the posterior distribution over the weights as:

$$\begin{aligned} &\log p(w\mid x_{1:n},y_{1:n}) \\ &=\log p(w)+\log p(y_{1:n}\mid x_{1:n},w)+\mathrm{const} \\ &=\log p(w)+\sum_{i=1}^n\log p(y_i\mid x_i,w)+\mathrm{const}\\ &=-\frac1{2\sigma_p^2}\|w\|_2^2-\frac1{2\sigma_n^2}\sum_{i=1}^n(y_i-w^\top x_i)^2+\mathrm{const} \\ &=-\frac1{2\sigma_p^2}\|w\|_2^2-\frac1{2\sigma_n^2}\|y-Xw\|_2^2+\mathrm{const} \\ &=-\frac1{2\sigma_p^2}w^\top w-\frac1{2\sigma_n^2}\Big(w^\top X^\top Xw-2y^\top Xw+y^\top y\Big)+\mathrm{const.}\\ &=-\frac12(w-\mu)^\top\Sigma^{-1}(w-\mu)+\mathrm{const} \end{aligned}$$

As the above is a quadratic form in $w$, it follows that the posterior distribution is a Gaussian.
This also shows that Gaussians with known variance and linear likelihood are self-conjugate. A distribution is said to be self-conjugate (or a conjugate prior is self-conjugate) if, when used as a prior distribution, it results in a posterior distribution that belongs to the same family of distributions. It can be shown more generally that Gaussians with known variance are self-conjugate to any Gaussian likelihood. For general distributions the posterior will not be closed-form. This is a very special property of Gaussians.
We can compute the MAP estimate for the weights,

$$\begin{gathered} \hat{w}_{\mathrm{MAP}} =\underset{w}{\operatorname*{\arg\max}}\log p(y_{1:n}\mid x_{1:n},w)+\log p(w) \\ =\arg\min_w\|y-Xw\|_2^2+\frac{\sigma_n^2}{\sigma_p^2}\|w\|_2^2. \end{gathered}$$

we can find that this is simply the MLE loss with an additional $l_2$ regularization term and this coincides with the optimization objective of ridge regression with weight decay $\lambda=\frac{\sigma_n^2}{\sigma_p^2}$. Also, recall that the MAP estimate corresponds to the mode of the posterior distribution, which in the case of a Gaussian is simply its mean.(The mode of a probability distribution is the value where the distribution reaches its maximum point. In the context of the posterior distribution, the mode corresponds to the most probable value of the parameter given the observed data.).

A commonly used alternative to ridge regression is the least absolute shrinkage and selection operator (or lasso), which uses $l_1$ regularization.It turns out that lasso can also be viewed as Bayesian learning, using a Laplace prior.

Note that using point estimates like the MAP estimate does not quantify uncertainty in the weights. The MAP estimate simply collapses all mass of the posterior around its mode.This is especially harmful when we are unsure about the best model.

Recursive Bayesian Updates

As data arrives online (i.e., in “real-time”), we can obtain the new posterior and use it to replace our prior.

Non-linear Regression

We can use linear regression not only to learn linear functions. The trick is to apply a non-linear transformation $\phi:\mathbb{R}^d\to\mathbb{R}^e$.
In Polynomial regression, to learn polynomials of degree m in d input dimensions, we need to apply the nonlinear transformation

$$\begin{aligned} \boldsymbol{\phi}(x)& =[1,x_1,\ldots,x_d,x_1^2,\ldots,x_d^2,x_1\cdot x_2,\ldots,x_{d-1}\cdot x_d, \\ &x_{d-m+1}\cdot\cdot\cdot\cdot\cdot x_d]. \end{aligned}$$

the dimension of the feature space grows exponentially in the degree of polynomials and input dimensions. Even for relatively small m and d, this becomes completely unmanageable.

Function-space View

Previously, we have been interpreting it as a distribution over the weights $w$ of a linear function $\hat{f}=\boldsymbol{\Phi}w.$, we can equivalently consider a distribution directly over the estimated function values. Instead of considering a prior over the weights${w\sim\mathcal{N}}(0,\sigma_p^2I)$, we now impose a prior directly on the values of our model at the observations.Using that Gaussians are closed under linear maps, we obtain the equivalent prior:

$$f\mid X\sim\mathcal{N}(\boldsymbol{\Phi}\mathbb{E}[\boldsymbol{w}],\boldsymbol{\Phi}\mathrm{Var}[\boldsymbol{w}]\boldsymbol{\Phi}^\top)=\mathcal{N}(\boldsymbol{0},\underbrace{\sigma_p^2\boldsymbol{\Phi}\boldsymbol{\Phi}^\top}_{\boldsymbol{\kappa}})$$

K is the kernel matrix.The kernel function is:
$k(x,x^{\prime})\doteq\sigma_p^2\cdot\phi(x)^\top\phi(x^{\prime})$
The kernel matrix is a covariance matrix and the kernel function measures the covariance of the function values $k(x,x^{\prime})=\mathrm{Cov}\big[f(x),f(x^{\prime})\big].$
Moreover, note that we have reformulated the learning algorithm such that the feature space is now implicit in the choice of kernel, and the kernel is defined by inner products of (nonlinearly transformed) inputs. In other words, the choice of kernel implicitly determines the class of functions that f is sampled from, and encodes our prior beliefs. This is known as the kernel trick.

reference

[1] A. Krause, “Probabilistic Artificial Intelligence”.