生成式模型、判别式模型、指数族分布、和广义线性模型

朴素贝叶斯 Naive Bayes classifier

逻辑回归 Logistic regression

指数族分布 Exponential family distributions

广义线性模型 Generalized linear models

本文分析生成式概率分类模型朴素贝叶斯和判别式概率分类模型逻辑回归的根源，关系，以及他们在广义线性模型框架下的统一。

（基于朱军老师的PPT和教科书《概率机器学习》）

特征空间分布

对于以上数据，为了尝试预测一个新数据属于红蓝两类分别的概率，我们先尝试获取两类数据的分布，当然，基于数据我们只能获取到两类的经验分布（empirical distribution），这里有两种方法：

1. 直方图（Histogram）
2. 核密度估计（Kernel density estimation, KDE），直方图（也就是经验数据）的平滑化（Smoothing）

核密度估计（KDE）

图中真实分布可以看作是两个正态分布的加权和

\[p(x) = 0.4 \times \mathcal{N}(0, 1) + 0.6 \times \mathcal{N}(6, 1)\]

要估计此分布，我们可以假设每个数据都是一个正态分布\(\mathcal{N}(0, 1)\)的期望值，叠加\(N\)个正态分布\(\mathcal{N}(0, 1)\)然后取平均：

\[p(\mathbf{x}) = \frac{1}{N} \sum_{n=1}^{N} k\left(\frac{\mathbf{x} - \mathbf{x}_n}{h}\right)\]

\[k(\mathbf{u}) = \mathcal{N}(\mathbf{0}, \mathbf{I})\]

可以看到，\(h\)（带宽，Bandwidth）决定了对于经验数据的拟合程度。

• \(h\)越小，越拟合数据，抖动越多
• \(h\)越大，越平滑，导致欠拟合

然后，有了特征空间的分布（的预测），就可以对新数据用上图占比分析两类的概率了。

生成式模型和判别式模型的基本介绍

这节先简单介绍两种模型是什么，与他们的相关性，后面会详细介绍他们的比较与优劣。

考虑垃圾邮件的例子，类别为垃圾邮件和正常邮件，特征为关键词。

生成式模型（Generative Model）

生成式模型学习（不一定真的进行生成这个行为，而是学习生成）每个类别的数据是如何“生成”的。

如，生成式模型朴素贝叶斯（符号参考下一节）会学习\(p(x|y)\)，给定类别\(y\)，每个特征\(x_i\)的概率分布是怎么样的。

这样，当考虑新输入（邮件）时，模型会判断这封邮件更像用垃圾邮件的写作模式写作（生成）出来的，还是正常邮件。

严格来讲，生成式模型假设联合分布\(P(X,Y)\)（或者分布\(P(Y)\)和\(P(X|Y)\)）并从训练数据尝试拟合分布，然后用最有可能的类别作为输出： \[ \hat{y} = \underset{y}{\mathrm{argmax}} \, P(\mathbf{x}, Y = y) \\ \hat{y} = \underset{y}{\mathrm{argmax}} \, P(Y = y \mid \mathbf{x}) \]

判别式模型（Discriminative Model）

判别式模型学习一个决策边界。直接对后验\(p(y|x)\)建模（“后验”参考下一小节），直接判断特征输入什么时候更像垃圾邮件，什么时候更像写作模式。

严格来讲，与其像生成式模型一样给出间接输出，判别式模型直接假设分布\(P(Y|X)\)并训练，最后直接给出输出。

相关性

根据贝叶斯公式，

\[p(y \mid x) = \frac{p(y,x)}{p(x)}=\frac{p(y) }{p(x)}p(x \mid y)\]

可以看到判别式模型直接对左边进行建模，而生成式模型对右边进行建模（\(p(x)\)固定，不需要建模）。

生成式模型：朴素贝叶斯 （Naive Bayes, NB）

上图中\(Y\)为类别，\(x_i\)为特征。NB假设同类下的特征彼此条件性独立（给定类标签）：

\[p(x_1,x_2,x_3,x_4|y)=p(x_1|y)p(x_2|y)p(x_3|y)p(x_4|y)\]

这也是朴素贝叶斯中朴素的来源。

给定所有特征\(x_i\)，朴素贝叶斯（NB）反推出这个数据属于哪个类别。也就是，NB通过后验概率\(p(Y|x)\)“生成”所属类别：

\[p(y \mid x) = \frac{p(y,x)}{p(x)}=\frac{p(y)\times p(x \mid y) }{p(x)}\]

\[y^* = \arg\max_{y\in\mathcal{Y}} p(y|x)\]

• \(p(y)\)是先验（prior）概率，也就是考虑具体特征\(x=x_1,...,x_4\)之前，对于类别\(y\)发生的可能性的“固有/先前”信念。
• \(p(x\mid y)\)是似然（likelihood），也是NB生成能力的体现（通过训练模型）。似然就是给定类别，各特征分别的分布是什么（根据朴素贝叶斯假设类内特征彼此条件独立）。
• \(p(y | x)\)是后验（posterior）概率，也就是在考虑特征后，对于\(y\)发生的可能性的修正概率。

因为我们只考虑两个\(p(y \mid x)\)哪个概率更高，所以可以忽略\(p(x)\)。

也就是，NB通过以上贝叶斯公式描述了一个信念（Belief）更新的过程：

1. 先学习一个先验信念\(p(y)\)，也就是垃圾邮件和正常邮件各自的占比。
2. 然后学习似然\(p(x|y)\)。
3. 当新数据\(X\)出现时，通过贝叶斯公式求出后验概率\(p(y|x)\)。
4. 由于我们只考虑哪个\(p(y|x)\)概率最大，所以可以忽略相同系数\(1/p(x)\)。

更概括的去说，NB（生成式分类器）学习联合分布\(p(y,x)\)来得到\(p(y \mid x)\)。

贝叶斯错误

\[p(error|\mathbf{x}) = \begin{cases} p(y=1|\mathbf{x}) & \text{if we decide } y=0 \\ p(y=0|\mathbf{x}) & \text{if we decide } y=1 \end{cases}\]

\[p(error) = \int_{-\infty}^{\infty} p(error|\mathbf{x})p(\mathbf{x})d\mathbf{x}\]

可以看到NB最大化后验概率等价于最小化错误率（对于每个\(x\)都最小化错误率），所以在所有分类其中NB必然是总（平均）错误率\(p(error)\)最低的。

我们定义贝叶斯错误率为此错误率，也就是最优分类器能达到的最低（平均）错误率。

可以看到，贝叶斯错误率的来源不是模型不够好，而是数据本身具有固有的不确定性。

类别和特征的概率分布

然而，类别和特征的真实分布模型并不知道，所以需要通过假设某些分布来近似真实分布。

二元特征NB

假设类别和特征（\(d\)个）都是伯努利分布（二元，值只取\(0,1\)）。

通过朴素贝叶斯假设，

\[p(x_1,x_2,x_3,x_4|y)=p(x_1|y)p(x_2|y)p(x_3|y)p(x_4|y)\]

学习联合分布\(p(x,y)=p(y)p(x|y)\)，相当于学习\(p(y)\)和每个\(p(x_i|y)\)。

那么，设定似然为：

\[p(x_j|y=0, \mathbf{q}) = \begin{cases} q_{0j} & \text{if } x_j = 1 \\ 1 - q_{0j} & \text{otherwise} \end{cases}\]

\[p(x_j|y=1, \mathbf{q}) = \begin{cases} q_{1j} & \text{if } x_j = 1 \\ 1 - q_{1j} & \text{otherwise} \end{cases}\]

先验为：

\[p(y|\pi) = \begin{cases} \pi & \text{if } y = 1 \text{ (i.e., bird)} \\ 1 - \pi & \text{otherwise} \end{cases}\]

总共需要学习\(2d+1\)个参数。

如果我们没有朴素贝叶斯假设的话，可以看到需要学习\(2^{d+1}\)指数级个参数。

最大似然估计（Maximum likelihood estimation, MLE）

我们可以用经验数据，通过最大似然估计（MLE）来学习调整上小节的参数。

假设有（\(N\)个数据的）数据集\(\{x_i,y_i\}_{i\le N}\)，我们的目标是让我们的伯努利参数\(\pi, q\)最大化经验数据发生的概率。也就是，寻找哪组参数让样本看起来最合理，最有可能发生。

假设所有数据都是i.i.d.的，那么

\[p(\{\mathbf{x}_i, y_i | \pi, \mathbf{q}\}) = \prod_{i=1}^{N} p(\mathbf{x}_i, y_i | \pi, \mathbf{q})\]

最佳参数为：

\[(\hat{\pi}, \hat{\mathbf{q}}) = \arg \max_{\pi, \mathbf{q}} p(\{\mathbf{x}_i, y_i\} | \pi, \mathbf{q})\]

但是，乘积远没有和好计算。于是，根据\(\log\)函数的单调递增性质，最大化上面乘积等价于最大化其log（假设不存在出现概率为0的数据）。

\[\log p(\{\mathbf{x}_i, y_i | \pi, \mathbf{q}\}) = \sum_{i=1}^{N} \log p(\mathbf{x}_i, y_i | \pi, \mathbf{q})\]

\[(\hat{\pi}, \hat{\mathbf{q}}) = \arg \max_{\pi, \mathbf{q}} \log p(\{\mathbf{x}_i, y_i\} | \pi, \mathbf{q})\]

于是，通过求导并等于0的方式，可以求出最佳参数为：

\[\hat{\pi} = \frac{N_1}{N} \quad \hat{q}_{0j} = \frac{N_0^j}{N_0} \quad \hat{q}_{1j} = \frac{N_1^j}{N_1}\]

\[N_k = \sum_{i=1}^{N} \mathbf{I}(y_i = k) \ : \ \text{类别为k的数据量} \]

\[N_k^j = \sum_{i=1}^{N} \mathbf{I}(y_i = k, x_{ij} = 1) \ : \ \text{类别k内特征为j的数据量}\]

\(\hat{\pi}\)的推导。

考虑\(\sum_{i=1}^{N} \log p(\mathbf{x}_i, y_i | \pi, \mathbf{q})\)，对于每个联合分布，拆分成\(p(\mathbf{x}_i, y_i | \pi, \mathbf{q})=p(y_i | \pi, \mathbf{q})p(\mathbf{x}_i | y_i, \pi, \mathbf{q})\)。

剔除不相关参数，变为\(p(y_i | \pi)p(\mathbf{x}_i | y_i, \mathbf{q})\)。

取\(\log\)，乘积拆分为和，所有样本的第一项第二项分别组在一起，由于第二项对于\(\pi\)是常数，对\(\pi\)求偏导其消失，我们只考虑第一项。

第一项因\(p(y_i|\pi)\)为伯努利，等于\(\pi^{y_i}(1-\pi)^{1-y_i}\)，取\(\log\)为\(y_i\log\pi+(1-y_i)\log(1-\pi)\)。

于是，对\(\sum_{i=1}^{N} \log p(\mathbf{x}_i, y_i | \pi, \mathbf{q})\)求\(\pi\)的偏导并设为0，变为 \[ \begin{aligned} \frac{\partial}{\partial \pi} \sum_{i=1}^{N} \left[ y_i \log \pi + (1-y_i) \log(1-\pi) \right] &= 0 \\ \sum_{i=1}^{N} \frac{\partial}{\partial \pi} \left[ y_i \log \pi + (1-y_i) \log(1-\pi) \right] &= 0 \\ \sum_{i=1}^{N} \left[ y_i \cdot \frac{1}{\pi} + (1-y_i) \cdot \frac{1}{1-\pi} \cdot (-1) \right] &= 0 \\ \sum_{i=1}^{N} \left[ \frac{y_i}{\pi} - \frac{1-y_i}{1-\pi} \right] &= 0 \\ \frac{1}{\pi} \sum_{i=1}^{N} y_i - \frac{1}{1-\pi} \sum_{i=1}^{N} (1-y_i) &= 0 \\ \frac{N_1}{\pi} - \frac{N_0}{1-\pi} &= 0 \\ \frac{N_1}{\pi} &= \frac{N_0}{1-\pi} \\ N_1 (1-\pi) &= N_0 \pi \\ N_1 - N_1 \pi &= N_0 \pi \\ N_1 &= N_0 \pi + N_1 \pi \\ N_1 &= (N_0 + N_1) \pi \\ N_1 &= N \pi \\ \hat{\pi} &= \frac{N_1}{N} \end{aligned} \]

MLE等价于最小化KL散度

KL散度的概念参考这里。

MLE可以看到就是在最小化经验分布和模型分布之间的KL散度：

首先考虑三个相关分布；

1. 数据的真实分布。模型无法直接得知，需要通过数据来近似，模拟。
2. 经验分布。观测到的数据代表的分布，当数据没有偏见时，随着样本量增加，经验分布逼近真实分布。
3. 模型分布。用参数\(\pi, q\)来近似，表示经验分布。

数据固定的情况真实分布和经验分布无法改变，只能让模型分布逼近经验分布来尝试对真实分布进行逼近。

在信息理论中，衡量分布是否“接近”就可以通过KL散度来判断。

我们现在证明MLE的最优参数等价于最小化KL散度的参数，也就是：

\[\theta^* = \arg\min_{\theta} D_{KL}(P_{true} \| Q_{\theta})\]

这里\(P_{true}\)是真实分布，\(Q_\theta\)是模型分布，\(\theta=(\pi,q)\)。

---

首先，

\[D_{KL}(P_{true} \| Q_{\theta}) = \int P_{true}(x) \log\left(\frac{P_{true}(x)}{Q_{\theta}(x)}\right) dx\]

拆分\(\log\)，

\[D_{KL}(P_{true} \| Q_{\theta}) = \underbrace{\int P_{true}(x) \log(P_{true}(x)) dx}_{\text{Term 1}} - \underbrace{\int P_{true}(x) \log(Q_{\theta}(x)) dx}_{\text{Term 2}}\]

第一项是\(P_{true}\) 的负熵 (negative entropy)，对于参数\(\theta\)是常数，忽略。

第二项是他们的交叉熵，也就是

\[E_{x \sim P_{true}}[\log(Q_{\theta}(x))]\]

即 \(\log(Q_{\theta}(x))\) 在真实分布下的期望。关于交叉熵信息论的定义参考上面超链接。

所以，最小化KL散度等价于最大化

\[E_{x \sim P_{true}}[\log(Q_{\theta}(x))]\]

我们不知道真实分布（true distribution），担可以通过经验分布（empirical distribution）模拟、近似真实分布。

\[E_{x \sim P_{true}}[\log(Q_{\theta}(x))]\approx E_{x \sim P_{emp}}[\log(Q_{\theta}(x))]\]

但，对于经验分布来说，每个数据点\(x\)发生的概率就是\(\frac1n\)，假设有\(n\)个数据点。

\[E_{x \sim P_{emp}}[\log(Q_{\theta}(x))] = \sum_{i=1}^n P_{emp}(x_i) \log(Q_{\theta}(x_i)) = \sum_{i=1}^n \frac{1}{n} \log(Q_{\theta}(x_i))\]

于是，最大化原始目标又变为了最大化（不考虑常数\(\frac1n\)）

\[\sum_{i=1}^n\log(Q_{\theta}(x_i))\]

这刚好对应了MLE最大化的log-似然（log-likelihood）：

\[\log p(\{\mathbf{x}_i, y_i\} | \pi, \mathbf{q})\]

贝叶斯决策里的似然与MLE的似然不一样，语境不同。

所以，他们等价。

---

拉普拉斯平滑

MLE的最佳参数之前提到了为

\[\hat{\pi} = \frac{N_1}{N} \quad \hat{q}_{0j} = \frac{N_0^j}{N_0} \quad \hat{q}_{1j} = \frac{N_1^j}{N_1}\]

但是，\(N_0, N_1\)有可能为0（数据稀疏性，也就是零计数问题）。

为了保证分母不为0，我们用拉普拉斯平滑

\[\hat{q}_{0j} = \frac{N_0^j + \alpha}{N_0 + 2\alpha}\]

\[\hat{q}_{1j} = \frac{N_1^j + \alpha}{N_1 + 2\alpha}\]

其中\(\alpha>0\)。

拉普拉斯平滑看起来是凭空捏造的，但其实是下一小节贝叶斯估计的一个特例。

从MLE，拉普拉斯平滑到最大后验估计

首先，正如上小节拉普拉斯平滑中提到的，MLE的缺陷就在于其本质上就是“数频率”，对于数据稀疏（部分特征一次都没出现过）的情况，MLE会崩溃。

为了解决\(q_{ij}\)无法被定义的问题，我们可以给所有\(q_{ij}\)一个”初始信念”（假设一个先验分布）。然后，根据观测数据来更新我们的信念，得到一个后验分布。最后，从后验分布中取概率最大的点为\(q_{ij}\)的估计（也叫做最大后验估计，Maximum a Posteriori Estimate，MAP)）。 想法类似于NB信念更新的过程。

也就是，我们需要一个“先验”，是为了避免 MLE 在数据不足时做出“非 0 即 1”的极端（过拟合）估计。

Beta分布

我们选取Beta分布作为参数的先验（对于参数\(q_{1j}\)）。

\[p_0(q_{1j} | \alpha_1, \alpha_2) = \text{Beta}(\alpha_1, \alpha_2) = \frac{\Gamma(\alpha_1 + \alpha_2)}{\Gamma(\alpha_1)\Gamma(\alpha_2)} q_{1j}^{\alpha_1 - 1} (1 - q_{1j})^{\alpha_2 - 1}\]

选Beta有几个原因：

1. \(q_{1j}\)作为伯努利的参数取值\([0,1]\)，而Beta的定义域（Domain）刚好匹配。

2. 共轭先验（Conjugate Prior）。正如上面一条提到的，\(q_{0j}\)来自伯努利过程（Bernoulli process），多个伯努利过程的叠加（求和）就是二项分布（Binomial distribution），似然的数学形式对应了\(L(q) \propto q^k (1-q)^{n-k}\)。

   这与上方Beta分布的概率密度函数（PDF）展开后的形式完全相同。

   \[ \begin{aligned} \text{Posterior} &\propto \text{Likelihood} \times \text{Prior} \\ p(q \mid \text{Data}) &\propto \left( q^{N_1^j} (1-q)^{N_1 - N_1^j} \right) \times \left( q^{\alpha_1-1} (1-q)^{\alpha_2-1} \right) \\ p(q \mid \text{Data}) &\propto q^{(N_1^j + \alpha_1) - 1} \cdot (1-q)^{(N_1 - N_1^j + \alpha_2) - 1} \end{aligned}\]

   先验和似然相乘得到的后验仍然是Beta分布。当先验分布（Beta）和似然函数（Binomial/Bernoulli）结合后，能得到一个与先验同类型（也是Beta）的后验分布时，我们称这个先验（Beta）是这个似然（Bernoulli）的“共轭先验”。

Beta作为共轭先验最大的好处就是信念更新仅需要做加法：

• 先验信念是 \(Beta(\alpha_1, \alpha_2)\)。
• 观测到 \(N_1^j\) 次成功（\(x_j=1\)）和 \((N_1 - N_1^j)\) 次失败（\(x_j=0\)）。
• 后验信念自动更新为 \(Beta(\alpha_1 + N_1^j, \alpha_2 + (N_1 - N_1^j))\)。

贝叶斯估计总结

考虑更新参数\(q_{1j}\)：

1. 用 \(Beta(\alpha_1, \alpha_2)\) 作为先验。
2. 观测 \(N_1^j\) 次 \(x_j=1\) 和 \((N_1 - N_1^j)\) 次 \(x_j=0\)（总共 \(N_1\) 次观测）。
3. 后验分布更新为 \(Beta(N_1^j + \alpha_1, N_1 - N_1^j + \alpha_2)\)。
4. MAP估计为新Beta分布的众数（mode，即峰值点）\(\frac{A-1}{A+B-2}\)。
5. 代入 \(A = N_1^j + \alpha_1\) 和 \(B = N_1 - N_1^j + \alpha_2\)得到最终后验估计：

\[\hat{q}_{1j} = \frac{(N_1^j + \alpha_1) - 1}{(N_1^j + \alpha_1) + (N_1 - N_1^j + \alpha_2) - 2} = \frac{N_1^j + \alpha_1 - 1}{N_1 + \alpha_1 + \alpha_2 - 2}\]

类似的，参数\(q_{0j}\)更新为：

\[\hat{q}_{0j} = \frac{N_0^j + \alpha_1 - 1}{N_0 + \alpha_1 + \alpha_2 - 2}\]

MLE和拉普拉斯平滑作为MAP估计的特例

MLE：\(\hat{q}_{1j} = \frac{N_{1}^{j}}{N_{1}}\)

拉普拉斯平滑：\(\hat{q}_{1j} = \frac{N_{1}^{j}+\alpha}{N_{1}+2\alpha}\)

MAP估计： \(\hat{q}_{1j} = \frac{N_{1}^{j}+\alpha_{1}-1}{N_{1}+\alpha_{1}+\alpha_{2}-2}\)

令\(\alpha_1=\alpha_2=1\)（对应了上图的绿线），MAP估计的最优参数退化为MLE。

解方程组：

\(\alpha = \alpha_1 - 1\)

\(2\alpha = \alpha_1 + \alpha_2 - 2\)

得到\(\alpha_1 = \alpha + 1\)， \(\alpha_2 = \alpha + 1\)。于是可以看到拉普拉斯平滑也为MAP估计的一个特例。

---

尽管MAP估计理论很泛用，现实中一般用的还是拉普拉斯平滑并仅考虑一个超参数\(\alpha\)。

高斯朴素贝叶斯（Gaussian Naive Bayes, GNB）

上面假设了类别和特征都是二元（Binary）的，通过简单的扩展可以处理所有离散（Discrete）类型的特征和类别。不仅如此，NB也可以扩展并处理连续的特征、类别。

高斯朴素贝叶斯假设了先验继续服从伯努利，而似然服从正态（高斯）分布： \[ Y \sim \text{Bernoulli}(\pi) \\ P(X_i \mid Y = y) = \mathcal{N}(\mu_{iy}, \sigma_{iu}^2) \] 每个参数可以用不同的均值和方差，但参数数量可能过多。

为此，可以假设：

• 方差与类别无关。\(\sigma_{i0}^2 = \sigma_{i1}^2 = \sigma_i^2\)。
• 方差与特征无关。\(\sigma_{1y}^2 = \sigma_{2y}^2 = \sigma_y^2\)。
• 方差与两者都无关。\(\sigma\)。

MLE估计参数

\[ \hat{\mu}_{ik} = \frac{1}{\sum_{n} \mathbb{I}(y_n = k)} \sum_{n} x_{ni} \mathbb{I}(y_n = k) \\ \hat{\sigma}_{ik}^2 = \frac{1}{\sum_{n} \mathbb{I}(y_n = k)} \sum_{n} (x_{ni} - \hat{\mu}_{ik})^2 \mathbb{I}(y_n = k) \\ h(\mathbf{x}) = \underset{y}{\mathrm{argmax}} \, P(y) \prod_{i} P(x_i \mid y) \]

和前面的MLE估计类似，可以看到估计的均值就是类别\(k\)和特征\(i\)样本的均值；估计的方差就是类别\(k\)和特征\(i\)样本的方差。

树增广朴素贝叶斯（Tree augmented NB）

除了类别和特征分布上的扩展，对于朴素贝叶斯假设也可以通过扩展模型结构来松化。如，使用树增强朴素贝叶斯（Tree augmented NB）：

我们松弛了朴素贝叶斯假设，此模型的参数也从原NB的\(2d+1\)变成了\(4d+1\)，还是线性，远好于指数。

贝叶斯模型的决策边界

NB看起来这么好用，错误率这么低，那么决策边界是否变为非线性了呢？

（同方差）高斯朴素贝叶斯 GNB

假设所有参数的正态分布方差与类别无关\(\sigma_{i0}^2 = \sigma_{i1}^2 = \sigma_i^2\)。

让\(X=X_1,...,X_d\)，决策边界是\(P(Y=0|X) = P(Y=1|X)\) 的地方，通过贝叶斯公式，这等价于\(P(X|Y=0)P(Y=0) = P(X|Y=1)P(Y=1)\)。

取对数几率，这等价于 \[ \log \frac{P(X|Y=0)P(Y=0)}{P(X|Y=1)P(Y=1)} = 0\\ \sum_{i=1}^d \left( \log P(X_i|Y=0) - \log P(X_i|Y=1) \right) + \log \frac{P(Y=0)}{P(Y=1)} = 0\\ \sum_{i=1}^d \left( -\frac{(x_i - \mu_{i0})^2}{2\sigma_i^2} - \left(-\frac{(x_i - \mu_{i1})^2}{2\sigma_i^2}\right) \right) + C = 0\\ \sum_{i=1}^d \frac{1}{2\sigma_i^2} \left( (x_i - \mu_{i1})^2 - (x_i - \mu_{i0})^2 \right) + C = 0\\ \sum_{i=1}^d \underbrace{\left(\frac{\mu_{i0} - \mu_{i1}}{\sigma_i^2}\right)}_{w_i} x_i + \underbrace{\left( \sum_{i=1}^d \frac{\mu_{i1}^2 - \mu_{i0}^2}{2\sigma_i^2} + C \right)}_{w_0} = 0 \] 上面代入了\(P(X_i|Y=y)=\mathcal{N}(\mu_{iy},\sigma_i^2)\)，\(\mathcal{N}(\mu, \sigma^2) \propto \exp\left(-\frac{(x - \mu)^2}{2\sigma^2}\right)\)，以及平方和展开 \[ \left( (x_i - \mu_{i1})^2 - (x_i - \mu_{i0})^2 \right) = (x_i^2 - 2x_i\mu_{i1} + \mu_{i1}^2) - (x_i^2 - 2x_i\mu_{i0} + \mu_{i0}^2)\\ = -2x_i\mu_{i1} + 2x_i\mu_{i0} + \mu_{i1}^2 - \mu_{i0}^2= 2x_i(\mu_{i0} - \mu_{i1}) + (\mu_{i1}^2 - \mu_{i0}^2) \] 哇，竟然是线性决策边界（无慈悲）。

一般多元高斯（非朴素）贝叶斯

假设我们抛弃朴素贝叶斯特征之间的条件独立性假设。那么，可以将（类内）所有特征合在一起看成一个多元高斯分布： \[ P_1 = P(Y = 0), \quad P_2 = P(Y = 1) \\ p_1(X) = p(X \mid Y = 0) = \mathcal{N}(M_1, \Sigma_1) \\ p_2(X) = p(X \mid Y = 1) = \mathcal{N}(M_2, \Sigma_2) \] 这里，我们需要用协方差矩阵\(\Sigma\)对所有特征对（\(X_i\), \(X_j\)）直接的关系进行建模。这也是对之前的朴素贝叶斯假设的一种泛化，可以看到假设所有特征彼此类内独立时，协方差矩阵\(\Sigma\)就退化为了对角线矩阵，并且联合分布\(X\)可以对于每个特征拆分成一个单元正态分布。

我们将证明如上图所示，如果两类的协方差矩阵相等，那么决策边界就是线性的；如果不相等，那么就是非线性的（二次的）。上一小节的同方差高速朴素贝叶斯就对应了协方差矩阵\(\Sigma_1=\Sigma_2=\Sigma\)相等的场景，我们也的确推导出了决策边界为线性。

和上一小节一样，我们考虑决策边界的等式，然后取\(\ln\)（而非\(\log\)） \[ P(Y=0|X) = P(Y=1|X)\\ P(X|Y=0)P(Y=0) = P(X|Y=1)P(Y=1)\\ \ln P(X|Y=0) + \ln P(Y=0) = \ln P(X|Y=1) + \ln P(Y=1) \] 然后，将\(P(X|Y=y)\)展开为多元高斯分布而非\(d\)个正态分布的乘积，并取\(\ln\)： \[ p(X | M_y, \Sigma_y) = \frac{1}{(2\pi)^{d/2} |\Sigma_y|^{1/2}} \exp\left(-\frac{1}{2}(X-M_y)^T \Sigma_y^{-1} (X-M_y)\right)\\ \ln p(X | M_y, \Sigma_y) = -\frac{1}{2}(X-M_y)^T \Sigma_y^{-1} (X-M_y) - \frac{1}{2}\ln|\Sigma_y| - C \] 已经可以看到第一项是\(X\)的二次项而非一次项。展开第一项： \[ (X-M_y)^T \Sigma_y^{-1} (X-M_y)= (X^T - M_y^T) \Sigma_y^{-1} (X - M_y)\\ = X^T \Sigma_y^{-1} X - X^T \Sigma_y^{-1} M_y - M_y^T \Sigma_y^{-1} X + M_y^T \Sigma_y^{-1} M_y\\ = X^T \Sigma_y^{-1} X - 2M_y^T \Sigma_y^{-1} X + M_y^T \Sigma_y^{-1} M_y \] 代入到刚开始我们推导的决策边界等式 \[ \ln P(X|Y=0) - \ln P(X|Y=1) + \text{Const} = 0 \] 我们得到： \[ -\frac{1}{2} \left( X^T \Sigma_0^{-1} X - 2M_0^T \Sigma_0^{-1} X + \dots \right)- \\ \left( -\frac{1}{2} \left( X^T \Sigma_1^{-1} X - 2M_1^T \Sigma_1^{-1} X + \dots \right) \right) + \text{Const} = 0 \] 整理一下： \[ \underbrace{-\frac{1}{2} X^T (\Sigma_0^{-1} - \Sigma_1^{-1}) X}_{\text{二次项 (Quadratic Term)}} + \underbrace{(M_0^T \Sigma_0^{-1} - M_1^T \Sigma_1^{-1}) X}_{\text{一次项 (Linear Term)}} + \underbrace{C'}_{\text{所有常数}} = 0 \] 这证明了一般多元高斯（非朴素）贝叶斯模型的决策边界为二次而非线性。

（同方差）高斯朴素贝叶斯等价于（逻辑回归的）Sigmoid函数

逻辑（Logistic）函数和Sigmoid函数的定义参考Logit, Logistic, Sigmoid, Softmax在MLP的区别 - alittlebear’s blog。

首先，一般多元高斯贝叶斯对于\(y=1\)的后验概率是： \[ p(y=1 \mid \mathbf{x}, \mu, \Sigma, \pi) = \frac{p(y=1, \mathbf{x} \mid \mu, \Sigma, \pi)}{p(\mathbf{x} \mid \mu, \Sigma, \pi)} \] 考虑朴素贝叶斯的情况，也就是协方差矩阵\(\Sigma\)为对角矩阵的情况，多元高斯分解为了单元高斯的乘积。

将\(p(\mathbf{x} \mid \mu, \Sigma, \pi)\)通过

\[p(x) = p(x|y=1)p(y=1) + p(x|y=0)p(y=0)\]

分解，分子分母同时除以\(p(x|y=1)p(y=1)\)，然后结合上面关于（同方差）高斯朴素贝叶斯的决策边界推导，可以得到 \[ \begin{aligned} p(y=1 \mid \mathbf{x}, \mu, \Sigma, \pi) &= \frac{1}{1 + \frac{p(y=0, \mathbf{x} \mid \mu, \Sigma, \pi)}{p(y=1, \mathbf{x} \mid \mu, \Sigma, \pi)}} \\ &= \frac{1}{1 + \frac{(1-\pi) \prod_i \mathcal{N}(x_i \mid \mu_{i0}, \sigma_i^2)}{\pi \prod_i \mathcal{N}(x_i \mid \mu_{i1}, \sigma_i^2)}} \\ &= \frac{1}{1 + \exp(-\mathbf{w}^\top \mathbf{x} - w_0)} \end{aligned} \] 神奇的变为了（下一小节逻辑回归拟合的）Sigmoid函数。

也就是说，（同方差）高斯朴素贝叶斯在形式上（模型表达能力）等价于下一节讲的逻辑回归模型（直接拟合逻辑Logistic函数，或者说Sigmoid函数作为输出），并且可以看到多元高斯（非朴素）贝叶斯的推广在模型表达能力上优于这两者。

判别式模型：逻辑回归（Logistic Regression，LR）

正如上一小节和开始关于判别式模型的描述，逻辑回归模型直接考虑输出概率为Logistic函数，或者说Sigmoid函数（参考上一节超链接）的输出： \[ P(y = 1 \mid \mathbf{x}) = \frac{1}{1 + \exp(-(w_0 + \mathbf{w}^\top \mathbf{x}))} \] 更一般的，可以通过添加参数\(\alpha\) \[ \psi_\alpha(v) = \frac{1}{1 + \exp(-\alpha v)} \] 来控制Sigmoid函数的陡峭程度。当\(\alpha\to\infty\)时Sgmoid函数变为阶梯函数（step function），阶梯函数对于一些神经网络可能更好，这里不细讲。

逻辑回归的决策边界

正如上一节提到的，逻辑回归的决策边界就是（同方差）高斯朴素贝叶斯的决策边界，线性的。

多类别逻辑回归

首先，为了简化表示，对于Sigmoid

\[ P(y = 1 \mid \mathbf{x}) = \frac{1}{1 + \exp(-(w_0 + \mathbf{w}^\top \mathbf{x}))} \]

我们增加输入和权重一个维度：

\[ \text{输入: } \begin{pmatrix} 1 \\ \mathbf{x} \end{pmatrix} \quad \text{模型权重: } \begin{pmatrix} w_0 \\ \mathbf{w} \end{pmatrix} \]

来省略\(w_0\)的显示：

\[P(y = 1 \mid \mathbf{x}) = \frac{1}{1 + \exp(-\mathbf{w}^\top \mathbf{x})}\]

现在，我们考虑对于多类别数据，如何扩展。

Softmax

扩展二元Sigmoid为Softmax： \[ P(Y=k|x) = \frac{\exp(w_k^\top x)}{\sum_{j=1}^{K} \exp(w_j^\top x)} \] 这保证了所有\(K\)个类别的概率之和为1。

最大条件似然估计（Maximum Conditional Likelihood Estimate, MCLE）

这里我们仅考虑二元类别的情况（Sigmoid而不是Softmax）。

因为逻辑回归不考虑联合分布\(P(X,Y)\)（或者\(P(Y)\)和\(P(X|Y)\)），无法通过MLE的方法找到最优权重\(w\)。

所以，我们直接考虑\(P(Y|X)\)并用其找到最优权重： \[ \hat{w} = \arg\max_w \prod_{i=1}^N P(y_i | x_i, w) \] 于NB的MLE相同，乘积不好优化，我们对似然（Likelihood）取Log并找对数似然的最优解，他们是相等的。

通过 \[ P(y|x) = P(y=1|x)^y \cdot P(y=0|x)^{1-y} \] 我们得到： \[ \begin{aligned} \mathcal{L}(\mathbf{w}) &= \log \prod_{i=1}^{N} P(y_i \mid \mathbf{x}_i, \mathbf{w}) \\ &= \sum_{i} [y_i \mathbf{w}^\top \mathbf{x}_i - \log(1 + \exp(\mathbf{w}^\top \mathbf{x}_i))]\\ \hat{\mathbf{w}} &= \underset{\mathbf{w}}{\mathrm{argmax}} \, \mathcal{L}(\mathbf{w}) \end{aligned} \] 然而，这个函数不存在闭式解（或者说解析解），但是是凹函数（某些中文教材的凹凸定义相反），所以可以用凸优化的方法（迭代）找到（全局）最优解。

优化算法

凸优化最常用的算法就是梯度上升（下降）。

梯度上升（Gradient Ascent）

由于 \(\mathcal{L}(w)\) 是凹函数，我们可以使用梯度上升法来寻找最大值。

考虑梯度： \[ \nabla_{\mathbf{w}} \mathcal{L}(\mathbf{w}) = \begin{pmatrix} \frac{\partial \mathcal{L}(\mathbf{w})}{\partial w_0} \\ \vdots \\ \frac{\partial \mathcal{L}(\mathbf{w})}{\partial w_d} \end{pmatrix} \] 然后用梯度更新权重： \[ \mathbf{w}_{t+1} = \mathbf{w}_t + \eta \nabla_{\mathbf{w}} \mathcal{L}(\mathbf{w}) |_{\mathbf{w}_t} \] 这里模型超参数\(\eta\)过小会导致收敛速度过慢（残余错误 residual error不大）；过大收敛会变快但残余错误过大，可能导致迭代时一直震荡无法收敛。

更具体的，Sigmoid函数的梯度通过推导可以看到十分简洁， \[ \nabla_w \mathcal{L}(w) = \nabla_w \sum_{i=1}^{N} \left[ \dots \right] = \sum_{i=1}^{N} \nabla_w \left[ y_i w^{\top}x_i - \log(1+\exp(w^{\top}x_i)) \right] \] 第一项的梯度为 \[ \nabla_w (y_i w^{\top}x_i) = y_i \nabla_w (w^{\top}x_i) = y_i x_i \] 第二项的梯度通过链式法则求出。让\(u = w^{\top}x_i\)，定义\(\mu_i = P(y=1|x_i, w_t)\)，得出 \[ \begin{aligned} \frac{d}{du} \left( -\log(1+\exp(u)) \right) &= - \frac{1}{1+\exp(u)} \cdot \frac{d}{du}(1+\exp(u))\\ &= - \frac{1}{1+\exp(u)} \cdot \exp(u)\\ &= - \frac{\exp(u)}{1+\exp(u)}\\ \mu_i := P(y=1|x_i, w) &= \frac{1}{1+\exp(-w^{\top}x_i)} = \frac{1}{1+\exp(-u)}\\ &= \frac{1 \cdot \exp(u)}{(1+\exp(-u)) \cdot \exp(u)} = \frac{\exp(u)}{\exp(u)+1}\\ \frac{d}{du} \left( -\log(1+\exp(u)) \right)&=-\mu_i\\ \nabla_w (-\log(1+\exp(w^{\top}x_i))) &= \left( - \mu_i \right) \cdot \nabla_w u = - \mu_i x_i \end{aligned} \] 合并两项，梯度简化为 \[ \nabla_w \mathcal{L}(w) = \sum_i (y_i - \mu_i) x_i \] 这也十分符合直觉。\((y_i - \mu_i^t)\)是第\(i\)个样本的真实值与预测概率之间的误差。算法让权重\(w\)沿着误差的方向更新，以减小误差。

梯度上升的问题

梯度上升（下降）作为凸优化最简单的算法，其收敛速度不是最优的。

如，下面几小节会提到的，可以通过牛顿法（Newton method），共轭梯度上升（conjugate gradient ascent），迭代重加权最小二乘（iterative reweighted least squares，IRLS，一会可以看到对于逻辑回归这等价于牛顿法）来加快收敛速度，减少所需迭代次数。

此外，权重\(w\)某些数值过大往往代表了模型过拟合，需要通过正则化来惩罚过大权重数值。

牛顿法（Newton’s Method，Newton-Raphson Method）

经典牛顿法通过一下更新来迭代找到多项式的根\(f(x)=0\)： \[ x_{t+1}=x_t-\frac{f(x_t)}{f'(x_t)} \] 对于我们的问题，为了最大化对数似然 \[ \mathcal{L}(\mathbf{w}) = \sum_{i} [y_i \mathbf{w}^\top \mathbf{x}_i - \log(1 + \exp(\mathbf{w}^\top \mathbf{x}_i))] \] 我们需要找到权重\(w^\star\)最大化\(\mathcal{L}(w)\)，也就是，让梯度为零（函数为凹函数所以就是全局最大）： \[ \nabla \mathcal{L}(\mathbf{w}^*) = 0 \] 把经典牛顿法套过来，迭代过程就变为了 \[ \mathbf{w}_{t+1} \leftarrow \mathbf{w}_t - H^{-1} \nabla_{\mathbf{w}} \mathcal{L}(\mathbf{w}) |_{\mathbf{w}_t} \] 其中\(H\)是Hessian矩阵（二阶导数矩阵） \[ H = \nabla_{\mathbf{w}}^2 \mathcal{L}(\mathbf{w}) |_{\mathbf{w}_t} \] 沿用梯度上升小节的定义\(\mu_i\)，可以得到具体形式（推导已省略）： \[ \nabla_{\mathbf{w}} \mathcal{L}(\mathbf{w}) |_{\mathbf{w}_t} = \sum_{i} (y_i - \mu_i) \mathbf{x}_i = X(\mathbf{y} - \mathbf{\mu}) \\ \mu_i = \psi(\mathbf{w}_t^\top \mathbf{x}_i)\\ H = \nabla_{\mathbf{w}}^2 \mathcal{L}(\mathbf{w}) |_{\mathbf{w}_t} = - \sum_{i} \mu_i (1 - \mu_i) \mathbf{x}_i \mathbf{x}_i^\top = - X R X^\top \\ R_{ii} = \mu_i (1 - \mu_i) \] 可以看到牛顿法（以及下面的ILRS）都是二阶方法：使用二阶导数（Hessian矩阵）来更快地逼近最优点。

逻辑回归的牛顿法也是迭代重加权最小二乘法（iterative reweighted least squares，ILRS）

上面的牛顿法可以看作是最大化分类任务的条件似然估计（MCLE），通过简单的变换，我们可以发现牛顿法某种意义上也是在最小化平方误差（某种意义上的最小二乘法），也是回归任务的优化目标。

这也给了我们一些启发，分类模型逻辑回归（或者说决策边界为线性的模型）是否类似回归模型线性回归（两者的区别可以参考这篇文章）。本文最后的广义线性模型部分就探讨了此问题。

标准线性回归最小二乘法的闭式解

线性回归模型最小化均方误差（mean square error, MSE），等价于最小化平方误差（squared error）

\[L(w) = ||y - \hat{y}||^2 = ||y - X^\top w||^2\]

展开变为 \[ \begin{aligned} L(w) &= (y - X^\top w)^\top (y - X^\top w) \\ &= (y^\top - w^\top X) (y - X^\top w) \\ &= y^\top y - y^\top X^\top w - w^\top X y + w^\top X X^\top w \\ &= y^\top y - 2 w^\top X y + w^\top (X X^\top) w\end{aligned} \] 求梯度并设其为0来找到最优参数（权重）\(w\)（可以看到我们可以直接求出闭式解/解析解） \[ \nabla_w L(w) = 0 - 2 X y + 2 (X X^\top) w=0\\ - 2 X y + 2 (X X^\top) w = 0\\ 2 (X X^\top) w = 2 X y\\ (X X^\top) w = X y\\ (X X^\top)^{-1} (X X^\top) w = (X X^\top)^{-1} X y\\ w = (X X^\top)^{-1} X y \] 这就是最优参数。

逻辑回归的牛顿法是迭代重加权最小二乘法

沿用牛顿法小节的符号，变换牛顿法权重迭代公式， \[ \begin{aligned} \mathbf{w}_{t+1} &= \mathbf{w}_t - H^{-1} \nabla_{\mathbf{w}} \mathcal{L}(\mathbf{w}_t) \\ &= \mathbf{w}_t - (X R X^\top)^{-1} X (\mathbf{\mu} - \mathbf{y}) \\ &= (X R X^\top)^{-1} \{ X R X^\top \mathbf{w}_t - X (\mathbf{\mu} - \mathbf{y}) \} \\ &= (X R X^\top)^{-1} X R \mathbf{z}\\ \mathbf{z} &= X^\top \mathbf{w}_t - R^{-1}(\mathbf{\mu} - \mathbf{y}) \end{aligned} \] 可以看到，标准最小二乘的闭式解是 \[ w = (XX^\top)^{-1} X y \] 变形后的牛顿法迭代是 \[ w_{t+1} = (X \mathbf{R} X^\top)^{-1} X \mathbf{R} z \] 形式完全就是：

• 迭代（iterative）。\(w_t \to w_{t+1} \to w_{t+2} \dots\)。
• 重加权（reweighted）。\(R\)每个迭代都更新。
• 最小二乘（least squares）。形式如上是最小二乘问题。

所以，虽然牛顿法是在最大化条件似然估计，但是对于逻辑回归来说，也是在最小化平方误差，等价于迭代重加权最小二乘法。

梯度上升和牛顿法/ILRS对比

上图横轴为迭代次数，纵轴为错误率；红线为IRLS（也就是牛顿法），蓝线为梯度上升（实际为梯度下降）。

可以看到，IRLS收敛比梯度下降更快（因为是二阶方法）。

但图中没给出的是，IRLS（作为二阶方法）由于需要计算\(d\times d\)的Hessian矩阵（\(d\)为特征数量），计算成本比梯度下降高很多。

更具体的，IRLS每轮迭代时间复杂度为\(O(N+d^3)\)，\(N\)为样本数量，\(d\)为特征数量；而梯度下降（上升）的每轮迭代时间复杂度为\(O(Nd)\)。

下面，我们讨论一些对牛顿法（IRLS法）或梯度下降法的时间复杂度进行改进的算法。

拟牛顿法（Quasi-Newton Method）

https://www.osti.gov/biblio/4252678

牛顿法（或者说IRLS法）最耗时的操作就是计算\(d\times d\)的Hessian矩阵。

拟牛顿法（更具体的，BFGS算法）通过维护一个近似的Hessian矩阵（迭代过程中仅需要\(O(d^2)\)来更新这个近似矩阵）来替代Hessian矩阵。

共轭梯度法（Conjugate Gradient Method）

https://nvlpubs.nist.gov/nistpubs/jres/049/jresv49n6p409_a1b.pdf

经典的梯度下降算法有很多时间都浪费在了“走回头路”上，于是可以控制算法，让其每步都尽量是往新的方向走。

共轭梯度法通过结合当前梯度和上一轮方向，来计算更适合的（不走回头路的）梯度。

相比于牛顿法，共轭梯度法不需要维护\(d\times d\)的矩阵，仅存储几个向量。对于特征维度\(d\)特别大的情况更适合。

随机梯度下降法（Stochastic Gradient Descent, SGD）

由于现代分类任务的样本数据量过于庞大（\(N\)十分大），每次更新都考虑所有\(N\)个数据在硬件上十分困难。

为此，现在基本都用随机梯度下降法（及其变种）来迭代。

更具体的，考虑优化问题 \[ \min_{w} \frac{1}{N} \sum_{i=1}^{N} f(\mathbf{w}; \mathbf{x}_i, y_i) + \psi(\mathbf{w}) \] 随机梯度下降法（SGD）每次更新\(w\)仅考虑\(N\)个数据中随机的一小部分（random mini-batch） \[ \mathbf{g}_t = \frac{1}{|B_t|} \sum_{i \in B_t} \nabla_w f(\mathbf{w}; \mathbf{x}_i, y_i) + \nabla_w \psi(\mathbf{w}) \] 并用其更新 \[ \mathbf{w}_{t+1} = \mathbf{w}_t - \eta_t \mathbf{g}_t; \]

（来自 The many faces of integration by parts – II : Randomized smoothing and score functions – Machine Learning Research Blog）

对于特定的学习率\(\eta_t\)（Robbins-Monro条件），可以证明SGD最终一定收敛到最优： \[ \lim_{t \to \infty} \eta_t = 0, \quad \sum_{t \ge 1} \eta_t = +\infty, \quad \sum_{t \ge 1} \eta_t^2 < +\infty \]

• 对于凸函数，SGD一定收敛到全局最优。
• 对于非凸函数，SGD一定收敛到局部最优。

证明参考（证明基本依赖于Martingale或Lyapunov函数）：

A Stochastic Approximation Method

Handbook of Convergence Theorems for (Stochastic) Gradient Methods

正如上面提到的，现在数据量这么庞大，基本都是用随机梯度下降算法来优化，因此存在很多其不同的优化算法来优化经典SGD算法的两大基本缺陷：

1. SGD在“峡谷”型（一个方向很陡，另一个方向很平）函数表面会来回震荡，收敛很慢。
2. 所有参数共享一个学习率，且这个学习率需要手动精细调整（调参）。

下面列出几个SGD变种的特点。

算法	核心思想	学习率	动量	优点	缺点
SGD	沿梯度反方向更新	固定的，全局	否	简单，内存小	收敛慢，易震荡
Momentum	增加“惯性”（速度）	固定的，全局	是	加速收敛，减少震荡	可能会“冲过头”
NAG	“预判”的动量	固定的，全局	是 (更智能)	比 Momentum 抖动小，收敛快	略微复杂
Adagrad	历史梯度平方和	自适应，单调递减	否	适合稀疏数据	学习率会过早“死亡”
Rmsprop	历史梯度平方的均值	自适应，非单调	否	解决 Adagrad 的死亡问题	-
Adam	Momentum + Rmsprop	自适应，非单调	是 (自适应)	结合两者优点，稳健，常用	内存占用略大

具体公式可以参考朱军老师的《概率机器学习》、

可以看到，这些变种基本围绕两个核心思想：动量法（momentum）和自适应学习率（adaptive learning rate）。

动量法（Momentum）

通过引用类似于物理的“速度”概念来加速SGD。

动量法的基本思想是在SGD基础上增加一个动量（“速度”向量），其基于过去所有更新方向的指数衰减平均。

对于SGD第一个缺陷，动量可以抵消掉那些来回震荡的梯度。

不过，正如现实中存在惯性，动量可能会导致迭代过程中在（局部或者全局）最低点附近来回震荡。

为此，SGD变种NAG（Nesterov Accelerated Gradient）给动量加了某种“预判”来减少过冲（overshooting）问题。

更具体的，标准的动量是先计算当前梯度，再叠加上一步的动量。NAG则是先按照上一步的动量“预跳”一下，在那个“未来”的位置计算梯度，再用这个梯度来修正最终的更新方向。

自适应学习率（Adaptive learning rate）

通过为每一个参数自动调整其学习率，来对更需要学习的参数进行更多调整。

Adagrad（adaptive gradient）为每个参数分配一个单独的学习率，学习率动态的基于过去所有梯度平方和的平方根。于是，Adagrad会给那些不经常更新的参数一个较大的学习率。但问题是，梯度平方和随着迭代次数单调递增，这导致学习率会不断衰减，影响后续学习。

为此，Rmsprop（root mean square propagation）在Adagrad的基础上使用指数衰减移动平均来计算梯度平方和的平方根；这使得Rmsprop只关注近期的梯度大小。

---

目前最流行的变种是Adam（adaptive moment estimation），其结合了动量法和Rmsprop的自适应学习率。

下面是SGD及其几个变种对于同一个任务的迭代动画：

（来自[1609.04747] An overview of gradient descent optimization algorithms，An overview of gradient descent optimization algorithms）

朴素贝叶斯 v.s. 逻辑回归

在2001，Ng和Jordan就通过朴素贝叶斯和逻辑回归考察了生成式模型和判别式模型的一些区别：

On Discriminative vs. Generative Classifiers: A comparison of logistic regression and naive Bayes

概念上的不同

通过前面的关于同方差高斯朴素贝叶斯模型和逻辑回归的等价关系，可以看到他们的模型表达能力在特殊情况相等。

不过对于更一般的情况，因为他们考虑的分布不同，且优化的目标也不同，往往有一些差异。

数据无限多

假设数据无限多。如果朴素贝叶斯假设（特征条件独立）确实成立的话，那么基于上面论文，有 \[ \epsilon_{\text{Dis}, \infty}\sim\epsilon_{\text{Gen}, \infty} \] 其中Dis为判别discriminative的缩写，Gen为生成generative的缩写。

如果不成立的话， \[ \epsilon_{\text{Dis}, \infty}<\epsilon_{\text{Gen}, \infty} \] 跟直观上一样，朴素贝叶斯的表现就不如逻辑回归了。

数据有限

对于现实情况（尤其是数据量不多的时候），假设有\(N\)个数据点，\(d\)个特征： \[ \epsilon_{\text{Dis}, n} \leq \epsilon_{\text{Dis}, \infty} + O\left(\sqrt{\frac{d}{N}}\right) \\ \epsilon_{\text{Gen}, n} \leq \epsilon_{\text{Gen}, \infty} + O\left(\sqrt{\frac{\log d}{N}}\right) \] 朴素贝叶斯仅需要\(N=O(\log d)\)个数据来保证错误率收敛；而逻辑回归需要\(N=O(d)\)个数据。

根本原因在于，朴素贝叶斯假设导致了每个特征之间独立，并且参数之间非耦合（not coupled，不相互依赖）。所以朴素贝叶斯模型所需的样本数\(N\)只需要克服总参数的不确定性，而这种不确定性只随特征维度\(d\)对数级增长。

而逻辑回归是直接学习维度的联合估计（joint estimation），参数之间耦合（相互依赖）。为了准确稳定地估计一个\(d\)维的向量\(w\)，所需的样本量\(N\)必须线性地依赖于特征维度\(d\)。这也体现了维度的诅咒（Curse of Dimensionality）这一概念。

实验比较

在经典数据集上，可以看到，他们实际中各有千秋。

Revisiting Discriminative vs. Generative Classifiers: Theory and Implications

这篇论文比较了朴素贝叶斯和逻辑回归作用在深度学习提取的特征上时的表现（结果如上面两个图）。可以看到关于两者错误率收敛所需数据量的结果依旧成立。

指数族（Exponential Family）

指数族分布是统计学中一类重要的概率分布，它包含了许多常见的分布，如伯努利、多项、高斯、泊松、伽马等。

对于随机变量\(X\)，如果一个概率分布有一下形式，那么就称其为指数族分布： \[ \begin{aligned} p(\mathbf{x} \mid \mathbf{\eta}) &= h(\mathbf{x}) \exp \left(\mathbf{\eta}^\top T(\mathbf{x}) - A(\mathbf{\eta})\right) \\ &= \frac{1}{Z(\mathbf{\eta})} h(\mathbf{x}) \exp \left(\mathbf{\eta}^\top T(\mathbf{x})\right) \end{aligned} \] 其中：

• \(\eta\)是自然参数（natural/canonical parameter）。
• \(T(x)\)是充分统计量（sufficient statistic）。
• \(h(x)\)是基础函数。
• \(A(\eta)=\log Z(\eta)\)是对数归一化函数（log normalizer），用来保证概率综合为1。

可以从指数族的角度给不同分布的参数优化（MLE等）一个统一的框架。

指数族分布参数转换例子

多元高斯分布

考虑高斯分布 \[ \begin{aligned} p(\mathbf{x} \mid \mathbf{\mu}, \mathbf{\Sigma}) &= \frac{1}{(2\pi)^{d/2}|\mathbf{\Sigma}|^{1/2}} \exp \left(-\frac{1}{2}(\mathbf{x}-\mathbf{\mu})^\top \mathbf{\Sigma}^{-1}(\mathbf{x}-\mathbf{\mu})\right) \\ &= \frac{1}{(2\pi)^{d/2}|\mathbf{\Sigma}|^{1/2}} \exp \left(-\frac{1}{2}\text{tr}(\mathbf{\Sigma}^{-1}\mathbf{x}\mathbf{x}^\top) + \mathbf{\mu}^\top \mathbf{\Sigma}^{-1}\mathbf{x} - \frac{1}{2}\mathbf{\mu}^\top \mathbf{\Sigma}^{-1}\mathbf{\mu} - \frac{1}{2}\log|\mathbf{\Sigma}|\right) \end{aligned} \] 通过下面定义可以将其转换为指数族分布的形式： \[ \begin{align} \mathbf{\eta} &= \left[\mathbf{\Sigma}^{-1}\mathbf{\mu}; -\frac{1}{2}\text{vec}(\mathbf{\Sigma}^{-1})\right] = [\mathbf{\eta}_1; \text{vec}(\mathbf{\eta}_2)], \quad \mathbf{\eta}_1 = \mathbf{\Sigma}^{-1}\mathbf{\mu} \quad \mathbf{\eta}_2 = -\frac{1}{2}\mathbf{\Sigma}^{-1}\\ T(\mathbf{x}) &= \left[\mathbf{x}; \text{vec}(\mathbf{x}\mathbf{x}^\top)\right]\\ A(\mathbf{\eta}) &= \frac{1}{2}\mathbf{\mu}^\top \mathbf{\Sigma}^{-1}\mathbf{\mu} + \log|\mathbf{\Sigma}| = -\frac{1}{2}\text{tr}(\mathbf{\eta}_2\mathbf{\eta}_1\mathbf{\eta}_1^\top) - \frac{1}{2}\log(-2\mathbf{\eta}_2)\\ h(\mathbf{x}) &= (2\pi)^{-d/2} \end{align} \]

多项分布

考虑多项分布 \[ \begin{aligned} p(\mathbf{x} \mid \mathbf{\pi}) &= \prod_{i=1}^{d} \pi_i^{x_i} = \exp \left(\sum_{i} x_i \ln \pi_i\right) \\ &= \exp \left(\sum_{i=1}^{d-1} x_i \ln \pi_i + \left(1-\sum_{i=1}^{d-1} x_i\right) \ln \left(1-\sum_{i=1}^{d-1} \pi_i\right)\right) \\ &= \exp \left(\sum_{i=1}^{d-1} x_i \ln \frac{\pi_i}{1-\sum_{i=1}^{d-1} \pi_i} + \ln \left(1-\sum_{i=1}^{d-1} \pi_i\right)\right) \end{aligned} \] 下组定义转换其为指数族分布： \[ \begin{align} \mathbf{\eta} &= [\ln(\pi_i/\pi_d); 0]\\ T(\mathbf{x}) &= \mathbf{x}\\ A(\mathbf{\eta}) &= -\ln \left(1-\sum_{i=1}^{d-1} \pi_i\right)= \ln \left(\sum_{i=1}^{d} e^{\eta_i}\right)\\ h(\mathbf{x}) &= 1 \end{align} \]

矩生成（Moment generating）

指数族分布最强大的特性之一是其对数归一化函数\(A(\eta)\)与分布的矩（如均值和方差）之间存在直接关系：

• 一阶导数（均值）：对 \(A(\eta)\) 求一阶导数，得到充分统计量 \(T(x)\) 的期望（即矩参数 \(\mu\)）。 \[ \nabla_{\eta}A(\eta) = \mathbb{E}_{p(x|\eta)}[T(x)] \triangleq \mu \]

• 二阶导数（方差）：对 \(A(\eta)\) 求二阶导数，得到 \(T(x)\) 的方差。 \[ \nabla_{\eta}^{2}A(\eta) = Var[T(x)]>0 \]

并且，由于方差在连续分布通常大于0（不然就退化成了单点分布，或成为狄拉克\(\delta\)函数），\(A(\eta)\)二阶导严格大于0，所以\(A(\eta)\)是严格凸函数。这也表明了\(\mu=\nabla_\eta A(\eta)\)作为\(A(\eta)\)的一阶导严格单调递增，于是我们就有了自然参数\(\eta\)和矩参数\(\mu\)之间的一对一关系： \[ \eta \triangleq \psi(\mu) \] 上面的\(\psi\)是一个一对一映射（保证了参数之间转换的可逆性）。

独立同分布的联合分布也是指数族

\[ \begin{aligned} p(D \mid \mathbf{\eta}) &= \prod_{i} h(\mathbf{x}_i) \exp \left(\mathbf{\eta}^\top T(\mathbf{x}_i) - A(\mathbf{\eta})\right) \\ &= \left(\prod_{i} h(\mathbf{x}_i)\right) \exp \left(\mathbf{\eta}^\top \sum_{i} T(\mathbf{x}_i) - N A(\mathbf{\eta})\right) \end{aligned} \]

并且，充分统计量\(T\)为所有样本的充分统计量总和： \[ \sum_{i} T(\mathbf{x}_i) \]

指数族分布的最大似然估计MLE对应了统计学的矩匹配（Moment matching）

矩匹配是一种参数估计方法 。它的核心思想是：

• 理论矩（Theoretical Moment）：即分布的统计矩（如均值、方差等，通常是关于模型参数的函数）。

• 经验矩（Empirical Moment）：即从样本数据中计算出来的统计矩（如样本均值、样本方差）。

矩匹配的过程就是设置理论矩等于经验矩，然后解方程组来求得最优参数的估计值。

对于独立同分布数据，考虑其指数族形式的对数似然 \[ \mathcal{L}(\mathbf{\eta}; D) = \sum_{n} \log h(\mathbf{x}_n) + \left(\mathbf{\eta}^\top \sum_{n} T(\mathbf{x}_n)\right) - N A(\mathbf{\eta}) \] 取对于\(\eta\)的梯度，设为0，可以得到最优矩参数 \[ \begin{align} \nabla_{\mathbf{\eta}} \mathcal{L}(\mathbf{\eta}; D) &= \sum_{n} T(\mathbf{x}_n) - N \nabla_{\mathbf{\eta}} A(\mathbf{\eta}) = 0\\ \nabla_{\mathbf{\eta}} A(\mathbf{\eta}) &= \frac{1}{N} \sum_{n} T(\mathbf{x}_n)\\ \hat{\mathbf{\mu}}_{\text{MLE}} &= \frac{1}{N} \sum_{n} T(\mathbf{x}_n) \end{align} \] 可以看到最优矩参数\(\hat{\mathbf{\mu}}_{\text{MLE}}\)（仅需要充分统计量\(T\)来计算）刚好对应了样本数据的最优经验矩。所以，MLE就是矩匹配。

从这个角度来看，对于指数族分布，用矩匹配的方式，先根据样本充分统计量\(T\)找到最优矩参数\(\hat{\mathbf{\mu}}_{\text{MLE}}\)，然后（用一对一函数\(\psi\)）匹配到自由自然参数\(\eta\)。这么计算避免了对许多参数进行复杂的迭代优化。

最优矩参数\(\hat{\mathbf{\mu}}_{\text{MLE}}\)的例子

高斯分布

对于高斯分布，充分统计量\(T(x)\)是一个向量，包含了\(x\)和 \(xx^\top\)的向量化形式\([x;vec(xx^{\top})]\)。 \[ \begin{align*} \mathbf{\eta} &= \left[\mathbf{\Sigma}^{-1}\mathbf{\mu}; -\frac{1}{2}\text{vec}(\mathbf{\Sigma}^{-1})\right]\\ T(\mathbf{x}) &= \left[\mathbf{x}; \text{vec}(\mathbf{x}\mathbf{x}^\top)\right] \\ A(\mathbf{\eta}) &= \frac{1}{2}\mathbf{\mu}^\top \mathbf{\Sigma}^{-1}\mathbf{\mu} + \log|\mathbf{\Sigma}| \\ h(\mathbf{x}) &= (2\pi)^{-d/2}\\ &\Rightarrow \hat{\mathbf{\mu}}_{\text{MLE}}= \frac{1}{N} \sum_{n} T_1(\mathbf{x}_n) = \frac{1}{N} \sum_{n} \mathbf{x}_n \end{align*} \]

多项分布

对于多项分布和下面的泊松分布，充分统计量\(T(x)\)就是特征\(x\)本身。 \[ \begin{align*} \mathbf{\eta} &= [\ln(\pi_i/\pi_d); 0] \\ T(\mathbf{x}) &= \mathbf{x} \\ A(\mathbf{\eta}) &= -\ln \left(1-\sum_{i=1}^{d-1} \pi_i\right) \\ h(\mathbf{x}) &= 1\\ &\Rightarrow \hat{\mathbf{\mu}}_{\text{MLE}} = \frac{1}{N} \sum_{n} \mathbf{x}_n \end{align*} \]

泊松分布

\[ \begin{align*} \eta &= \log \lambda \\ T(x) &= x \\ A(\eta) &= \lambda = e^\eta \\ h(x) &= \frac{1}{x!}\\ \Rightarrow \hat{\mu}_{\text{MLE}} &= \frac{1}{N} \sum_{n} x_n \end{align*} \]

广义线性模型（Generalized Linear Models, GLIMs）

上一节的指数族尝试将不同分布（的参数优化）囊括在同一个框架下；而这一节的广义线性模型如之前提到的，基于指数族分布的框架，尝试将线性回归、逻辑回归等模型囊括在同一个框架下。

基本定义

广义线性模型包含了三个核心假设：

响应变量的条件分布属于指数族

给定输入\(X\)，输出（相应变量）\(y\)的条件概率分布\(P(y|X)\)假设为指数族分布。

这个假设决定了要解决的问题类型。

• 如果\(y\)是连续值（房价预测），选择高斯分布（对应线性回归）。
• 如果\(y\)是二元值（是/否），选择伯努利分布（对应逻辑回归）。
• 如果\(y\)是计数（网页点击量），选择泊松分布。

预测器为线性

模型的输入 \(X\) 假设是通过一个线性组合来影响输出： \[ \xi = \theta^{\top}X \] 这里的\(\xi:=\theta^\top X\)被称为线性预测器（linear predictor）。

这个假设也是名称广义线性模型中线性的由来。

链接函数的存在

假设存在响应函数\(f\)链接线性预测器\(\xi\)和响应变量\(y\)的期望\(\mu\)： \[ \mathbb{E}_p[y]=\mu = f(\xi) = f(\theta^{\top}X) \] 响应函数\(f\)确保了线性预测器\(\xi\)能够映射输入\(X\)到正确的\(\mu\)值。

• 对于逻辑回归，响应函数\(f(\xi) = \frac{1}{1+e^{-\xi}}\)把输入的线性组合\(\xi=\theta^\top X=w^\top X+w_0\)映射到区间\((0,1)\)。

模型	分布	线性预测器 ξ	自然参数 η	矩参数 μ (E[y])	响应函数 f(ξ)
线性回归	高斯分布	\(\theta^{\top}X\)	\(\eta = \mu\)	\(\mu\)	\(f(\xi) = \xi\)（恒等函数）
逻辑回归	伯努利分布	\(\theta^{\top}X\)	\(\eta = \log\left(\frac{\mu}{1-\mu}\right)\)	\(\mu\)	\(f(\xi) = \frac{1}{1+e^{-\xi}}\)（Sigmoid）

标准响应函数

我们进一步增加对于链接函数的假设；一会可以看到限制下的广义线性模型的参数求解都极其简洁（和指数族求解最优矩参数\(\mu\)一样）。

考虑上一节中矩生成定义的一对一映射函数\(\psi\)： \[ \eta \triangleq \psi(\mu) \] 现在，假设响应函数\(f\)为标准响应函数\(f=\psi^{-1}\)。

从上图可以看到，自然参数\(\eta=\psi(\mu)=\psi(f(\xi))\)，代入假设得到\(\eta=\xi\)。

那就是，响应函数假设等价于假设线性预测器\(\xi\)等于自然参数\(\eta\)。

当然，最基本的线性回归和逻辑回归的响应函数都满足此假设：

• 对于线性回归，\(\eta\)恰好等于\(\mu\)，所以\(\eta=\xi=\theta^\top X\)，满足标准响应函数假设。
• 对于逻辑回归，\(\eta=\log\left(\frac{\mu}{1-\mu}\right)\)直接等于线性预测器\(\xi:=\theta^\top X\)。

广义线性模型的最大似然估计MLE

考虑指数组分布为：
\[ \begin{gather} p(x \mid \eta) = h(x) \exp \left(\eta^\top T(x) - A(\eta)\right) \end{gather} \]

对数似然为
\[ \begin{aligned} \mathcal{L}(\mathbf{\theta}; D) &= \sum_{n} \log h(y_n) + \sum_{n} (\eta_n y_n - A(\eta_n)) \\ \eta_n &= \psi(\mu_n), \quad \mu_n = f(\xi_n), \quad \xi_n = \mathbf{\theta}^\top \mathbf{x}_n \end{aligned} \] 导数为 \[ \begin{aligned} \nabla_{\mathbf{\theta}} \mathcal{L} &= \sum_{n} \left(y_n \nabla_{\mathbf{\theta}}\eta_n - \frac{d A(\eta_n)}{d\eta_n} \nabla_{\mathbf{\theta}}\eta_n\right)\\ &= \sum_{n} (y_n - \mu_n) \nabla_{\mathbf{\theta}}\eta_n \end{aligned} \] 可以看到求解让梯度为0的\(\theta\)并不简单。

标准响应函数假设下的广义线性模型的最大似然估计MLE

但是如果（根据标准响应函数假设）\(\eta_n=\theta^\top x_n\)，那么\(\nabla_\theta \eta_n=x_n\)大大简化对数似然的导数： \[ \begin{align} \nabla_{\theta}\mathcal{L} &= \sum_{n}(y_{n}-\mu_{n})x_{n}\\ &= X(y-\mu) \end{align} \] 这完完全全就是上文逻辑回归的优化算法小节需要优化的形式。

于是，可以用各种各样的算法，如牛顿法（也是IRLS法）、随机梯度下降（上升）等方法来在线学习并优化最优模型参数\(\theta\)； \[ \mathbf{\theta}_{t+1} = \mathbf{\theta}_t + \rho(y_n - \mu_n^t)\mathbf{x}_n\\ \mu_n^t = f(\mathbf{\theta}_t^\top \mathbf{x}_n) \]

---