MIT 6.790 — Lecture 5

Bayesian Linear Regression

§ Problems with Linear Regression

Recall that θ^=(X⊤X)−1X⊤Y\hat{\theta} = (X^{\top}X)^{-1}X^{\top}Y is the minimizer of RSS(θ):=(Xθ−Y)⊤(Xθ−Y)\mathrm{RSS}(\theta) := (X\theta - Y)^{\top}(X\theta - Y) (the “Residual Sum of Squares”). One possible problem is that XX might not be full rank; two or more features might be perfectly collinear, making θ^\hat{\theta} non-unique.

Maybe in practice, there's some tiny amount of noise in our observations, so XX is “almost always” guaranteed to be full rank. That doesn't actually fix the issue at hand, though—the effect of that noise is still to effectively pick one of the infinitely-many maximum-likelihood choices for θ^\hat{\theta} at random.

Remark. To clarify, collinearity of features isn't uncommon at all; take “credit card limit” and “credit card rating” as a practical example of two nearly-collinear features.

The issue is that linear regression is overconfident; it always returns a single point estimate θ^\hat{\theta}. Rather than decisively estimating the exact value of θ^\hat{\theta}, perhaps we should take a Bayesian approach: let incoming data {(x(n),y(n))}n=1N\{(x^{(n)}, y^{(n)})\}_{n = 1}^N update our priors on the distribution of θ\theta instead.

§ Bayesian Linear Regression: Theory

Here's the model; refer to Lecture 21 of 18.650 for a reminder on Bayesian inference.

A prior and a model for our data are enough for us to determine the posteriors that result from reading in data.

Theorem. (Linear Regression Posterior) Suppose we collect NN data points, represented by a matrix XX and a column vector YY. Then the posterior is [θ∣X,Y]∼N(μN,ΣN)[\theta \mid X, Y] \sim \mathcal{N}(\mu_N, \Sigma_N), where μN\mu_N and ΣN\Sigma_N satisfy:

ΣN−1=Σ0−1+1σ2X⊤X    and    ΣN−1μN=Σ0−1μ0+1σ2X⊤Y\Sigma_N^{-1} = \Sigma_0^{-1} + \frac{1}{\sigma^2}X^{\top}X ~~~ \text{ and } ~~~ \Sigma_N^{-1}\mu_N = \Sigma_0^{-1}\mu_0 + \frac{1}{\sigma^2}X^{\top}Y

Remark. Goes without saying, but… even though these aren't explicit formulas for ΣN\Sigma_N and μN\mu_N, they do provide enough information to fully compute them if desired.

Proof: Plain computation, with some simplifying tricks. Recall that p(θ∣X,Y)∝p(θ)⋅LN(θ)p(\theta \mid X, Y) \propto p(\theta) \cdot L_N(\theta), where p(θ)p(\theta) is the PDF of a multivariable Gaussian, and LN(θ)L_N(\theta) is a product of NN single-variable Gaussian PDFs.

p(θ)=1(2π)Ddet⁡(Σ0)exp⁡(−12(θ−μ0)⊤Σ0−1(θ−μ0))    and    LN(θ)=1(2πσ2)Nexp⁡(−12σ2∥Y−Xθ∥2).p(\theta) = \frac{1}{\sqrt{(2\pi)^D \det(\Sigma_0)}} \exp\left(-\frac{1}{2}(\theta - \mu_0)^{\top}\Sigma_0^{-1}(\theta - \mu_0)\right) ~~~ \text{ and } ~~~ L_N(\theta) = \frac{1}{\sqrt{(2\pi\sigma^2)^N}} \exp\left(-\frac{1}{2\sigma^2}\|Y - X\theta\|^2\right).

See Lecture 6 of 18.650 for a derivation of the multivariable case.

For ease of computation, it helps to instead work with log⁡p(θ∣X,Y)\log p(\theta \mid X, Y). Also recall that we only care about p(θ∣X,Y)p(\theta \mid X, Y) up to a proportionality constant (as indicated by the ∝\propto symbol), so after taking logs, we can throw away any +C+C terms.

We'll compute the logarithm in two steps, starting with log⁡p(θ)\log p(\theta). (Note that all terms in log⁡p(θ)\log p(\theta) are real numbers.)

log⁡p(θ)=−12(θ−μ0)⊤Σ0−1(θ−μ0)+C=−12[θ⊤Σ0−1θ−μ0⊤Σ0−1θ−θ⊤Σ0−1μ0+μ0⊤Σ0−1μ0]+C=−12[θ⊤Σ0−1θ−μ0⊤Σ0−1θ−θ⊤Σ0−1μ0]+C(★1)=θ⊤(−12Σ0−1)θ+(μ0⊤Σ0−1)θ+C(★2)\begin{align*}\log p(\theta) & = -\frac{1}{2}(\theta - \mu_0)^{\top}\Sigma_0^{-1}(\theta - \mu_0) \textcolor{#008080}{+C} \\ & = -\frac{1}{2}\left[\theta^{\top}\Sigma_0^{-1}\theta - \mu_0^{\top}\Sigma_0^{-1}\theta - \theta^{\top}\Sigma_0^{-1}\mu_0 \textcolor{#008080}{+ \mu_0^{\top}\Sigma_0^{-1}\mu_0}\right] \textcolor{#008080}{+C} \\ & = -\frac{1}{2}\left[\theta^{\top}\Sigma_0^{-1}\theta \textcolor{#bf0140}{- \mu_0^{\top}\Sigma_0^{-1}\theta} \textcolor{#bf0140}{- \theta^{\top}\Sigma_0^{-1}\mu_0}\right] \textcolor{#008080}{+C} & (\bigstar_1) \\ & = \textcolor{blue}{\theta^{\top}\left( -\frac{1}{2}\Sigma_0^{-1} \right)\theta} \textcolor{#bf0140}{+ \left(\mu_0^{\top}\Sigma_0^{-1}\right)\theta} \textcolor{#008080}{+C} & (\bigstar_2)\end{align*}

Two subtle simplifying tricks just took place:

Using some similar tricks, we can compute log⁡LN(θ)\log L_N(\theta) as well.

log⁡LN(θ)=−12σ2(Y−Xθ)⊤(Y−Xθ)+C=−12σ2[Y⊤Y−θ⊤X⊤Y−Y⊤Xθ+θ⊤X⊤Xθ]+C=θ⊤(−12σ2X⊤X)θ+(1σ2Y⊤X)θ+C.\begin{align*}\log L_N(\theta) & = -\frac{1}{2\sigma^2}(Y - X\theta)^{\top}(Y - X\theta) \textcolor{#008080}{+C} \\ & = -\frac{1}{2\sigma^2}\left[ \textcolor{#008080}{Y^{\top}Y} \textcolor{#bf0140}{-\theta^{\top}X^{\top}Y} \textcolor{#bf0140}{-Y^{\top}X\theta} + \theta^{\top}X^{\top}X\theta \right] \textcolor{#008080}{+C} \\ & = \textcolor{blue}{ \theta^{\top}\left( -\frac{1}{2\sigma^2} X^{\top}X \right)\theta } \textcolor{#bf0140}{ +\left(\frac{1}{\sigma^2}Y^{\top}X\right)\theta } \textcolor{#008080}{+C}.\end{align*}

So altogether, we have the following:

(♣1)   log⁡p(θ∣X,Y)=θ⊤(−12[Σ0−1+1σ2X⊤X])θ+(μ0⊤Σ0−1+1σ2Y⊤X)θ+C.(\clubsuit_1) ~~~ \log p(\theta \mid X, Y) = \textcolor{blue}{ \theta^{\top}\left( -\frac{1}{2}\left[ \Sigma_0^{-1} + \frac{1}{\sigma^2}X^{\top}X \right] \right) \theta } \textcolor{#bf0140}{+ \left( \mu_0^{\top}\Sigma_0^{-1} + \frac{1}{\sigma^2}Y^{\top}X \right)\theta} \textcolor{#008080}{+C}.

We just need to compare this with the log of the PDF of N(μN,ΣN)\mathcal{N}(\mu_N, \Sigma_N)… well, we already computed what the logarithm looks like—it's just line (★2)(\bigstar_2) from earlier!

(♣2)   log⁡[PDF of N(μN,ΣN)]=θ⊤(−12ΣN−1)θ+(μN⊤ΣN−1)θ+C(\clubsuit_2) ~~~ \log \left[\text{PDF of } \mathcal{N}(\mu_N, \Sigma_N)\right] = \textcolor{blue}{\theta^{\top}\left( -\frac{1}{2}\Sigma_N^{-1} \right)\theta} \textcolor{#bf0140}{+ \left(\mu_N^{\top}\Sigma_N^{-1}\right)\theta} \textcolor{#008080}{+C}

Pattern-matching (♣1)(\clubsuit_1) against (♣2)(\clubsuit_2) (and taking the transpose of the latter equation…) yields the desired result.   ■\blacksquare

Remark. The constant simplification of everything into the form θ⊤Aθ+bθ+c\theta^{\top}A\theta + b\theta + c isn't a coincidence—every multivariable quadratic polynomial is expressible in such a form! See Lecture 16 of 18.701 for further detail.

§ Bayesian Linear Regression: Performance

The whole point of Bayesian linear regression is to avoid overconfidence. To what extent have we been successful?

Example Set-Up. We draw features x(n)∼Unif[−1,1]x^{(n)} \sim \mathrm{Unif}[-1, 1] and generate responses via y(n)∼N(θ1∗x(n)+θ0∗,σ2)y^{(n)} \sim \mathcal{N}(\theta_1^*x^{(n)} + \theta_0^*, \sigma^2), where θ0∗=−0.3\theta_0^* = -0.3, θ1∗=0.5\theta_1^* = 0.5, and σ=0.2\sigma = 0.2. The algorithm assumes the prior θ∼N(0,I2)\theta \sim \mathcal{N}(0, I_2).

Under this set-up, Bayesian linear regression performs expectedly well; as more data comes in, the posterior becomes more focused around the correct value of θ\theta.

But ordinary linear regression would do just as well here; there is no feature collinearity here.

Altered Set-Up. Everything held constant, except the only feature value that is observed is x(n)=1x^{(n)} = 1.

This is equivalent to feature collinearity. (Recall that, for affine linear regression, all feature vectors have an extra “1” appended to their front to allow for “y-intercepts”; see Lecture 24 of 18.650.) Bayesian linear regression still performs well here!

The regression is able to be confident about θ∗\theta^* in some parameter directions, but not all parameter directions. There's still a one-dimensional space of potential values of θ∗=(θ0∗,θ1∗)\theta^* = (\theta_0^*, \theta_1^*) that the regression sees as possible—still a strong improvement over the two-dimensional prior!

Further Altered Set-Up. The Bayesian regression believes responses are generated via y(n)∼N(θ1x(n)+θ0,σ2)y^{(n)} \sim \mathcal{N}(\theta_1 x^{(n)} + \theta_0, \sigma^2), but they're actually generated by y(n)∼N(0.5−(x(n))2,σ2)y^{(n)} \sim \mathcal{N}(0.5 - (x^{(n)})^2, \sigma^2).

Bayesian regression can't save you from model misspecification. As more data arrives, the posterior becomes more and more resolute in its estimate for θ\theta, even though the regression hasn't come close to actually understanding the response generation.