NO.30.tip： 线性回归

背景

线性回归(linear regression)是广义线性模型最简单的情况。线性回归指的是以下模型：

$$p(y|\mathbf{x},\mathbf{\theta })=\mathcal{N}(y|w_0+{\mathbf{w}}^{\top }\mathbf{x},{\sigma }^2)\text{\hspace{1em}(1)}$$

其中，$\mathbf{\theta }=(w_0,\mathbf{w},{\sigma }^2)$是模型的所有参数。（在统计学中，参数$w_0$和$\mathbf{w}$通常使用${\beta }_0$和$\mathbf{\beta }$表示。）

我们将简要讨论最大似然估计，然后重点讨论贝叶斯分析。

源与流

普通最小二乘法

根据式(1)，我们可以得出数据的负对数似然，如下所示：

$$\text{NLL}(\mathbf{w},{\sigma }^2)=-\sum _{n=1}^N\log \left[{\left(\frac{1}{2\pi {\sigma }^2}\right)}^{\frac{1}{2}}\exp \left(-\frac{1}{2{\sigma }^2}(y_n-{\mathbf{w}}^{\top }{\mathbf{x}}_n{)}^2\right)\right]\text{\hspace{1em}(2)}$$

$$=\frac{1}{2{\sigma }^2}\sum _{n=1}^N(y_n-{\hat{y}}_n{)}^2+\frac{N}{2}\log (2\pi {\sigma }^2)\text{\hspace{1em}(3)}$$

其中，我们定义了预测响应 ${\hat{y}}_n\triangleq {\mathbf{w}}^{\top }{\mathbf{x}}_n$。我们证明了最大似然估计由下式给出：

$${\hat{\mathbf{w}}}_{\text{mle}}=({\mathbf{X}}^{\top }\mathbf{X}{)}^{-1}{\mathbf{X}}^{\top }\mathbf{y}\text{\hspace{1em}(4)}$$

这被称为普通最小二乘(Ordinary Least Squares, OLS)解。

观测噪声的最大似然估计由下式给出：

$${\hat{\sigma }}_{\text{mle}}^2=\arg \underset{{\sigma }^2}{\min }\text{NLL}(\hat{\mathbf{w}},{\sigma }^2)=\frac{1}{N}\sum _{n=1}^N(y_n-{\mathbf{x}}_n^{\top }\hat{\mathbf{w}}{)}^2\text{\hspace{1em}(5)}$$

这正是残差的均方误差，这是一个直观的结果。

共轭先验

我们使用共轭先验推导参数的后验。首先考虑只有$\mathbf{w}$未知的情况（因此观测噪声方差参数${\sigma }^2$是固定的），然后考虑一般情况，其中${\sigma }^2$和$\mathbf{w}$均是未知的。

噪声方差已知的情况

线性回归的共轭先验具有以下形式：

$$p(\mathbf{w})=\mathcal{N}(\mathbf{w}|\mathbf{w},\mathbf{\Sigma })\text{\hspace{1em}(6)}$$

我们经常使用$\mathbf{w}=0$作为先验均值，$\mathbf{\Sigma }={\tau }^2{\mathbf{I}}_D$作为先验协方差。（假设偏置项包含在权重向量中，但通常使用较弱的先验，因为我们通常不想对输出的总体平均水平进行正则化。）

为了推导后验，让我们首先根据多元高斯分布重写似然，如下所示：

$$\ell (\mathbf{w})=p(\mathcal{D}|\mathbf{w},{\sigma }^2)=\prod _{n=1}^{N}p(y_n|{\mathbf{w}}^{\top }\mathbf{x},{\sigma }^2)=\mathcal{N}(\mathbf{y}|\mathbf{Xw},{\sigma }^2{\mathbf{I}}_N)\text{\hspace{1em}(7)}$$

其中，${\mathbf{I}}_N$是$N\times N$单位矩阵。然后，可以使用高斯的贝叶斯规则来推导出后验，如下所示：

$$p(\mathbf{w}|\mathbf{X},\mathbf{y},{\sigma }^2)\propto \mathcal{N}(\mathbf{w}|\mathbf{w},\mathbf{\Sigma })\mathcal{N}(\mathbf{y}|\mathbf{Xw},{\sigma }^2{\mathbf{I}}_N)=\mathcal{N}(\mathbf{w}|\hat{\mathbf{w}},\hat{\mathbf{\Sigma }})\text{\hspace{1em}(8)}$$

$$\hat{\mathbf{w}}\triangleq \hat{\mathbf{\Sigma }}\left({\mathbf{\Sigma }}^{-1}\mathbf{w}+\frac{1}{{\sigma }^2}{\mathbf{X}}^{\top }\mathbf{y}\right)\text{\hspace{1em}(9)}$$

$$\hat{\mathbf{\Sigma }}\triangleq {\left({\tilde{\mathbf{\Sigma }}}^{-1}+\frac{1}{{\sigma }^2}{\mathbf{X}}^{\top }\mathbf{X}\right)}^{-1}\text{\hspace{1em}(10)}$$

其中，$\hat{\mathbf{w}}$是后验均值，$\hat{\mathbf{\Sigma }}$是后验协方差。

现在，假设$\mathbf{w}=0$，$\mathbf{\Sigma }={\tau }^2\mathbf{I}$。在这种情况下，后验均值变为

$$\hat{\mathbf{w}}=\frac{1}{{\sigma }^2}\hat{\mathbf{\Sigma }}{\mathbf{X}}^{\top }\mathbf{y}={\left(\frac{{\sigma }^2}{{\tau }^2}\mathbf{I}+{\mathbf{X}}^{\top }\mathbf{X}\right)}^{-1}{\mathbf{X}}^{\top }\mathbf{y}\text{\hspace{1em}(11)}$$

如果我们定义$\lambda =\frac{{\sigma }^2}{{\tau }^2}$，就会发现这相当于岭回归，用于优化下式：

$$\mathcal{L}(\mathbf{w})=\text{RSS}(\mathbf{w})+\lambda \Vert \mathbf{w}{\Vert }^2\text{\hspace{1em}(12)}$$

其中，RSS是残差平方和：

$$\text{RSS}(\mathbf{w})=\frac{1}{2}\sum _{n=1}^N(y_n-{\mathbf{w}}^{\top }{\mathbf{x}}_n{)}^2=\frac{1}{2}\Vert \mathbf{Xw}-\mathbf{y}{\Vert }_2^2=\frac{1}{2}(\mathbf{Xw}-\mathbf{y}{)}^{\top }(\mathbf{Xw}-\mathbf{y})\text{\hspace{1em}(13)}$$

噪声方差未知的情况

我们假设$\mathbf{w}$和${\sigma }^2$都是未知的。似然由下式给出：

$$\ell (\mathbf{w},{\sigma }^2)=p(\mathcal{D}|\mathbf{w},{\sigma }^2)\propto ({\sigma }^2{)}^{-N/2}\exp \left(-\frac{1}{2{\sigma }^2}\sum _{n=1}^N(y_n-{\mathbf{w}}^{\top }{\mathbf{x}}_n{)}^2\right)\text{\hspace{1em}(14)}$$

由于回归权重现在取决于似然中的${\sigma }^2$，因此$\mathbf{w}$的共轭先验具有以下形式：

$$p(\mathbf{w}|{\sigma }^2)=\mathcal{N}(\mathbf{w}|\mathbf{w},{\sigma }^2\mathbf{\Sigma })\text{\hspace{1em}(15)}$$

对于噪声方差${\sigma }^2$，共轭先验基于逆伽马分布，其形式为：

$$\text{IG}({\sigma }^2|a,b)=\frac{b^a}{\Gamma (a)}({\sigma }^2{)}^{-(a+1)}\exp \left(-\frac{\tilde{b}}{{\sigma }^2}\right)\text{\hspace{1em}(16)}$$

将这两者结合在一起，我们发现联合共轭先验是正态逆伽马分布：

$$\text{NIG}(\mathbf{w},{\sigma }^2|\mathbf{w},\mathbf{\Sigma },a,b)\triangleq \mathcal{N}(\mathbf{w}|\mathbf{w},{\sigma }^2\mathbf{\Sigma })\text{IG}({\sigma }^2|a,b)\text{\hspace{1em}(17)}$$

$$=\frac{b^a}{(2\pi {)}^{D/2}|\mathbf{\Sigma }{|}^{\frac{1}{2}}\Gamma (a)}({\sigma }^2{)}^{-(a+(D/2)+1)}\times \exp \left[-\frac{(\mathbf{w}-\mathbf{w}{)}^{\top }{\mathbf{\Sigma }}^{-1}(\mathbf{w}-\mathbf{w})+2\tilde{b}}{2{\sigma }^2}\right]\text{\hspace{1em}(18)}$$

这产生以下的后验：

$$p(\mathbf{w},{\sigma }^2|\mathcal{D})=\text{NIG}(\mathbf{w},{\sigma }^2|\hat{\mathbf{w}},\hat{\mathbf{\Sigma }},\hat{a},\hat{b})\text{\hspace{1em}(19)}$$

$$\hat{\mathbf{w}}=\hat{\mathbf{\Sigma }}({\mathbf{\Sigma }}^{-1}\mathbf{w}+{\mathbf{X}}^{\top }\mathbf{y})\text{\hspace{1em}(20)}$$

$$\hat{\mathbf{\Sigma }}=({\tilde{\mathbf{\Sigma }}}^{-1}+{\mathbf{X}}^{\top }\mathbf{X}{)}^{-1}\text{\hspace{1em}(21)}$$

$$\hat{a}=\tilde{a}+N/2\text{\hspace{1em}(22)}$$

$$\hat{b}=b+\frac{1}{2}({\mathbf{w}}^{\top }{\mathbf{\Sigma }}^{-1}\mathbf{w}+{\mathbf{y}}^{\top }\mathbf{y}-{\hat{\mathbf{w}}}^{\top }{\hat{\mathbf{\Sigma }}}^{-1}\hat{\mathbf{w}})\text{\hspace{1em}(23)}$$

$\hat{\mathbf{w}}$和$\hat{\mathbf{\Sigma }}$的表达式类似于${\sigma }^2$已知的情况。$\hat{a}$的表达式也是直观的，因为该表达式只是更新计数。$\hat{b}$的表达式可以解释如下：该表达式是先验平方和$\tilde{b}$，加上经验平方和${\mathbf{y}}^{\top }\mathbf{y}$，再加上由于$\mathbf{w}$上先验误差所得到的项。

后验边缘如下所示。对于方差，我们有：

$$p({\sigma }^2|\mathcal{D})=\int p(\mathbf{w}|{\sigma }^2,\mathcal{D})p({\sigma }^2|\mathcal{D})\mathrm{d}\mathbf{w}=\text{IG}({\sigma }^2|\hat{a},\hat{b})\text{\hspace{1em}(24)}$$

对于回归权重，可以证明：

$$p(\mathbf{w}|\mathcal{D})=\int p(\mathbf{w}|{\sigma }^2,\mathcal{D})p({\sigma }^2|\mathcal{D})\mathrm{d}{\sigma }^2=\mathcal{T}(\mathbf{w}|\hat{\mathbf{w}},\frac{\hat{b}}{\hat{a}}\hat{\mathbf{\Sigma }},2\hat{a})\text{\hspace{1em}(25)}$$

15.2.2.3 后验预测分布

在机器学习中，我们通常更关心预测的不确定性（以及准确率），而不是参数估计。幸运的是，我们可以推导出封闭形式的后验预测分布。具体来说，可以证明，给定$N^{\prime }$个新的测试输入$\tilde{\mathbf{X}}$，我们有：

$$p(\mathbf{y}|\mathbf{X},\mathcal{D})=\iint p(\mathbf{y}|\mathbf{X},\mathbf{w},{\sigma }^2)p(\mathbf{w},{\sigma }^2|\mathcal{D})\mathrm{d}\mathbf{w}\mathrm{d}{\sigma }^2\text{\hspace{1em}(26)}$$

$$=\iint \mathcal{N}(\mathbf{y}|\mathbf{X}\mathbf{w},{\sigma }^2{\mathbf{I}}_{N^{\prime }})\text{NIG}(\mathbf{w},{\sigma }^2|\hat{\mathbf{w}},\hat{\mathbf{\Sigma }},\hat{a},\hat{b})\mathrm{d}\mathbf{w}\mathrm{d}{\sigma }^2\text{\hspace{1em}(27)}$$

$$=\mathcal{T}(\mathbf{y}|\mathbf{X}\hat{\mathbf{w}},\frac{\hat{b}}{\hat{a}}({\mathbf{I}}_{N^{\prime }}+\mathbf{X}\hat{\mathbf{\Sigma }}{\mathbf{X}}^{\top }),2\hat{a})\text{\hspace{1em}(28)}$$

后验预测均值相当于"正态"线性回归，但其中我们插入$\hat{\mathbf{w}}=\mathbb{E}[\mathbf{w}|\mathcal{D}]$而不是最大似然估计。后验预测方差有两个分量：由于测量噪声导致的$\frac{\hat{b}}{\hat{a}}{\mathbf{I}}_{N^{\prime }}$，以及由于$\mathbf{w}$的不确定性导致的$\frac{\hat{b}}{\hat{a}}\mathbf{X}\hat{\mathbf{\Sigma }}{\mathbf{X}}^{\top }$。后一项根据测试输入与训练数据的接近程度而变化。结果类似于使用高斯先验（具有固定的${\hat{\sigma }}^2$）的情况，只是预测分布更宽，因为我们考虑了${\sigma }^2$的不确定性。

无信息性先验

对贝叶斯推理的一个常见非议是贝叶斯推理需要使用先验知识。这种现象有时被认为"污染"了人们从数据中所得出的推论。我们可以通过使用无信息性先验来最小化先验的影响。下面将讨论线性回归的各种无信息性先验。

Jeffreys先验

我们知道对于位置参数的Jeffreys先验，其形式为$p(\mathbf{w})\propto 1$。我们知道对于尺度因子的Jeffreys先验，其形式是$p(\sigma )\propto {\sigma }^{-1}$。我们可以使用不适当的正态逆伽马分布先验来模拟这些先验，其中$\mathbf{w}=0$，$\mathbf{\Sigma }=\infty \mathbf{I}$，$a=-D/2$，$b=0$。对应的后验由下式给出：

$$\hat{\mathbf{w}}={\hat{\mathbf{w}}}_{\text{MLE}}=({\mathbf{X}}^{\top }\mathbf{X}{)}^{-1}{\mathbf{X}}^{\top }\mathbf{y}\text{\hspace{1em}(29)}$$

$$\hat{\mathbf{\Sigma }}=({\mathbf{X}}^{\top }\mathbf{X}{)}^{-1}\triangleq \mathbf{C}\text{\hspace{1em}(30)}$$

$$\hat{a}=\frac{\nu }{2}\text{\hspace{1em}(31)}$$

$$\hat{b}=\frac{s^2\nu }{2}\text{\hspace{1em}(32)}$$

$$s^2\triangleq \frac{\Vert \mathbf{y}-\hat{\mathbf{y}}{\Vert }^2}{\nu }\text{\hspace{1em}(33)}$$

$$\nu =N-D\text{\hspace{1em}(34)}$$

因此，权重的后验分布由下式给出：

$$p(\mathbf{w}|\mathcal{D})=\mathcal{T}(\mathbf{w}|\hat{\mathbf{w}},s^2\mathbf{C},\nu )\text{\hspace{1em}(35)}$$

其中，$\hat{\mathbf{w}}$是最大似然估计。因此，每个权重的边缘概率具有以下形式：

$$p(w_d|\mathcal{D})=\mathcal{T}(w_d|{\hat{w}}_d,s^2{C}_{dd},\nu )\text{\hspace{1em}(36)}$$

与频率学派统计学的联系

有趣的是，当使用Jeffreys先验时，后验在形式上等价于最大似然估计的频率采样分布，其形式为：

$$p({\hat{w}}_d|{\mathcal{D}}^{\ast })=\mathcal{T}({\hat{w}}_d|w_d,s^2{C}_{dd},\nu )\text{\hspace{1em}(37)}$$

其中，${\mathcal{D}}^{\ast }=(\mathbf{X},{\mathbf{y}}^{\ast })$是在给定（固定）输入$\mathbf{X}$的情况下，从真实模型生成的假设数据。在关于频率学派统计学的书籍中，这通常被写成以下等价形式：

$$\frac{{\hat{w}}_d-w_d}{s\sqrt{C_{dd}}}\sim t_{N-D}\text{\hspace{1em}(38)}$$

采样分布在数值上与式(39)中的后验分布相同，因为$\mathcal{T}(w|\mu ,{\sigma }^2,\nu )=\mathcal{T}(\mu |w,{\sigma }^2,\nu )$。然而，采样分布在语义上大不相同，因为采样分布不以观察到的数据为条件，而是基于从模型中提取的假设数据。更多有关使用无信息性先验时简单线性模型的贝叶斯分析和频率学派统计学分析之间等价性的讨论。

Zellner的g-先验

假设${\sigma }^2$上有一个无信息性先验，这通常是合理的假设，因为${\sigma }^2$只是一个对结果没有太大影响的标量，但对$\mathbf{w}$使用无信息性先验可能是危险的，正如我们从岭回归中所知的，先验的强度控制着模型的正则化程度。

一个常见的折衷方案是使用正态逆伽马分布先验，$a=-D/2$，$b=0$（以确保$p({\sigma }^2)\propto 1$），$\mathbf{w}=0$，$\mathbf{\Sigma }=g({\mathbf{X}}^{\top }\mathbf{X}{)}^{-1}$，其中$g>0$的作用类似于岭回归中的$1/\lambda$。这被称为Zellner的g-先验。我们发现，先验协方差与$({\mathbf{X}}^{\top }\mathbf{X}{)}^{-1}$成正比，而不是与$\mathbf{I}$成正比；这确保了后验对输入的缩放是不变的。

有了这个先验，后验就变成了：

$$p(\mathbf{w},{\sigma }^2|g,\mathcal{D})=\text{NIG}(\mathbf{w},{\sigma }^2|{\mathbf{w}}_N,{\mathbf{V}}_N,a_N,b_N)\text{\hspace{1em}(40)}$$

注意，这个先验是以输入$\mathbf{X}$为条件的，而不是以输出$\mathbf{y}$为条件的；这在条件中是完全有效的，其中所有计算都以$\mathbf{X}$为条件，它被视为固定常数输入。

$${\mathbf{V}}_N=\frac{g}{g+1}({\mathbf{X}}^{\top }\mathbf{X}{)}^{-1}\text{\hspace{1em}(41)}$$

$${\mathbf{w}}_N=\frac{g}{g+1}{\hat{\mathbf{w}}}_{\text{mle}}\text{\hspace{1em}(42)}$$

$$a_N=N/2\text{\hspace{1em}(43)}$$

$$b_N=\frac{s^2}{2}+\frac{1}{2(g+1)}{\hat{\mathbf{w}}}_{\text{mle}}^{\top }{\mathbf{X}}^{\top }\mathbf{X}{\hat{\mathbf{w}}}_{\text{mle}}\text{\hspace{1em}(44)}$$

已经提出了各种方法来设置$g$，这些方法包括交叉验证、经验贝叶斯、层次贝叶斯等。

“尖峰和平板”先验

在执行预测时，通常需要能够选择输入特征的子集，以减少过拟合，或者提高模型的可解释性。如果我们确保权重向量$\mathbf{w}$是稀疏的（即，具有许多零元素），就可以实现这一点，因为如果$w_d=0$，则$x_d$在内积${\mathbf{w}}^{\top }\mathbf{x}$中就不起作用。

当使用贝叶斯推理时，实现稀疏性的规范方法是使用"尖峰和平板"(spike-and-slab)先验，该方法具有两个分量的混合模型形式，其中一个分量是在0处的"尖峰"，另一个分量是$-a$和$a$之间的均匀"平板"：

$$p(\mathbf{w})=\prod _{d=1}^D(1-\pi )\delta (w_d)+\pi \text{Unif}(w_d|-a,a)\text{\hspace{1em}(45)}$$

其中，$\pi$是每个系数为非零的先验。因此，相应的系数对数先验为：

$$\log p(\mathbf{w})=\Vert \mathbf{w}{\Vert }_0\log (1-\pi )+(D-\Vert \mathbf{w}{\Vert }_0)\log \pi =-\lambda \Vert \mathbf{w}{\Vert }_0+\text{常数}\text{\hspace{1em}(46)}$$

其中，$\lambda =\log \frac{\pi }{1-\pi }$控制模型的稀疏性，$\Vert \mathbf{w}{\Vert }_0={\sum }_{d=1}^D\mathbb{I}(w_d\ne 0)$是权重的${\ell }_0$范数(norm)。

因此，具有"尖峰和平板"先验的最大后验估计等价于${\ell }_0$正则化(regularization)；这惩罚了非零系数的数量。有趣的是，后验样本也是稀疏的。

相比之下，考虑使用拉普拉斯先验。有趣的是，使用具有"尖峰和平板"先验的后验均值比使用具有拉普拉斯先验的后验众数具有更好的预测准确率。

在实践中，我们经常使用宽的高斯分布来近似均匀的"平板"，

$$p(\mathbf{w})=\prod _d(1-\pi )\delta (w_d)+\pi \mathcal{N}(w_d|0,{\sigma }_w^2)\text{\hspace{1em}(47)}$$

当${\sigma }_w^2\to \infty$时，第二项在$[-1,+1]$上接近均匀分布。我们可以通过将二元随机变量与每个系数相关联来实现混合模型$s_d\sim \text{Ber}(\pi )$，以指示该系数是"打开"还是"关闭"状态。

遗憾的是，具有这种离散混合先验的最大后验估计（更不用说完全贝叶斯推理）在计算上非常困难。研究学者已经提出了各种近似推理方法，包括贪婪搜索或马尔可夫链蒙特卡罗方法。

拉普拉斯先验（贝叶斯套索）

为了实现稀疏性，一种具有较低计算成本的方法是，通过最小化带惩罚项的负对数似然并利用拉普拉斯先验执行最大后验估计：

$$\text{PNLL}(\mathbf{w})=-\log p(\mathcal{D}|\mathbf{w})-\log p(\mathbf{w}|\lambda )=\Vert \mathbf{Xw}-\mathbf{y}{\Vert }_2^2+\lambda \Vert \mathbf{w}{\Vert }_1\text{\hspace{1em}(48)}$$

其中，$\Vert \mathbf{w}{\Vert }_1\triangleq {\sum }_{d=1}^D|w_d|$是$\mathbf{w}$的${\ell }_1$范数。这种方法被称为套索，代表"最小绝对收缩和选择算子"(least absolute shrinkage and selection operator)。

我们将讨论具有该先验的后验推理；这被称为贝叶斯套索(Bayesian lasso)。具体来说，我们假设以下先验分布：

$$p(\mathbf{w}|{\sigma }^2)=\prod _j\frac{\lambda }{2\sqrt{{\sigma }^2}}{e}^{-\lambda |w_j|/\sqrt{{\sigma }^2}}\text{\hspace{1em}(49)}$$

（注意，为了确保完全后验是单模态的，将先验条件设置为依赖于${\sigma }^2$是很重要的。）

为了简化推理，将拉普拉斯先验表示为高斯缩放的混合。特别地，可以证明拉普拉斯分布是高斯的无限加权和，其中精度来自伽马分布：

$$\text{Laplace}(w|0,\lambda )=\int \mathcal{N}(w|0,{\tau }^2)\text{Ga}\left({\tau }^2|1,\frac{{\lambda }^2}{2}\right)\mathrm{d}{\tau }^2\text{\hspace{1em}(50)}$$

因此，我们可以将贝叶斯套索模型表示为一个层次潜在变量模型，如图15-3a所示。对应的联合分布具有以下形式：

$$p(\mathbf{y},\mathbf{w},\mathbf{\tau },{\sigma }^2|\mathbf{X})=\mathcal{N}(\mathbf{y}|\mathbf{Xw},{\sigma }^2{\mathbf{I}}_N)\left[\prod _j\mathcal{N}(w_j|0,{\sigma }^2{\tau }_j^2)\text{Ga}({\tau }_j^2|1,{\lambda }^2/2)\right]p({\sigma }^2)\text{\hspace{1em}(51)}$$

我们还可以创建一个高斯缩放的混合模型来匹配组套索先验，该模型同时将多个系数设置为零：

$${\mathbf{w}}_g|{\sigma }^2,{\tau }_g^2\sim \mathcal{N}(0,{\sigma }^2{\tau }_g^2{\mathbf{I}}_{d_g})\text{\hspace{1em}(52)}$$

$${\tau }_g^2\sim \text{Ga}\left(\frac{d_g+1}{2},\frac{{\lambda }^2}{2}\right)\text{\hspace{1em}(53)}$$

其中，$d_g$是组$g$的大小。因此，我们看到每个组有一个方差项，每个方差项都来自一个伽马先验，其形状参数取决于组的大小，其速率参数由$\gamma$控制。

给定这些层次模型，我们可以很容易地推导出一个吉布斯采样算法来从后验分布进行采样。遗憾的是，即使最大后验估计是稀疏的，这些后验样本也不是稀疏的。这是因为先验对每个系数为零的事件赋予了无穷小的概率。

“马蹄”先验

拉普拉斯先验不适用于稀疏贝叶斯模型，因为后验样本不是稀疏的。"尖峰和平板"先验则没有这个问题，但这个方法通常速度太慢。幸运的是，可以设计出既稀疏又计算高效的连续先验（没有离散的潜在变量）。这种类型的一个流行先验是"马蹄"先验(horseshoe prior)，它因为其密度函数的形状而得名。

在"马蹄"先验中，我们使用以下高斯缩放混合，而不是对每个权重使用拉普拉斯先验：

$$w_j\sim \mathcal{N}(0,{\lambda }_j^2{\tau }^2)\text{\hspace{1em}(54)}$$

$${\lambda }_j\sim {\mathcal{C}}_{+}(0,1)\text{\hspace{1em}(55)}$$

$${\tau }^2\sim {\mathcal{C}}_{+}(0,1)\text{\hspace{1em}(56)}$$

其中，${\mathcal{C}}_{+}(0,1)$是半柯西分布（见2.2.2.4节），${\lambda }_j$是局部收缩因子，${\tau }^2$是全局收缩因子。柯西分布具有非常粗的尾部，因此${\lambda }_j$可能为0或离0很远，这模仿了之前的"尖峰和平板"先验，但采取连续的方式。

自动相关性确定

使用促进稀疏性的先验进行后验推理的另一种选择是，使用高斯先验$w_j\sim \mathcal{N}(0,1/{\alpha }_j)$进行后验推理，但其中我们使用经验贝叶斯来优化精度${\alpha }_j$。也就是说，我们首先计算$\hat{\mathbf{\alpha }}=\arg {\max }_{\mathbf{\alpha }}p(\mathbf{y}|\mathbf{X},\mathbf{\alpha })$，然后计算$\hat{\mathbf{w}}=\arg {\max }_{\mathbf{w}}\mathcal{N}(\mathbf{w}|0,{\hat{\mathbf{A}}}^{-1})$。也许令人惊讶的是，由于在15.2.8.2节中解释的原因，我们发现结果会导致稀疏的估计。

这种技术被称为稀疏贝叶斯学习(sparse Bayesian learning)或自动相关性确定(Automatic Relevance Determination, ARD)。该技术也被称为NUV估计，代表"具有未知方差的正态先验"(normal prior with unknown variance)。这个技术最初是为神经网络开发的（其中，稀疏性应用于第一层权重），但此处我们将其应用于线性模型。

线性模型的自动相关性确定

我们通过将自动相关性确定应用于线性回归来更详细地解释这个概念。似然为$p(\mathbf{y}|\mathbf{x},\mathbf{w},\beta )=\mathcal{N}(\mathbf{y}|{\mathbf{w}}^{\top }\mathbf{x},1/\beta )$，其中$\beta =1/{\sigma }^2$。先验是$p(\mathbf{w})=\mathcal{N}(\mathbf{w}|0,{\mathbf{A}}^{-1})$，其中$\mathbf{A}=\text{diag}(\mathbf{\alpha })$。边缘似然可以解析计算，如下所示：

$$p(\mathbf{y}|\mathbf{X},\mathbf{\alpha },\beta )=\int \mathcal{N}(\mathbf{y}|\mathbf{Xw},(1/\beta ){\mathbf{I}}_N)\mathcal{N}(\mathbf{w}|0,{\mathbf{A}}^{-1})\mathrm{d}\mathbf{w}\text{\hspace{1em}(57)}$$

$$=\mathcal{N}(\mathbf{y}|0,{\beta }^{-1}{\mathbf{I}}_N+\mathbf{X}{\mathbf{A}}^{-1}{\mathbf{X}}^{\top })\text{\hspace{1em}(58)}$$

$$=\mathcal{N}(\mathbf{y}|0,{\mathbf{C}}_{\mathbf{\alpha }})\text{\hspace{1em}(59)}$$

其中，${\mathbf{C}}_{\mathbf{\alpha }}\triangleq {\beta }^{-1}{\mathbf{I}}_N+\mathbf{X}{\mathbf{A}}^{-1}{\mathbf{X}}^{\top }$。这与"尖峰和平板"先验下的边缘似然非常相似，由下式给出：

$$p(\mathbf{y}|\mathbf{X},\mathbf{s},{\sigma }_y^2,{\sigma }_w^2)=\int \mathcal{N}(\mathbf{y}|{\mathbf{X}}_s{\mathbf{w}}_s,{\sigma }_y^2\mathbf{I})\mathcal{N}({\mathbf{w}}_s|0_s,{\sigma }_w^2\mathbf{I})\mathrm{d}{\mathbf{w}}_s=\mathcal{N}(\mathbf{y}|0,{\mathbf{C}}_s)\text{\hspace{1em}(60)}$$

其中，${\mathbf{C}}_s={\sigma }_y^2{\mathbf{I}}_N+{\sigma }_w^2{\mathbf{X}}_s{\mathbf{X}}_s^{\top }$。（此处，${\mathbf{X}}_s$指的是设计矩阵，而我们只选择$\mathbf{X}$中满足$s_d=1$的那些列。）不同之处在于，我们将二元变量$s_j\in \{0,1\}$替换为连续变量${\alpha }_j\in {\mathbb{R}}^{+}$，这使得优化问题变得更加简单。

目标是对数边缘似然，由下式给出：

$$\ell (\mathbf{\alpha },\beta )=-\frac{1}{2}\log p(\mathbf{y}|\mathbf{X},\mathbf{\alpha },\beta )=\log |{\mathbf{C}}_{\mathbf{\alpha }}|+{\mathbf{y}}^{\top }{\mathbf{C}}_{\mathbf{\alpha }}^{-1}\mathbf{y}\text{\hspace{1em}(61)}$$

可以使用各种优化算法来优化$\ell (\mathbf{\alpha },\beta )$。

自动相关性确定可以用作${\ell }_1$正则化的替代方案。尽管自动相关性确定目标不是凸函数，但它往往会给出更稀疏的结果。此外，可以证明，自动相关性确定目标的局部最优值比${\ell }_0$正则化目标少得多，因此更容易优化。

为什么自动相关性确定会导致稀疏解

一旦我们估计了$\mathbf{\alpha }$和$\beta$，就可以使用高斯的贝叶斯规则来计算参数的后验，得到$p(\mathbf{w}|\mathcal{D},\hat{\mathbf{\alpha }},\hat{\beta })=\mathcal{N}(\mathbf{w}|\hat{\mathbf{w}},\hat{\mathbf{\Sigma }})$，其中${\hat{\mathbf{\Sigma }}}^{-1}=\hat{\beta }{\mathbf{X}}^{\top }\mathbf{X}+\mathbf{A}$和$\hat{\mathbf{w}}=\hat{\beta }\hat{\mathbf{\Sigma }}{\mathbf{X}}^{\top }\mathbf{y}$。如果我们有${\hat{\alpha }}_d\to \infty$，那么${\hat{w}}_d\to 0$，因此解向量将是稀疏的。

现在，我们说明在什么情况下这种稀疏解可能是最优的。为了简单起见，我们假设$\beta =1/{\sigma }^2$是固定的。考虑具有两个训练样例的一维线性回归，因此$\mathbf{X}=\mathbf{x}=(x_1,x_2)$，$\mathbf{y}=(y_1,y_2)$。我们可以将$\mathbf{x}$和$\mathbf{y}$绘制为平面中的向量。假设该特征与预测响应无关，因此$\mathbf{x}$指向与$\mathbf{y}$几乎正交的方向。让我们看看当改变$\alpha$时，边缘似然会发生什么变化。边缘似然由$p(\mathbf{y}|\mathbf{x},\alpha ,\beta )=\mathcal{N}(\mathbf{y}|0,{\mathbf{C}}_{\alpha })$给出，其中${\mathbf{C}}_{\alpha }=\frac{1}{\beta }\mathbf{I}+\frac{1}{\alpha }\mathbf{x}{\mathbf{x}}^{\top }$。如果$\alpha$是有限的，后验分布将沿着$\mathbf{x}$方向拉长。然而，如果$\alpha =\infty$，我们有${\mathbf{C}}_{\alpha }=\frac{1}{\beta }\mathbf{I}$，这是球形的，如果$|{\mathbf{C}}_{\alpha }|$保持不变，后者将更高的概率密度分配给观察到的响应向量$\mathbf{y}$，因此这是优选的解。换言之，当${\alpha }_d$较小但${\mathbf{X}}_{:,d}$不相关时，边缘似然会"惩罚"这些解，因为它们浪费了概率质量。（从贝叶斯奥卡姆剃刀的角度来看）消除冗余维度是更简约的。

理解自动相关性确定稀疏性的另一种方法是在层次贝叶斯模型中进行近似推理。特别是，假设我们将共轭先验放置在每个精度${\alpha }_d\sim \text{Ga}(a,b)$，和观测精度$\beta \sim \text{Ga}(c,d)$上。由于学生先验的精确推理是难以处理的，我们可以使用变分贝叶斯，使用以下形式的因子后验分布：

$$q(\mathbf{w},\mathbf{\alpha })=q(\mathbf{w})q(\mathbf{\alpha })\approx \mathcal{N}(\mathbf{w}|\mathbf{\mu },\mathbf{\Sigma })\prod _d\text{Ga}({\alpha }_d|{\hat{a}}_d,{\hat{b}}_d)\text{\hspace{1em}(62)}$$

自动相关性确定通过点估计来近似$q(\mathbf{\alpha })$。然而，在变分贝叶斯中，我们对$\mathbf{\alpha }$进行积分，由此得到的权重上的后验边缘$q(\mathbf{w})$由下式给出：

$$p(\mathbf{w}|\mathcal{D})=\int \mathcal{N}(\mathbf{w}|0,\text{diag}(\mathbf{\alpha }{)}^{-1})\prod _d\text{Ga}({\alpha }_d|{\hat{a}}_d,\hat{b})\mathrm{d}\mathbf{\alpha }\text{\hspace{1em}(63)}$$

这是一个高斯缩放的混合模型，可以证明，该混合模型与多变量学生分布相同，具有非对角协方差。需要注意的是，学生分布在0处有一个大尖峰，这直观地解释了为什么后验均值（对于学生分布，等价于后验众数）是稀疏的。

最后，我们还可以将自动相关性确定视为具有非因子先验(non-factorial prior)的最大后验估计问题。直观地说，尽管使用了对角高斯先验分布，但$w_j$参数之间的依赖性仍然会出现，因为先验精度${\alpha }_j$是在边缘化所有$\mathbf{w}$之后估计的，因此依赖于所有特征。有趣的是，我们证明了具有非因子先验的最大后验估计在以下意义上严格优于具有任何可能的因子先验的最大后验估计：非因子目标总是比因子目标具有更少的局部极小值，同时仍然满足非因子目标的全局最优与${\ell }_0$目标的全局优化相对应的性质——${\ell }_0$正则化没有局部极小值，因此不具有这种性质。

自动相关性确定算法

可以使用多种算法来优化$\ell (\mathbf{\alpha },\beta )$。一种方法是使用期望最大化，其中我们在期望步骤中计算$p(\mathbf{w}|\mathcal{D},\mathbf{\alpha },\beta )$，然后在最大化步骤中最大化$\mathbf{\alpha }$。在变分贝叶斯中，我们同时推理$\mathbf{w}$和$\mathbf{\alpha }$。

最近，我们发现，这些方法执行的嵌套迭代计算可以通过循环神经网络来进行模拟（见16.3.4节）。此外，通过训练该模型，可以实现比手动设计的优化算法更快的收敛。

相关向量机

假设我们创建一个形式为$p(y|\mathbf{x};\mathbf{\theta })=\mathcal{N}(y|{\mathbf{w}}^{\top }\mathbf{\varphi }(\mathbf{x}),{\sigma }^2)$的线性回归模型，其中$\mathbf{\varphi }(\mathbf{x})=[\mathcal{K}(\mathbf{x},{\mathbf{x}}_1),\dots ,\mathcal{K}(\mathbf{x},{\mathbf{x}}_N)]$，$\mathcal{K}()$是核函数，并且${\mathbf{x}}_1,\dots ,{\mathbf{x}}_N$是$N$个训练点。这被称为核基函数展开(kernel basis function expansion)，它将输入从$\mathbf{x}\in \mathcal{X}$转换为$\mathbf{\varphi }(\mathbf{x})\in {\mathbb{R}}^N$。显然，这个模型具有$O(N)$个参数，因此是非参数化模型。然而，我们可以使用自动相关性确定来选择样本的一小部分。这种技术被称为相关向量机(Relevance Vector Machine, RVM)。

多元线性回归

我们将考虑多元线性回归模型，其形式为：

$$\mathbf{Y}=\mathbf{WX}+\mathbf{E}\text{\hspace{1em}(64)}$$

其中，$\mathbf{W}\in {\mathbb{R}}^{N_y\times N_x}$是回归系数的矩阵，$\mathbf{X}\in {\mathbb{R}}^{N_x\times N}$是输入特征的矩阵（每行对应一个输入变量，每列对应一种观察值），$\mathbf{Y}\in {\mathbb{R}}^{N_y\times N}$是响应矩阵（每行对应一个输出变量，每列对应一种观察值），$\mathbf{E}=[{\mathbf{e}}_1,\dots ,{\mathbf{e}}_N]$是残差矩阵，其中${\mathbf{e}}_i\stackrel{\text{iid}}{\sim }\mathcal{N}(0,\mathbf{\Sigma })$。从定义中可以看出，如果给定$\mathbf{\Sigma }$、$\mathbf{W}$以及$\mathbf{X}$，则$\mathbf{Y}$的列是遵循多元正态分布的独立随机变量。所以观察值的似然为：

$$p(\mathbf{Y}|\mathbf{W},\mathbf{X},\mathbf{\Sigma })=\frac{1}{(2\pi {)}^{N_{y}N/2}|\mathbf{\Sigma }{|}^{N/2}}\exp \left(\sum _{i=1}^N-\frac{1}{2}({\mathbf{y}}_i-\mathbf{W}{\mathbf{x}}_i{)}^{\top }{\mathbf{\Sigma }}^{-1}({\mathbf{y}}_i-\mathbf{W}{\mathbf{x}}_i)\right)\text{\hspace{1em}(65)}$$

$$=\frac{1}{(2\pi {)}^{N_{y}N/2}|\mathbf{\Sigma }{|}^{N/2}}\exp \left(-\frac{1}{2}\text{tr}((\mathbf{Y}-\mathbf{WX}{)}^{\top }{\mathbf{\Sigma }}^{-1}(\mathbf{Y}-\mathbf{WX}))\right)\text{\hspace{1em}(66)}$$

$$=\mathcal{M}\mathcal{N}(\mathbf{Y}|\mathbf{WX},\mathbf{\Sigma },{\mathbf{I}}_{N\times N})\text{\hspace{1em}(67)}$$

其共轭先验是矩阵正态逆威沙特分布：

$$\mathbf{W},\mathbf{\Sigma }\sim \text{MNIW}({\mathbf{M}}_0,{\mathbf{V}}_0,{\nu }_0,{\mathbf{\Psi }}_0)\text{\hspace{1em}(68)}$$

其中，矩阵正态逆威沙特分布被定义为：

$$\mathbf{W}|\mathbf{\Sigma }\sim \mathcal{M}\mathcal{N}({\mathbf{M}}_0,\mathbf{\Sigma },{\mathbf{V}}_0)\text{\hspace{1em}(69)}$$

$$\mathbf{\Sigma }\sim \text{IW}({\nu }_0,{\mathbf{\Psi }}_0)\text{\hspace{1em}(70)}$$

其中，${\mathbf{V}}_0\in {\mathbb{R}}_{++}^{N_x\times N_x}$，${\mathbf{\Psi }}_0\in {\mathbb{R}}_{++}^{N_y\times N_y}$，$\nu_0>N_x-1$是逆威沙特分布的自由度。

$\{\mathbf{W},\mathbf{\Sigma }\}$的后验分布仍然遵循矩阵正态逆威沙特分布。联合分布的密度为：

$$p(\mathbf{Y},\mathbf{W},\mathbf{\Sigma })\propto |\mathbf{\Sigma }{|}^{-({\nu }_0+N_y+1+N_x+N)/2}\times \exp \left\{-\frac{1}{2}\text{tr}({\mathbf{\Omega }}_0)\right\}\text{\hspace{1em}(71)}$$

$${\mathbf{\Omega }}_0\triangleq {\mathbf{\Psi }}_0{\mathbf{\Sigma }}^{-1}+(\mathbf{Y}-\mathbf{WX}{)}^{\top }{\mathbf{\Sigma }}^{-1}(\mathbf{Y}-\mathbf{WX})+(\mathbf{W}-{\mathbf{M}}_0{)}^{\top }{\mathbf{\Sigma }}^{-1}(\mathbf{W}-{\mathbf{M}}_0){\mathbf{V}}_0\text{\hspace{1em}(72)}$$

我们首先聚合指数中包括$\mathbf{W}$的项，使其采用矩阵正态分布的"完成平方"技术。这类似于我们推导多元正态分布的共轭后验时使用的技术。具体来说，即：

$$\text{tr}[(\mathbf{Y}-\mathbf{WX}{)}^{\top }{\mathbf{\Sigma }}^{-1}(\mathbf{Y}-\mathbf{WX})+(\mathbf{W}-{\mathbf{M}}_0{)}^{\top }{\mathbf{\Sigma }}^{-1}(\mathbf{W}-{\mathbf{M}}_0){\mathbf{V}}_0]\text{\hspace{1em}(73)}$$

$$=\text{tr}({\mathbf{\Sigma }}^{-1}[(\mathbf{Y}-\mathbf{WX})(\mathbf{Y}-\mathbf{WX}{)}^{\top }+(\mathbf{W}-{\mathbf{M}}_0){\mathbf{V}}_0(\mathbf{W}-{\mathbf{M}}_0{)}^{\top }])\text{\hspace{1em}(74)}$$

$$=\text{tr}({\mathbf{\Sigma }}^{-1}[\mathbf{W}{\mathbf{S}}_{xx}{\mathbf{W}}^{\top }-2{\mathbf{S}}_{yx}{\mathbf{W}}^{\top }+{\mathbf{S}}_{yy}])\text{\hspace{1em}(75)}$$

$$=\text{tr}({\mathbf{\Sigma }}^{-1}[(\mathbf{W}-{\mathbf{S}}_{yx}{\mathbf{S}}_{xx}^{-1}){\mathbf{S}}_{xx}(\mathbf{W}-{\mathbf{S}}_{yx}{\mathbf{S}}_{xx}^{-1}{)}^{\top }+{\mathbf{S}}_{y|x}])\text{\hspace{1em}(76)}$$

其中，

$${\mathbf{S}}_{xx}=\mathbf{X}{\mathbf{X}}^{\top }+{\mathbf{V}}_0,{\mathbf{S}}_{yx}=\mathbf{Y}{\mathbf{X}}^{\top }+{\mathbf{M}}_0{\mathbf{V}}_0\text{\hspace{1em}(77)}$$

$${\mathbf{S}}_{yy}=\mathbf{Y}{\mathbf{Y}}^{\top }+{\mathbf{M}}_0{\mathbf{V}}_0{\mathbf{M}}_0^{\top },{\mathbf{S}}_{y|x}={\mathbf{S}}_{yy}-{\mathbf{S}}_{yx}{\mathbf{S}}_{xx}^{-1}{\mathbf{S}}_{yx}^{\top }\text{\hspace{1em}(78)}$$

因此，从式(76)可以看出，如果给定$\mathbf{\Sigma }$，则$\mathbf{W}$遵循矩阵正态分布：

$$\mathbf{W}|\mathbf{\Sigma },\mathbf{X},\mathbf{Y}\sim \mathcal{M}\mathcal{N}({\mathbf{S}}_{yx}{\mathbf{S}}_{xx}^{-1},\mathbf{\Sigma },{\mathbf{S}}_{xx})\text{\hspace{1em}(79)}$$

通过将$\mathbf{W}$边缘化（即去除式(72)中包含$\mathbf{W}$的指数项），可以证明$\mathbf{\Sigma }$的后验分布是逆威沙特分布。实际上，通过将式(76)替换为式(72)中的相应项，可以看出积分出$\mathbf{W}$后只剩下${\mathbf{\Sigma }}^{-1}\mathbf{\Psi }$和${\mathbf{\Sigma }}^{-1}{\mathbf{S}}_{y|x}$，这表明后验逆威沙特分布的标度矩阵是${\mathbf{\Psi }}_0+{\mathbf{S}}_{y|x}$。

总之，对于给定的观测值，$\{\mathbf{W},\mathbf{\Sigma }\}$的联合后验分布为：

$$\mathbf{W},\mathbf{\Sigma }|\mathbf{X},\mathbf{Y}\sim \text{MNIW}({\mathbf{M}}_1,{\mathbf{V}}_1,{\nu }_1,{\mathbf{\Psi }}_1)\text{\hspace{1em}(80)}$$

$${\mathbf{M}}_1={\mathbf{S}}_{yx}{\mathbf{S}}_{xx}^{-1}\text{\hspace{1em}(81)}$$

$${\mathbf{V}}_1={\mathbf{S}}_{xx}\text{\hspace{1em}(82)}$$

$${\nu }_1=N+{\nu }_0\text{\hspace{1em}(83)}$$

$${\mathbf{\Psi }}_1={\mathbf{\Psi }}_0+{\mathbf{S}}_{y|x}\text{\hspace{1em}(84)}$$

$\mathbf{W}$和$\mathbf{\Sigma }$的最大后验估计是后验矩阵正态逆威沙特分布的模式。为了推导这一点，请注意，$\mathbf{W}$只出现在矩阵正态密度函数的后验分布中，因此最大化$\{\mathbf{W},\mathbf{\Sigma }\}$后验密度的矩阵$\mathbf{W}$，是使$\mathbf{W}$的矩阵正态后验最大化的矩阵。因此，$\mathbf{W}$的最大后验估计为$\hat{\mathbf{W}}={\mathbf{M}}_1={\mathbf{S}}_{yx}{\mathbf{S}}_{xx}^{-1}$，这适用于$\mathbf{\Sigma }$的任何值。通过将$\mathbf{W}=\hat{\mathbf{W}}$插入$\{\mathbf{W},\mathbf{\Sigma }\}$的联合后验分布，并取$\mathbf{\Sigma }$的导数，可以发现密度最大化的矩阵是$({\nu }_1+N_y+N_x+1{)}^{-1}{\mathbf{\Psi }}_1$。由于${\mathbf{\Psi }}_1$是正定的，因此${\mathbf{\Psi }}_1$是$\mathbf{\Sigma }$的最大后验估计。

总之，$\{\mathbf{W},\mathbf{\Sigma }\}$的最大后验估计为：

$$\hat{\mathbf{W}}={\mathbf{S}}_{yx}{\mathbf{S}}_{xx}^{-1}\text{\hspace{1em}(85)}$$

$$\hat{\mathbf{\Sigma }}=\frac{1}{{\nu }_1+N_y+N_x+1}({\mathbf{\Psi }}_0+{\mathbf{S}}_{y|x})\text{\hspace{1em}(86)}$$