5. Linear Least Squares Breakdown and Bayesian Linear Regression

A Motivating Example

Consider the system of equations

with

This is a minimal example from my favourite book on inverse problems, Hansen (2010). While this is a contrived example, it illustrates an very common issue with least squares estimates.

The least squares estimator for this equation is

Now that looks fine at first glance, but what if I told you that the right hand side $b$ was noisy?

Our vector of measurements $b$ had been perturbed by some noise $\varepsilon$, i.e. ${\bf b} = {\bf b}_{true} + \varepsilon$, with

This is not a huge amount of noise, but the least squares solution to ${\bf Ax} = {\bf b}_{true}$ is in fact ${\bf \hat{x}}_{true} = \begin{bmatrix} 1 \\ 1 \end{bmatrix}$! We can modify the original equation a little to take the noise into account:

with $\varepsilon \sim \mathcal N(0, \sigma)$, that is the signal $\bf b$ includes additive white Gaussian noise.

In this example, I knew exactly how much noise was added, so I could recover the true coefficients, but in reality you will not know the exact measurement errors. If you did, they wouldn't be errors! In the best case scenario, you will know the approximate variance of the noise.

How did this precariuous situation come to be? Well, as it turns out, $\bf A$ is what is commonly referred to as ill-conditioned, meaning that (certain) large pertubations in the input $\bf x$ do not change the output $\bf b$ much. And conversely, small pertubations in the measurements can lead to large pertubatons in the estimated coefficients.

Regularisation

A common approach to alleviate the above problem is ridge regression, also sometimes referred to as Tikhonov regularisation. The idea here is to add a shift to the diagonal of the moment matrix ${\bf A}^T {\bf A}$ before inverting it, leading to the ridge estimator

where $\lambda$ is the ridge parameter, which determines the "stength" of regularisation. For example, for $\lambda = 0.001$, we get $\hat{\bf x}_{\lambda=0.001} = \begin{bmatrix} 1.20 \\ 0.70 \end{bmatrix}$, which is much closer to the true parameters than the naïve least squares estimator.

This is equivalent of solving the regularised least squares problem

That is all fine and well, but it leaves an open question of how to choose $\lambda$? In practice, this is often done heuristically, and according to certain well established criteria, as one can study in much more detail in e.g. Hansen (2010). However, these methods still produce a single estimate, rather than a distribution of possible estimates. This is why we will now give this problem the Bayesian treatment.

A Bayesian Approach

To keep this section simple, we are going to make a few relatively strong assumptions. First, we will assume that the noise comes from a zero-mean Gaussian with known variance $\sigma_\varepsilon^2$, $\varepsilon \sim \mathcal N({\bf 0}, \sigma_\varepsilon^2 {\bf I})$. That allows is to write the likelihood as $p({\bf b}|{\bf x}) = \mathcal N({\bf b} - {\bf Ax}, \sigma_\varepsilon^2 {\bf I})$. Furthermore, we assume that the prior distribution of model parameters is also Gaussian, $p({\bf x}) = \mathcal N({\bf \mu}_0, {\bf \Sigma}_0)$. That alows us to write the posterior distribution as (you guessed it!) a Gaussian:

with

Please refer to e.g. Bishop (2006) for more details. For example, if we for the problem outlined above choose $p({\bf x}) = \mathcal N({\bf 0}, {\bf I})$, and $\sigma_\varepsilon^2 = 0.01$, we get

References

Hansen, Per Christian. Discrete Inverse Problems: Insight and Algorithms. Fundamentals of Algorithms. Philadelphia: Society for Industrial and Applied Mathematics, 2010.

Bishop, Christopher M. Pattern Recognition and Machine Learning. Information Science and Statistics. New York: Springer, 2006.