Book Notes on Pattern Recognition and Machine Learning (PRML)
我们假设输入为 n 个数据，每个数据为 p 维；输出为 n 个一维数据。\(X_i\), \(\hat{y}\), \(\hat{\beta}\) 均为列向量，\(X=[X_1^T,X_2^T,...,X_n^T]^T\) 那么此时线性模型有如下形式
\[f(X) = \beta_0 + \sum_{i=1}^pX_i\beta_i\]到后来我们可以发现，这里的截距 \(\beta_0\) 是用于补偿训练输入的标签均值和训练数据的加权平均之间的差异。为了方便表示，我们将 \(\beta_0\) 整合进模型，即此时 \(X=[1, X_1^T,X_2^T,...,X_n^T]^T\) 为 n*(p+1) 的矩阵，预测输出 \(\hat{y}\) 为 n*1 的列向量，线性模型 \(\hat{\beta}\) 为 (p+1)*1 的列向量，则此时我们有
\[Y^* = X\beta + \epsilon\]这里 \(y^*\) 表示ground-truth；\(\epsilon\) 为误差项，是 n*1 的列向量，假设服从正态分布；\(Y_{pred} = X\beta\)
值得注意的是，当我们需要同时输出多个值时（即 \(\boldsymbol{y_i}\) 为 K 维 vector），此时我们完全可以将其看成 K 个独立的回归问题。这是因为在高斯噪声分布下，我们的模型参数 \(\beta\) 只定义了噪声分布的均值，而它与方差独立（在极大似然的解下）。
在回归任务中，我们希望误差越小越好，也就是说，我们希望 \(min_{\beta}\ \epsilon\) 。在这其中，最常用的一种估计方法是最小二乘，用来最小化残差和，这里我们首先给出它的形式
\(RSS(\beta) = \sum_{i=1}^n (y_i^*-x_i^T\beta)^2\)
我们通过这个损失函数来优化我们的线性模型 \(\beta\) 。此时我们的目标转化为
因为这是一个二次函数，我们可以直接通过梯度得到其解析解，即
\[\frac{\partial RSS}{\partial \beta} = -2X^T(Y-X\beta) = 0 \Rightarrow \beta = (X^TX)^{-1}X^TY\]此时，我们可以得到输出预测值为
\[\hat{Y} = X\beta = X(X^TX)^{-1}X^TY\]预测值 \(\hat{Y}\) 可以看成是真实值 Y 在 X 的列向量所构成的子空间上的投影，而残差 \(\epsilon\) 垂直与这个子空间。矩阵 \(H = X(X^TX)^{-1}X^T\) 用于计算正交投影。有时 X 中的数据不是线性独立的，则此时我们的模型 \(\beta\) 不是唯一表示的，我们可以通过重编码或去掉冗余数据来解决这一问题。在图片和信号处理中经常出现这一问题，我们通常用滤波降维，或者正则化控制拟合。
以上是我们对最小二乘的介绍及求解，接下来我们分析最小二乘为什么是对的。或者说，为什么最小二乘法这么流行。
最小二乘法实际上与极大似然估计是相联系的。我们假设误差服从高斯分布，即 \(\epsilon \sim N(0,1)\) ，结合之前的关系，我们有
\[y_i^*-x_i^T\beta=\epsilon_i \Rightarrow P(y_i^*\vert x_i)=p(\epsilon_i)=\frac{1}{\sqrt{2\pi}}exp(-(y_i^*-x_i^T\beta)^2/2)\]此时，我们化简极大似然的公式
\[logP(y_i^*\vert x_i)=-\frac{1}{2}(y_i^*-x_i^T\beta)^2 + constant\]所以我们可以看出，当我们最小化最小二乘loss时，相当于我们在最大化极大似然。即当误差服从高斯分布时，最小二乘与最大似然等价。
那为什么极大似然又是正确的呢？
这里我们用 \(P_{data}(x,y)\) 表示世界上无穷的真实分布，我们的训练集 \((x_i,y_i)\) 是从 \(P_{data}(x,y)\) 采样得到。首先介绍大数定律，当采样数足够大时，我们有
根据极大似然的定义，结合大数定律，我们可以写出其数学形式
\[max_{\theta} L(\theta) \equiv max_{\theta}\frac{1}{n}\sum_{i=1}^nlogP_{\theta}(y_i\vert x_i) \equiv E_{P_{data}}[logP_{\theta}(y\vert x)]\]此时我们可以在优化目标中加上一项 \(E_{P_{data}}[logP_{data}(y\vert x)]\) ，这一项与 \(\theta\) 模型无关，对模型优化无影响。其中 \(P_{data}(y\vert x)\) 表示ground-truth 的条件概率（上帝角度）。加入该项后，我们对优化目标改写并化简可得到
\[max_{\theta} L(\theta) \equiv min_{\theta} \{ -E_{P_{data}}[logP_{\theta}(y\vert x)] + E_{P_{data}}[logP_{data}(y\vert x)] \}\]我们可以看到，这就是模型 \(P_{\theta}(y\vert x)\) 与 ground-truth \(P_{data}(y \vert x)\) 的 KL 距离。
对于判别分类模型来说，极大似然估计的本质是最小化模型与真实之间的 KL 距离
尽管我们在上面已经得到了最小二乘对应的模型解析解，但当数据集非常大时，这样的计算消耗很大，为此我们需要采用在线算法，逐步地去迭代更新我们的模型。在这里我们考虑 stochastic gradient descent SGD 方法，每次选取 n 个数据，对这 n 个数据的Loss \(L_n\)求和，并更新模型，更新方式如下
\[\beta^{(t+1)} = \beta^{(t)} - \eta\nabla L_n\]这里 t 为迭代次数，\(\eta\) 为学习率，模型一开始被设为初值 \(\beta^{(0)}\) ，具体到线性模型中，在采用最小二乘 Loss 时，对应的更新公式为
\[\beta^{(t+1)} = \beta^{(t)} + \eta \sum_{i=1}^n(y_i^*-x_i^T\beta^{(t)})x_i\]线性模型下的梯度下降算法又被称为 least-mean-squares or the LMS algorithm 。这里的学习率需要仔细选取来保证算法收敛。
为了抑制过拟合的现象，我们可以在损失函数中添加正则项。这里我们主要介绍两种正则方式，分别是 Ridge Regression 和 Least Absolute Shrinkage and Selection Operator (LASSO)
在 Ridge Regression 下，我们改写最小二乘的loss函数为
\[L(\beta) = ||Y-X\beta||^2 + \lambda||\beta||^2\]其中 \(\lambda\vert\vert\beta\vert\vert^2\) 为正则项。这种正则方式类似于枪打出头鸟，惩罚明显特征，迫使模型关注更多的其他特征。添加正则项后，损失函数仍然是二次方程，因此可以利用和上面相同的方法直接得到模型的解析解，即
\[\beta_{\lambda} = (X^TX+\lambda I_p)^{-1}X^TY\]在 Lasso 下，我们改写最小二乘的loss函数为
\[L(\beta) = ||Y-X\beta||^2 + \lambda||\beta||_{l1},\quad where \ ||\beta||_{l1}=\sum_{j=1}^n|\beta_j|\]这种正则方式的意义是让模型不能每个feature都用。当输入数据为高维时没有解析解，当其为一维时，可得
\[\beta_{\lambda} = sign(\gamma)max(0,|\gamma|-\lambda/||X||^2), \quad where \ \gamma = (Y,X)/||X||^2\]和 Ridge Regression 相比
那么为什么 L1 正则会导致稀疏呢？
这是因为利用 Lagrange 求解时，L1 对权重的限制空间为矩形，而 L2 为圆形。由于 Lasso 对模型 \(\beta\) 的限制是有棱角的，因此解更容易相切在某一维为0的点
在这里我们来分析一下 Loss function 中的 bias 与 variance
在 squared loss function 下，理论上最好的模型是条件分布的期望。我们用 \(h(x)^{\star}\) 表示理论上最好的模型，可以将其表示为
\[h(x)^* = E[y\vert x] = \int yP(y \vert x)dy\]接下来我们写出 squared loss 的期望值，记为 E[L]
\[E[L] = \int \{f(x)-h(x)^*\}^2p(x)dx +\int \{h(x)^* - y^*\}^2p(x,y)dxdy\]这里我们实际建模得到的模型是 \(f(x)\) ，首先关注上式中的第二项，它与我们建立的模型无关。它的意义是理论上最好的模型作出的预测误差，这来源于数据自身内部的噪声，我们对此无能为力。理论上第二项也是我们能得到的最小误差。
而第一项则取决于我们实际建模的模型，我们的优化也是为了找到使第一项最小的模型。换句话说，如果我们有无限的数据，无限的计算资源，那么我们就可以建模出理论最优模型 \(h(x)^*\) ，从而使第一项完全消除。
接下来我们进一步考虑第一项。在实际中，我们只能通过一个具体的数据集 D 来训练模型。所以在这里我们用 \(f(x;D)\) 代替上式中的 \(f(x)\) ，表示通过数据集 D 训练得到的模型。如果我们对世界上所有的数据集求期望，可改写第一项为
\[E_D[\{f(x;D)-h(x)^*\}^2] = \{E_D[f(x;D)]-h(x)^*\}^2+E_D[\{f(x;D)-E_D[f(x;D)]\}^2]\]此时，上式中的第一项即为 bias 的平方，表示在世界上所有数据集的平均预测结果与理论最优模型之间的差距。上式中的第二项即为 variance，表示在某个特定数据集上的表现结果与在所有数据集上表现结果的差距。我们将其代入 squared loss 的期望值，可以得到
\[expected\ loss = (bias)^2 + variance + noise\]bias 和 variance 之间存在一种tradeoff的关系。如果我们的模型在所有数据集上都能较好地表现（即bias低），那么当其应用于某一特定的数据集时，其表现就会略差一些（即variance高）。相反，如果我们的模型只专门针对某个特定的数据集，那么可以学习得很好（即variance低），此时作用在其他数据集上，表现就会差一些（即bias高）。所以我们需要找一个最佳的模型去tradeoff bias and variance
以上都是理论上的分析，帮助理解 bias 与 variance 之间的矛盾。在实际中，我们只能得到一个观测数据集，无法去分析具体的 bias 与 variance 值。即使有多个数据集，我们也可以将其合并为一个数据集以实现更好的预测结果。
在上面的内容中，我们主要是从频率学派的观点来看待，即认为参数是固定的、未知的常量，我们只需要按某种方法去达到它，最后可以看作是一个优化问题。而接下来我们将从贝叶斯学派的观点来看待，即参数是未知的变量，它自身也是遵循某个概率分布的，我们只有它的先验分布，需要根据观察到的数据来进行调整。
和上面保持一致，对于给定数据x，其对应的ground-truth \(y^*\) 可以表示为
\[y^* = x^T\beta + \epsilon\]在这里，噪声 \(\epsilon \sim N(0,\sigma^2)\) 。和上面相同，我们有 \(y^* \sim N(x^T\beta,\sigma^2)\) 。唯一不同的是，在贝叶斯中我们认为模型 \(\beta\) 自身也是有一个概率分布的，我们只知道其先验分布为 \(\beta \sim N(0,\Sigma)\)
根据贝叶斯法则，我们写出模型参数的后验分布 \(P(\beta\vert Y,X)\)
\[P(\beta\vert Y,X) = \frac{P(X,Y,\beta)}{P(X,Y)} = \frac{P(X\vert\beta)}{P(X,Y)}P(Y\vert X,\beta)P(\beta)\]可以看出，后验正比于先验与似然之积，且似然与先验都是高斯分布，所以后验也是高斯分布。我们首先假设 \(\beta\vert Y,X \sim N(\mu',\Sigma')\) ，在下面的部分我们将求出后验分布的具体均值及方差。
因为我们已经确定后验服从高斯分布，所以我们只考虑 \(P(Y\vert X,\beta)P(\beta)\) 的指数部分就可以确定未知量。
\[P(Y\vert X,\beta)P(\beta) = Cexp(-\frac{1}{2\sigma^2}(Y-X\beta)^T(Y-X\beta))exp(-\frac{1}{2}\beta^T\Sigma^{-1}\beta)\]化简可得
\[P(Y\vert X,\beta)P(\beta) = Cexp(-\frac{1}{2}\beta^T(\sigma^{-2}X^TX+\Sigma^{-1})\beta + \sigma^{-2}Y^TX\beta -\frac{1}{2}\sigma^{-2}Y^TY)\]而我们将假设的 \(\beta\vert Y,X \sim N(\mu',\Sigma')\) 展开写出来可得到
\[P(\beta\vert Y,X) = C'exp(-\frac{1}{2}\beta^T(\Sigma^{'-1})\beta + \mu^{'T}\Sigma^{'-1}\beta -\frac{1}{2}\mu^{'T}\Sigma^{'-1}\mu'\]我们已经知道了正比关系 \(P(\beta\vert Y,X) \propto P(Y\vert X,\beta)P(\beta)\) ，所以将上面两个形式对比可得到
\[\Sigma^{'-1} = \sigma^{-2}X^TX+\Sigma^{-1}\]为对称矩阵，\((\Sigma^{'-1})^T=\Sigma^{'-1}\)，我们进一步求均值 \(\mu'\)
\[\Sigma^{'-1}\mu' = (\mu^{'T}(\Sigma^{'-1})^T)^T = (\sigma^{-2}Y^TX)^T = \sigma^{-2}X^TY\] \[\mu' = \sigma^{-2}\Sigma'X^TY\]至此，我们得到了后验分布 \(\beta\vert Y,X \sim N(\mu',\Sigma')\) 的参数，接着就能够根据给定的数据 x 进行预测了。因为模型 \(\beta\) 不再是存在一个固定的最优值，而是服从一个分布，因此 \(y_{pred}\) 也可以看作是一个随机变量的函数，进而得到其均值和方差
\[y_{pred} = x^T\beta \Rightarrow y_{pred} \sim N(x^T\mu',x^T\Sigma'x)\]