Evidence approximation in linear regression: A method that produces “automatically regularized” solutions.


In this post, I look at a Bayesian treatment of the linear regression problem. Making use of basis functions allows you to model non-linear patterns in data, however taking this route usually requires that you regularize your solution. To find the best regularization parameter often requires cross validation, but by looking at a framework known as Evidence Approximation, this allows us to reformulate regression optimization into a form where we can analytically find the most likely value for the regularization parameter. What results is a lightweight algorithm that requires no parameter-tuning.

(Note: This post is based on my collection of reading notes from Section 3.5 of Chris Bishop’s book “Pattern Recognition and Machine Learning”. I highly recommend this text, and encourage you to look at this source for derivations and deeper insights.)

Some theory…

We start with a linear regression problem which has the form,

[latex] y(\mathbf{x}, \mathbf{w}) = \mathbf{w}^T \boldsymbol\phi(\mathbf{x}) [/latex]

where [latex]\boldsymbol\phi(\mathbf{x})[/latex] is a vector of basis function values. Using linear basis functions is powerful, as it lets you model non-linear relationships in your data. Unfortunately, selecting the best basis functions is an art and often times you need to take a “kitchen sink” approach and throw in as many basis functions as possible.

As soon as you do this though you run into overfitting problems as well as numerical problems as your design matrix can become ill-conditioned. The solution to this, of course, is to add a regularization penalty on $latex \mathbf{w}$, and the common practice is to use cross validation to find the right value of $latex \alpha$.

It turns out that by taking a Bayesian look at the linear regression problem, we can derive a solution for the value of the regularization parameter that does not require cross validation, rather it’s easily estimated from your training data. Let’s dig into this:

We start with a generative model that measures the likelihood of your observations $latex \{\mathbf{t},\mathbf{x}\}$,

[latex]p(\mathbf{t}|\mathbf{w},\alpha,\beta)=p(\mathbf{t}|\boldsymbol\Phi^T\mathbf{w},\beta) p(\mathbf{w}|\mathbf{0},\alpha\mathbf{I})[/latex]

Where the probability density functions on the right are both Gaussian, [latex]\boldsymbol\Phi[/latex] is the design matrix of basis functions evaluated at your input points [latex]\mathbf{x}[/latex], and $latex \beta$ and $latex \alpha$ are precision parameters for their respective distributions. Note that the last distribution on the r.h.s. is a prior on the distribution of weights $latex \mathbf{w}$ with zero mean and precision $latex \alpha \mathbf{I}$.

Our goal is to find the most likely value for [latex]\alpha[/latex] and [latex]\beta[/latex] given the training data, this requires that we integrate out $latex \mathbf{w}$. Hence, we have,

[latex] p(\mathbf{t}|\alpha,\beta ) =\int p(\mathbf{t} | \boldsymbol\Phi^T\mathbf{w}, \beta) p(\mathbf{w}|\alpha \mathbf{I})\delta\mathbf{w}[/latex]. (1)

As we’ll see shortly, this solution leads to an optimally regularized linear regression solution.

We can re-write Equation 1 as,

[latex](\frac{\beta}{2\pi})^{N/2} (\frac{\alpha}{2\pi})^{M/2}\int-E(\mathbf{w})\delta\mathbf{w}[/latex]

where [latex]E(\mathbf{w}) = \frac{\beta}{2}||\mathbf{t} – \boldsymbol\Phi^T\mathbf{w}||^2 + \frac{\alpha}{2}\mathbf{w}^T\mathbf{w} [/latex], N is the number of training points and M is the number of basis functions.

Here is where the connection to regularized linear regression comes out. Notice that [latex]E( \mathbf{w} )[/latex] is the same the regularized least squares problem where the regularization parameter [latex]\lambda = \alpha / \beta[/latex]. Also note that optimization of the log-likelihood of (1) w.r.t $latex \mathbf{w}$ reduces to the optimization of $latex E( \mathbf{w} )$.

Our goal at this point is to find the most likely [latex]\alpha[/latex] and [latex]\beta[/latex] within [latex]p(\mathbf{t}|\alpha,\beta )[/latex]. This gives us an implicit value for the regularization parameter, at which point we are free to estimate the most likely value for $latex \mathbf{w}$ using the least squares solution. There are a few methods for doing this, for example, EM or trying to find a direct solution. In this post we look at the direct solution. You can follow the derivation in Section 3.5 of Bishop’s text. The solutions for [latex]\alpha[/latex] and [latex]\beta[/latex] are,


[latex]\beta=\frac{1}{N-\gamma}||\mathbf{t}-\boldsymbol\Phi \mathbf{m}||^2[/latex]

where [latex]\gamma=\sum{\frac{\lambda_i}{\alpha+\lambda_i}}[/latex] where the  [latex]\lambda_i[/latex]’s are the eigenvalues of the matrix [latex]\beta\boldsymbol\Phi^T\boldsymbol\Phi[/latex], [latex]\mathbf{m}=\beta \mathbf{A}^{-1} \boldsymbol\Phi^T\mathbf{t}[/latex], and [latex]\mathbf{A}=\alpha\mathbf{I}+\beta \boldsymbol\Phi^T\boldsymbol\Phi[/latex].

These solutions are implicit, but we can estimate [latex]\alpha[/latex] and [latex]\beta[/latex] by plugging in initial estimates and iterating over the recurrence relations until convergence.

Testing it out…

I wrote a Matlab script implementing the solution (You can download the code below). Again, the nice thing about this solution is that all you need to provide is the training data and your design matrix and an optimally “regularized” solution just pops out.

I tested this method out on the classic “Series G”  data set of Box and Jenkins  which represents airline passenger totals from 1949 to 1960. The dataset is displayed as the solid blue curve in the image above. Noticing the obviously periodic trends I created a design matrix that contained sinusoids with various periods and offsets. Specifically I added sinusoids [latex]\phi_i(x) = sin(ax + b) [/latex] for[latex]a[/latex] between .25 and 10 in steps of .25 and [latex]b[/latex] between 0 and [latex]\pi[/latex] in steps of [latex]\pi/4[/latex]. Also, noticing the general uptrend in the data I added a single linear function to the design matrix.

The results are satisfying. The method is able to find the correct weights for [latex]\alpha[/latex]  and [latex]\beta[/latex], such that it picks up the periodic uptrend, all without needing to resort to cross validation to tune parameters.


  • Download a Matlab implementation of the code.

Leave a Reply

Your email address will not be published. Required fields are marked *