ESL4.3 线性判别分析（LDA&QDA）学习笔记

4.3 线性判别分析

这是一篇有关《统计学习基础》，原书名The Elements of Statistical Learning的学习笔记，该书学习难度较高，有很棒的学者将其翻译成中文并放在自己的个人网站上，翻译质量非常高，本博客中有关翻译的内容都是出自该学者的网页，个人解读部分才是自己经过查阅资料和其他学者的学习笔记，结合个人理解总结成的原创内容。
有关ESL更多的学习笔记的markdown文件，可在作者GitHub上查看下载。

原文	The Elements of Statistical Learning
翻译	szcf-weiya
时间	2018-08-21
注解	Hytn Chen
更新	2020-02-23

若想深入理解本节内容，首先需要对协方差矩阵有更深更直观的理解，阅读这篇总结后就明白了协方差矩阵为什么要进行特征值分解，其英文原文也可进行阅读。

翻译原文

分类的判别理论（2.4节）告诉我们，我们需要知道的是最优分类的类别后验概率 $\Pr(G\mid X)$Pr(G∣X)．假设 $f_k(x)$fk​(x) 是类别 $G=k$G=k 中 $X$X 的类别条件密度，并令 $\pi_k$πk​ 为类别 $k$k 的先验概率，满足 $\sum_{k=1}^K\pi_k=1$∑k=1K​πk​=1．简单地应用一下贝叶斯定理得到

$\Pr(G=k\mid X=x)=\dfrac{f_k(x)\pi_k}{\sum_{\ell=1}^Kf_{\ell}(x)\pi_\ell}\tag{4.7}$Pr(G=k∣X=x)=∑ℓ=1K​fℓ​(x)πℓ​fk​(x)πk​​(4.7)

我们看到就 判别能力(ability to classify) 而言，知道 $f_k(x)$fk​(x) 几乎等价于知道概率 $\Pr(G=k\mid X=x)$Pr(G=k∣X=x)．

---

Hytn注：

这里贝叶斯定理的应用过程如下，证明之前首先明确概念类别条件密度，即假定x是一个连续随机变量，其分布取决于类别状态，可表示成$P(X \mid G)$P(X∣G)的形式，也就是类别状态为G时输入特征X的概率密度函数。在上式中类条件概率密度函数就是$f_k(x)$fk​(x) ，类别$k$k先验概率$P(G)$P(G)就是$\pi_k$πk​。
$\begin{aligned} \Pr(G=k\mid X=x) &= \frac{P(GX)}{P(X)}\\ &=\frac{P(X \mid G)P(G)}{P(X)}\\ &=\dfrac{f_k(x)\pi_k}{\sum_{\ell=1}^Kf_{\ell}(x)\pi_\ell} \end{aligned}$Pr(G=k∣X=x)​=P(X)P(GX)​=P(X)P(X∣G)P(G)​=∑ℓ=1K​fℓ​(x)πℓ​fk​(x)πk​​​

---

许多方法是基于类别密度的模型：

• 采用高斯密度的线性和二次判别分析
• 更加灵活的混合的高斯密度，允许非线性判别边界（6.8 节）
• 对每个类别密度进行一般的非参数密度估计，允许最大的灵活性（6.6.2 节）
• 朴素贝叶斯(Naive Bayes) 模型是上个情形的变种，并且假设每个类密度是边缘密度的乘积，这也就是，假设输入在每一类中都是条件独立的（6.6.3节）

假设我们每个类别密度用多元高斯分布来建模

$f_k(x)=\frac{1}{(2\pi)^{p/2}\vert \boldsymbol\Sigma_k\vert^{1/2}}e^{-\frac{1}{2}(x-\mu_k)^T\boldsymbol\Sigma^{-1}_k(x-\mu_k)} \tag{4.8}$fk​(x)=(2π)p/2∣Σk​∣1/21​e−21​(x−μk​)TΣk−1​(x−μk​)(4.8)

假设类别有相同的协方差矩阵 $\Sigma_k=\Sigma\;\forall k$Σk​=Σ∀k 时，会导出 线性判别分析 (LDA)．在比较两个类别 $k$k 和 $\ell$ℓ 时，比较 log-ratio 就足够了，而且我们可以看到
$\begin{aligned} &\log\frac{\Pr(G=k\mid X=x)}{\Pr(G=\ell\mid X=x)}\\=&\log\frac{f_k(x)}{f_\ell(x)}+\log\frac{\pi_k}{\pi_\ell}\\ =&\log\frac{\pi_k}{\pi_\ell}-\frac{1}{2}(\mu_k+\mu_\ell)^T\boldsymbol\Sigma^{-1}(\mu_k-\mu_\ell)+x^T\boldsymbol\Sigma^{-1}(\mu_k-\mu_\ell) \end{aligned}\tag{4.9}$==​logPr(G=ℓ∣X=x)Pr(G=k∣X=x)​logfℓ​(x)fk​(x)​+logπℓ​πk​​logπℓ​πk​​−21​(μk​+μℓ​)TΣ−1(μk​−μℓ​)+xTΣ−1(μk​−μℓ​)​(4.9)

!!! note “weiya 注”
如果 $D$D 为对称矩阵，则
$a'Da-b'Db=(a+b)'D(a-b)\,.$a′Da−b′Db=(a+b)′D(a−b).
这是个关于 $x$x 的线性等式．协方差矩阵相等时会消除了正规化因子，以及指数中的二次项．这个线性的 log-odds 函数表明类别 $k$k 和类别 $\ell$ℓ 的判别边界为 $p$p 维超平面——$\Pr(G=k\mid X=x)=\Pr(G=\ell\mid X=x)$Pr(G=k∣X=x)=Pr(G=ℓ∣X=x) 的集合，这关于 $x$x 是线性的．对于任意类别对也是成立的，所以所有的判别边界都是线性的．如果我们把 $\rm{IR}^p$IRp 分成不同区域，记为类别 1，类别 2 等等，这些区域会被超平面所分离．图 4.5（左图）显示了有三个类别且 $p=2$p=2 的假想的例子．这里数据是从 3 个协方差矩阵相等的高斯分布中产生的．我们已经在图中画出了对应 95% 概率密度的等高图，以及三个类别的形心．注意到判别边界不是连接两个形心的垂直平分线．如果协方差矩阵 $\Sigma$Σ 是 $\sigma^2\mathbf I$σ2I，且类别的先验概率相等．从(4.9)我们可以看出线性判别函数为

$\delta_k(x)=x^T\boldsymbol\Sigma^{-1}\mu_k-\frac{1}{2}\mu_k^T\boldsymbol\Sigma^{-1}\mu_k+\log\pi_k\tag{4.10}$δk​(x)=xTΣ−1μk​−21​μkT​Σ−1μk​+logπk​(4.10)

这是判别规则的等价描述，$G(x)=\mathrm{argmax}_k \delta_k(x)$G(x)=argmaxk​δk​(x)．

图 4.5. 左图显示了三个高斯分布，有相同的协方差和不同的均值．图中画出了在包含每个类别 95% 可能性的等高线．两两类别之间的贝叶斯判别边界用虚线显示，并且贝叶斯判别边界将所有的三个类别分隔开用实线表示出来（前者的一个子集）．右图我们看到有来自每个高斯分布的 30 个样本点以及画出了拟合后的 LDA 判别边界．

实际应用中我们不知道高斯分布的参数，我们需要用我们的训练数据去估计它们：

• $\hat \pi_k=N_k/N$π^k​=Nk​/N，其中 $N_k$Nk​ 是第 $k$k 类观测值的个数；
• $\hat\mu_k=\sum_{g_i=k}x_i/N_k$μ^​k​=∑gi​=k​xi​/Nk​；
• $\hat{\boldsymbol\Sigma}=\sum_{k=1}^K\sum_{g_i=k}(x_i-\hat\mu_k)(x_i-\hat\mu_k)^T/(N-K)$Σ^=∑k=1K​∑gi​=k​(xi​−μ^​k​)(xi​−μ^​k​)T/(N−K)．

---

Hytn注：

观察式4.9，可以发现结果是关于x的一次函数$wx+w_0$wx+w0​，所以线性判别的说法由此得来。其中参数的计算就如上面训练数据估计式中写的那样。

关于LDA总结：LDA针对的就是数据服从高斯分布，且均值不同，方差相同的情况。得出式4.10的前提是假设了协方差矩阵 $\Sigma$Σ 是 $\sigma^2\mathbf I$σ2I，且类别的先验概率相等。

---

图 4.5（右图）显示了基于从 3 个高斯分布中取出的大小为 30 的样本集得到的判别边界的估计，图 4.1 是另外一个例子，但是那里类别不是高斯分布．
$\underset{B}{\min}\sum\limits_{i=1}^N\Vert y_i-[(1,x_i^T)\mathbf B]^T\Vert^2\tag{4.5}$Bmin​i=1∑N​∥yi​−[(1,xiT​)B]T∥2(4.5)
两个类别的情况下在线性判别分析和线性最小二乘之间有一个简单的对应，如(4.5)所示．如果满足下面条件则 LDA 分类规则分给第二类

$x^T\hat{\boldsymbol\Sigma}^{-1}(\hat\mu_2-\hat\mu_1)>\frac{1}{2}\hat\mu^T_2\hat\Sigma^{-1}\hat\mu_2-\frac{1}{2}\hat\mu_1^T\hat{\boldsymbol\Sigma}^{-1}\hat\mu_1+\log(N_1/N)-\log(N_2/N)\tag{4.11}$xTΣ^−1(μ^​2​−μ^​1​)>21​μ^​2T​Σ^−1μ^​2​−21​μ^​1T​Σ^−1μ^​1​+log(N1​/N)−log(N2​/N)(4.11)

$x^T\hat{\boldsymbol\Sigma}^{-1}(\hat\mu_2-\hat\mu_1)>\frac{1}{2}(\hat\mu_2+\hat\mu_1)^T\hat{\boldsymbol\Sigma}^{-1}(\hat\mu_2-\hat\mu_1)-\log(N_2/N_1)\tag{4.11}$xTΣ^−1(μ^​2​−μ^​1​)>21​(μ^​2​+μ^​1​)TΣ^−1(μ^​2​−μ^​1​)−log(N2​/N1​)(4.11)

否则分为第一类．假设我们将两个类别分别编码为 $+1$+1 和 $-1$−1．可以简单地证明得出最小二乘的系数向量与 4.11式给出的 LDA 方向成比例（练习 4.2）．【实际上，这种对应关系对于任意的编码都会有，见练习 4.2】．然而除了 $N_1=N_2$N1​=N2​ 的情形，截距不同因此得到的判别边界不一样．

---

Hytn注：

式(4.11)的条件等价于$\delta_{2}(x)>\delta_{1}(x)$δ2​(x)>δ1​(x)，如果满足这个那么就选择类别2，否则就选择类别1。其中贝叶斯判别函数表达式为
$\delta_{k}(x)=\ln \left(p\left(x | \omega_{k}\right)\right)+\ln \left(\pi_{k}\right)$δk​(x)=ln(p(x∣ωk​))+ln(πk​)

那么如何通过最小二乘得到LDA方向？可参考ESL4.2节的式(4.5)的批注，求解方式还是令均方误差最小化的求导，得到参数。

总结一下上面的线性判别函数就是描绘了各个特征属于多元高斯分布下时，两个类别之间边界的函数。高斯分布的参数由训练集估计得到。

---

因为通过最小二乘得到的 LDA 方向不需要对特征做高斯分布的假设，它的应用可以不局限于高斯分布的数据．然而，(4.11)式给出的特定的截距和分离点确实需要高斯分布的数据．因此对于给定的数据集凭经验地选择分割点最小化训练误差是有意义的．这也是我们在实际中发现效果很好，但是在理论中没有提到的原因．

多于两个类别的情形时，LDA 与类别指示矩阵的线性回归不是一样的，而且它避免了跟这个方法有关的掩藏问题 (Hastie et al., 1994¹)． 回归与 LDA 之间的对应可以通过在 12.5 节中讨论的最优评分来建立．

回到一般的判别问题(4.8)，如果没有假设 $\boldsymbol\Sigma_k$Σk​ 相等，则(4.9)中的抵消不会发生；特别地，关于 $x$x 的平方项保留了下来．于是我们得到了平方判别函数 QDA
$\delta_k(x)=-\frac{1}{2}\log\vert\Sigma_k\vert-\frac{1}{2}(x-\mu_k)^T\Sigma_k^{-1}(x-\mu_k)+\log\pi_k\tag{4.12}$δk​(x)=−21​log∣Σk​∣−21​(x−μk​)TΣk−1​(x−μk​)+logπk​(4.12)
每个 类别对 (pair of classes) $k$k 和 $\ell$ℓ 的判别边界由二次等式来描述 ${x:\delta_k(x)=\delta_\ell(x)\\}$x:δk​(x)=δℓ​(x)．

---

Hytn注：上式去掉第一项就变成了线性判别函数，而${x:\delta_k(x)=\delta_\ell(x)}$x:δk​(x)=δℓ​(x)判别边界是二次的。

QDA与LDA的区别是什么呢，当数据均值不同方差相同（LDA）的时候，面对两个数据的交汇交叠部分，直接进行线性的划分造成的误差不大，而如果交汇部分的数据对于各类的意义不一样（方差不同，QDA）直接线性划分就出问题了。有关该点描述图出自这篇博客：

---

图 4.6 显示了一个例子（图 4.1），这三个类别都是混合高斯分布（6.8 节），并且判别边界由关于 $x$x 的二次函数近似给出．这里我们描述两种拟合这些判别边界的方式．右图使用这里描述的 QDA，而左图是在增广的五维二次多项式空间中使用 LDA．两者之间的差别很小，QDA 是最好的方式，LDA 方法是更方便的替换．

图 4.6. 拟合二次边界的两种方法．左图显示了图 4.1 的二次判别边界（在五维空间 $X_1,X_2,X_1X_2,X_1^2,X_2^2$X1​,X2​,X1​X2​,X12​,X22​ 中运用 LDA 得到）．右图显示了通过 QDA 寻找到的二次判别边界．两者差别很小，通常也是这种情况．

QDA 的估计类似 LDA 的估计，除了协方差矩阵必须要按每一类来估计．当 $p$p 很大这意味着系数有显著性的增长．因为判别边界是系数密度的函数，参数的个数必须要考虑．对于 LDA，似乎有 $(K-1)\times(p+1)$(K−1)×(p+1) 个参数，因为我们仅仅需要判别函数之间的不同 $\delta_k(x)-\delta_K(x)$δk​(x)−δK​(x)，其中 $K$K 是一些预先选好的类别（这里我们已经选了最后一类），每个差异需要 $p+1$p+1 个参数．对于 QDA 类似地，我们会有 $(K-1)\times{p(p+3)/2+1\\}$(K−1)×p(p+3)/2+1 个参数．LDA 和 QDA 在非常大以及离散的数据集的分类上面表现得很好．举个例子，在 STATLOG 项目中 (Michie et al. 1994²) LDA 在 7 个数据集（总共 22 个数据集）中的表现排前三名，QDA 在四个数据集中排前三名，对于 10 个数据集两种方法的其中一种排前三名．两种方法都被广泛运用，整本书集中讨论 LDA．在各种外来方法风靡一时的今天，我们总是会有两种简单的方法可以使用．为什么 LDA 和 QDA 有那么好的效果？原因不可能是数据近似服从高斯分布，对于 LDA 协方差矩阵也不可能近似相等．很可能的一个原因是数据仅仅可以支持简单的判别边界比如线性和二次，并且通过高斯模型给出的估计是稳定的，这是一个偏差与方差之间的权衡——我们可以忍受线性判别边界的偏差因为它可以通过用比其它方法更低的方差来弥补．这个论点对于 QDA 更是不可想象，因为它自身有许多的参数，尽管或许比非参估计的参数要少．

正则化判别分析

Friedman (1989)³ 提出 LDA 和 QDA 之间的一个权衡，使得 QDA 的协方差阵向 LDA 中的共同协方差阵收缩．这些方法非常类似岭回归．正则化协方差矩阵有如下形式
$\hat{\boldsymbol\Sigma}_k(\alpha)=\alpha\hat{\boldsymbol\Sigma}_k+(1-\alpha)\hat{\boldsymbol\Sigma}\tag{4.13}$Σ^k​(α)=αΣ^k​+(1−α)Σ^(4.13)
其中，$\hat\Sigma$Σ^ 是和 LDA 一样用的联合协方差矩阵．这里 $\alpha\in[0,1]$α∈[0,1] 允许在 LDA 和 QDA 之间连续变化的模型，而且需要指定．实际中，$\alpha$α 可以基于在验证数据的表现上进行选择，或者通过交叉验证．

图 4.7. 对元音数据应用一系列 $\alpha\in[0,1]$α∈[0,1] 的正则化判别分析的测试和训练误差．测试数据的最优点发生在 $\alpha=0.9$α=0.9 附近，离二次判别分析很相近．

图 4.7 显示了将 RDA 运用到元音数据的结果．训练和测试误差都随着 $\alpha$α 的增大获得了改善，尽管在 $\alpha=0.9$α=0.9 后测试误差急剧上升．训练和测试误差最大的区别部分因为在小数量的个体上有很多重复观测，在训练和测试集上是不同的．

类似的修改使得 $\hat\Sigma$Σ^ 向着标量协方差收缩，对于 $\gamma\in[0,1]$γ∈[0,1]，有
$\hat{\boldsymbol\Sigma}(\gamma)=\gamma\hat{\boldsymbol\Sigma}+(1-\gamma)\hat{\sigma}^2\mathbf I.\tag{4.14}$Σ^(γ)=γΣ^+(1−γ)σ^2I.(4.14)
用 $\hat{\boldsymbol\Sigma}(\gamma)$Σ^(γ) 替换掉 (4.13) 的 $\hat{\boldsymbol\Sigma}$Σ^ 导出了一个更加一般的协方差阵族 $\hat{\boldsymbol\Sigma}(\alpha,\gamma)$Σ^(α,γ)，这是由一对参数来表示的．

在 12 章，我们将讨论 LDA 的另一种正则化版本，当数据来自数字化的相似的信号和图像时，这种方法会更加地合适．在这些情形下特征是高维的且相关的，LDA 系数可以通过正则化在原始信号域变得光滑和稀疏．这导出了更好的一般化并且对于这些系数有着简单的解释．在第 18 章中我们也处理非常高维的问题，举个例子如微阵列研究中基因表达测量的特征．那里方法集中在(4.14)中 $\gamma=0$γ=0 的情形，以及其它 LDA 的正则化版本．

LDA 的计算

作为下一主题的铺垫，我们简单地岔开去讨论 LDA，特别是 QDA 的计算．这些计算可以通过对角化 $\hat\sigma$σ^ 或 $\hat\sigma_k$σ^k​ 来简化．对于后者，假设我们对每一个计算特征值分解 $\hat\sigma_k=\mathbf U_k\mathbf D_k\mathbf U_k^T$σ^k​=Uk​Dk​UkT​，其中 $\mathbf U_k$Uk​ 是 $p\times p$p×p 的正交矩阵，$\mathbf D_k$Dk​ 是正的特征值 $d_{k\ell}$dkℓ​ 组成的对角矩阵．则 $\delta_k(x)$δk​(x) (4.12) 的组成成分是

• $(x-\hat\mu_k)^T\hat\sigma_k^{-1}(x-\hat\mu_k)=[U_k^T(x-\hat\mu_k)]^TD_k^{-1}[U_k^T(x-\hat\mu_k)]$(x−μ^​k​)Tσ^k−1​(x−μ^​k​)=[UkT​(x−μ^​k​)]TDk−1​[UkT​(x−μ^​k​)]
• $\log\vert \hat\sigma_k\vert=\sum_{\ell}\log d_{k\ell}$log∣σ^k​∣=∑ℓ​logdkℓ​

按照上面列出的计算步骤，LDA 分类器可以通过下面的步骤来实现

• 对数据关于协方差矩阵 $\hat{\boldsymbol\Sigma}$Σ^ 球面化：$X^*\leftarrow \mathbf D^{-\frac{1}{2}}\mathbf U^T\mathbf X$X∗←D−21​UTX，其中 $\hat{\boldsymbol\Sigma}=\mathbf U\mathbf D\mathbf U^T$Σ^=UDUT．$X^*$X∗ 的共同协方差矩阵变为单位阵．
• 考虑类别先验概率 $\pi_k$πk​ 的影响，在变换后的空间里面分到最近的类别形心．

---

Hytn注：

这里的$\hat\sigma_k$σ^k​可理解为类内协方差矩阵，是以一个类别内所有的训练样本为数据计算期望和方差后，再计算得出特征之间的协方差矩阵的。对这个协方差矩阵进行特征分解可得到最重要的几个特征（根据特征值大小衡量）的顺位排序。

$X^*$X∗ 为伴随矩阵，伴随矩阵的求法不再赘述，其性质如下
$A A^{*}=A^{*} A=|A| E_{n}$AA∗=A∗A=∣A∣En​
共同协方差矩阵其实就是组内协方差矩阵，因为是线性判别分析，方差相同，所以就变成了每个类的协方差矩阵都一致，也就是前文所述的$\Sigma$Σ 是 $\sigma^2\mathbf I$σ2I。

---

降维线性判别分析

至此我们已经讨论了限制为高斯分类器的 LDA．它受欢迎的部分原因是因为额外的限制使得我们可以看到数据在低维空间中富含信息的投影．

在 $p$p 维输入空间的 $K$K 个形心位于维数 $\le K-1$≤K−1 的超平面子空间中，如果 $p$p 比 $K$K 大很多，维数上会有显著的降低．更多地，在确定最近的形心时，我们可以忽略到子空间的垂直距离，因为它们对每个类的作用同样大．因此我们可能仅仅需要将数据 $X^*$X∗ 投射到形心张成的子空间 $H_{K-1}$HK−1​ 中，并且在子空间内比较距离．因此在 LDA 中存在一个基础维数的降低，换句话说就是，我们仅仅需要考虑在维数至多为 $K-1$K−1 的子空间的数据．如果 $K=3$K=3，举个例子，这允许我们可以在二维图中观察数据，对类别进行颜色编码．这样做我们不会丢失任何 LDA 分类需要的信息．

如果 $K > 3$K>3？我们可能在某种意义下要求一个最优的 $L < K-1$L<K−1 维子空间 $H_L\subseteq H_{K-1}$HL​⊆HK−1​．Fisher 定义的最优意思是投影形心关于方差要尽可能地分散．这意味着寻找形心的主成分空间（主成分在 3.5.1 节中有简短的描述，在 14.5.1 节中将详细讨论）．图 4.4 显示了对于元音数据的一个最优的二维子空间．这里是一个 $10$10 维的输入空间， 其中有 $11$11 个类别，每一类指不同的元音发音．形心在这种情形下要求全空间，因为 $K-1=p$K−1=p，但是我们已经展示了一个最优的二维子空间．维数是有序的，所以我们可以依次计算额外的维数．图 4.8 显示了 $4$4 个额外的坐标对，也被称作 典则 (canonical) 或者 判别 (discriminant) 变量．

图 4.8. 在成对典则变量上的四个投影．注意到当典则变量的秩增大，形心变得更不发散．右下角的图像看起来像是叠加上去的，类别也更加难以确定．

总结一下，寻找 LDA 的最优子空间序列涉及以下步骤：

• 计算 $K\times p$K×p 的类别形心矩阵 $\mathbf M$M 以及共同协方差矩阵 $\mathbf W$W（组内 (within-class) 协方差）
• 使用 $\mathbf W$W 特征值分解计算 $\mathbf M^*=\mathbf M\mathbf W^{-\frac{1}{2}}$M∗=MW−21​
• 计算 $\mathbf M^*$M∗ 的协方差矩阵 $\mathbf B^*$B∗，（$\mathbf B$B 是 组间 (between-class) 协方差），以及特征值分解 $\mathbf B^*=\mathbf V^*\mathbf D_B{\mathbf V^*}^T$B∗=V∗DB​V∗T．$\mathbf V^*$V∗ 的列 $v_\ell^*$vℓ∗​ 从第一个到最后一个依次定义了最优子空间的坐标．

结合上述的操作，第 $\ell$ℓ 个 判别变量 (discriminant variable) 由 $Z_\ell=v_\ell^TX$Zℓ​=vℓT​X 给出，其中 $v_\ell=W^{-\frac{1}{2}}v_\ell^*$vℓ​=W−21​vℓ∗​．

---

Hytn注：

首先解释特征值分解，该篇博客将特征值与特征向量的概念与应用实例解释得非常完美，可以将之前的许多点串在一起，详尽阅读非常有助理解。

有关矩阵的负二分之一次方的计算步骤：假如$A=Q^{-1}\Lambda Q$A=Q−1ΛQ，则
$A^{\frac{1}{2}}=Q^{-1}(\Lambda)^{\frac{1}{2}} Q$A21​=Q−1(Λ)21​Q
之后再求$(A^{\frac{1}{2}})^{-1}$(A21​)−1即可。

M中存储了K个类别的形心坐标，每个坐标是p维的。W代表了每个类内的协方差，

---

!!! note “weiya 注：”
结合 Ex. 4.1 的证明过程来理解上述算法．主要思想是先求 $\mathbf W^{-\frac 12}\mathbf B\mathbf W^{-\frac 12}$W−21​BW−21​ 的特征向量 $v_\ell^*$vℓ∗​，则 $\mathbf W^{-1}\mathbf B$W−1B 的特征向量为 $\mathbf W^{-\frac 12}v_\ell^*$W−21​vℓ∗​．

Fisher 通过不同的方式得到这个分解，完全没有引用高斯分布．他提出下面的问题：

寻找线性组合 $Z=a^TX$Z=aTX 使得组间方差相对于组内方差最大化．

再一次，组间方差是 Z 的均值的方差，而组内方差是关于均值的联合方差．图 4.9 显示了为什么这个准则是有意义的．尽管连接形心的方向能将均值尽可能分离开（比如，使得组间方差最大化），但是由于协方差的本性在投影类别上有很大重叠．同时考虑协方差，可以找到最小化重叠的方向．

图4.9. 尽管连接形心的直线定义了最大形心分散的方向，但是由于协方差（左图）投影数据会发生重叠．判别边界的方向使得高斯数据的重叠最小（右图）．

$Z$Z 的组间方差为 $a^T\mathbf Ba$aTBa，而组内方差为 $a^T\mathbf Wa$aTWa，$\mathbf W$W 是很早定义的，$\mathbf B$B 是类别形心矩阵 $\mathbf M$M 的协方差矩阵．注意到 $\mathbf {B+W=T}$B+W=T，其中 $\mathbf T$T 是 $\mathbf X$X 的总协方差矩阵，忽略掉了类别信息．

---

Hytn注：

这里 $\mathbf B$B（类间散度矩阵），$\mathbf W$W（类内散度矩阵） 和 $\mathbf T$T（全局散度矩阵） 均忽略掉了自由度．若有 $N$N 个观测，$K$K 个类别，且 $\bar x_i,\,\bar x,\,\bar x_{ij}$xˉi​,xˉ,xˉij​ 均为 $p$p 维向量，则有
$\begin{aligned} \mathbf B &= \sum\limits_{i=1}^Kn_i(\bar x_i-\bar x)(\bar x_i-\bar x)'\\ \mathbf W &= \sum\limits_{i=1}^K\sum\limits_{j=1}^{n_i}(x_{ij}-\bar x_i)(x_{ij}-\bar x_i)'\\ \mathbf T &= \sum\limits_{i=1}^K\sum\limits_{j=1}^{n_i}(x_{ij}-\bar x)(x_{ij}-\bar x)' \end{aligned}$BWT​=i=1∑K​ni​(xˉi​−xˉ)(xˉi​−xˉ)′=i=1∑K​j=1∑ni​​(xij​−xˉi​)(xij​−xˉi​)′=i=1∑K​j=1∑ni​​(xij​−xˉ)(xij​−xˉ)′​
它们自由度分别为 $K-1$K−1，$N-K$N−K 和 $N-1$N−1，易证 $\mathbf{T=B+W}$T=B+W．

---

Fisher 问题因此等价于最大化 Rayleigh quotient,
$\underset{a}{\max}\;\dfrac{a^T\mathbf Ba}{a^T\mathbf Wa}\tag{4.15}$amax​aTWaaTBa​(4.15)
或者等价地，
$\underset{a}{\max}\;a^T\mathbf Ba \; st \; a^T\mathbf Wa=1\tag{4.16}$amax​aTBastaTWa=1(4.16)
这是一个一般化的特征值问题，$a$a 是由 $\mathbf W^{-1}\mathbf B$W−1B 的最大特征值给出．不难证明（练习 4.1）最优 $a_1$a1​ 等于上面定义的 $v_1$v1​．类似地，可以找到下一个方向 $a_2$a2​，在 $W$W 中与 $a_1$a1​ 正交，使得 $a_2^TBa_2/a_2^TWa_2$a2T​Ba2​/a2T​Wa2​ 最大化；解为 $a_2=v_2$a2​=v2​，其余类似．$a_\ell$aℓ​ 被称作 判别坐标 (discriminant coordinates)，不会与判别函数相混淆．他们也被称作 典则变量 (canonical variables)，因为这些结果的一个变形是通过在预测变量矩阵 $X$X 上对指示响应矩阵 $Y$Y 进行典则相关分析得到的．这一点将在 12.5 节继续讨论．

!!! info “weiya 注：Ex. 4.1”
已解决，详见 Issue 142: Ex. 4.1．

个人解读

对多类LDA补充一下，a是低维空间基向量，可以以矩阵的形式也可以以向量的形式表示，有关更具体的LDA原理可参考这篇文章。

上面是多类LDA的原理，下面给出二类LDA原理：

原始第i类数据均值$\mu_i$μi​表达式如下
$\mu_{i}=\frac{1}{N_{i}} \sum_{x \in \omega_{i}} x$μi​=Ni​1​x∈ωi​∑​x
投影后的第i类数据均值表达式如下
$\begin{aligned} \widetilde{\mu}_{i}=\frac{1}{N_{i}} \sum_{y \in \omega_{i}} y &=\frac{1}{N_{i}} \sum_{x \in \omega_{i}} w^{T} x \\ &=w^{T} \frac{1}{N_{i}} \sum_{x \in \omega_{i}} x=w^{T} \mu_{i} \end{aligned}$μ​i​=Ni​1​y∈ωi​∑​y​=Ni​1​x∈ωi​∑​wTx=wTNi​1​x∈ωi​∑​x=wTμi​​
投影后的第i类数据方差表达式如下
$\begin{aligned} \widetilde{s}_{i}^{2}=\sum_{y \in \omega_{i}}\left(y-\widetilde{\mu}_{i}\right)^{2} &=\sum_{x \in \omega_{i}}\left(w^{T} x-w^{T} \mu_{i}\right)^{2} \\ &=\sum_{x \in \omega_{i}} w^{T}\left(x-\mu_{i}\right)\left(x-\mu_{i}\right)^{T} w \\ &=w^{T} S_{i} w \end{aligned}$si2​=y∈ωi​∑​(y−μ​i​)2​=x∈ωi​∑​(wTx−wTμi​)2=x∈ωi​∑​wT(x−μi​)(x−μi​)Tw=wTSi​w​
因为目标函数为
$J(w)=\frac{\left|\widetilde{\mu}_{1}-\widetilde{\mu}_{2}\right|^{2}}{\widetilde{S}_{1}^{2}+\widetilde{S}_{2}^{2}}$J(w)=S12​+S22​∣μ​1​−μ​2​∣2​
上式中的分子表示不同类别均值之差，分母表示不同类别方差之和，目标就是最大化目标函数，这和书中的表述也是一致的。

对于分母
$\widetilde{s_{1}}^{2}+\widetilde{s_{2}^{2}}=w^{T} S_{1} w+w^{T} S_{2} w=w^{T}\left(S_{1}+S_{2}\right) w=w^{T} W w=\widetilde{S}_{W}$s1​​2+s22​​=wTS1​w+wTS2​w=wT(S1​+S2​)w=wTWw=SW​
对于分子
$\begin{aligned} \left(\widetilde{\mu}_{1}-\widetilde{\mu}_{2}\right)^{2}=&\left(w^{T} \mu_{1}-w^{T} \mu_{2}\right)^{2} \\ &=w^{T} \left(\mu_{1}-\mu_{2}\right)\left(\mu_{1}-\mu_{2}\right)^{T} w \\ &=w^{T} B w=\widetilde{S}_{B} \end{aligned}$(μ​1​−μ​2​)2=​(wTμ1​−wTμ2​)2=wT(μ1​−μ2​)(μ1​−μ2​)Tw=wTBw=SB​​
所以目标函数，也就是广义瑞利商最终为
$J(w)=\frac{\left|\widetilde{\mu}_{1}-\widetilde{\mu}_{2}\right|^{2}}{\widetilde{s}_{1}^{2}+\widetilde{s}_{2}^{2}}=\frac{w^{T} B w}{w^{T} W w}$J(w)=s12​+s22​∣μ​1​−μ​2​∣2​=wTWwwTBw​
上面的表示是为了符合原书中的表示，这里的$B=S_B$B=SB​，$W=S_W$W=SW​那如何求解目标函数最大值呢？
$\begin{aligned} \frac{d}{d w} J(w) &=\frac{d}{d w}\left(\frac{w^{T} S_{B} w}{w^{T} S_{W} w}\right)=0 \\ & \Rightarrow\left(w^{T} S_{W} w\right) \frac{d}{d w}\left(w^{T} S_{B} w\right)-\left(w^{T} S_{B} w\right) \frac{d}{d w}\left(w^{T} S_{W} w\right)=0 \\ & \Rightarrow\left(w^{T} S_{W} w\right) 2 S_{B} w-\left(w^{T} S_{B} w\right) 2 S_{W} w=0 \end{aligned}$dwd​J(w)​=dwd​(wTSW​wwTSB​w​)=0⇒(wTSW​w)dwd​(wTSB​w)−(wTSB​w)dwd​(wTSW​w)=0⇒(wTSW​w)2SB​w−(wTSB​w)2SW​w=0​
对上式进行继续变换
$\begin{array}{l} {\Rightarrow\left(\frac{w^{T} S_{W} w}{w^{T} S_{W} w}\right) S_{B} w-\left(\frac{w^{T} S_{B} w}{w^{T} S_{W} w}\right) S_{W} w=0} \\ {\Rightarrow S_{B} w-J(w) S_{W} w=0} \\ {\Rightarrow S_{W}^{-1} S_{B} w-J(w) w=0} \end{array}$⇒(wTSW​wwTSW​w​)SB​w−(wTSW​wwTSB​w​)SW​w=0⇒SB​w−J(w)SW​w=0⇒SW−1​SB​w−J(w)w=0​
这样问题就转变为求解特征值问题了
$W^{-1} B w=\lambda w \quad \text { where } \quad \lambda=J(w)=\text {scalar}$W−1Bw=λw where λ=J(w)=scalar
其中
$w^{*}=\underset{w}{\arg \max } J(w)=\underset{w}{\arg \max }\left(\frac{w^{T} B w}{w^{T} W w}\right)=w^{-1}\left(\mu_{1}-\mu_{2}\right)$w∗=wargmax​J(w)=wargmax​(wTWwwTBw​)=w−1(μ1​−μ2​)
总结LDA算法总体思路如下

输入：数据集 $D = \{(x_1, y_1), \dots ,(x_m, y_m)\}$D={(x1​,y1​),…,(xm​,ym​)}，任意样本$x_i$xi​为n维向量，$\mathrm{y}_{\mathrm{i}} \in\{\mathrm{C}_1, \ldots, \mathrm{C_k}\}$yi​∈{C1​,…,Ck​}，共k个类别。现在要将其降维到d维；
输出：降维后的数据集D’。

（1）计算类内散度矩阵 B;

（2）计算类间散度矩阵 W;

（3）计算矩阵$W^{-1}B$W−1B;

（4）计算$W^{-1}B$W−1B的最大的d个特征值和对应的d个特征向量(w1,w2,…wd)，得到投影矩阵$w$w;

（5）对样本集中的每一个样本特征$x_i$xi​，转化为新的样本$z_i=w^Tx_i$zi​=wTxi​;

（6）得到输出样本集D′={(z1,y1),(z2,y2),…,((zm,ym))}。

继续原文

总结一下至今为止的发展：

• 有相同的协方差矩阵的高斯分类导出线性判别边界．分类可以通过对数据关于 $\mathbf W$W 球面化得到，并且划分到球空间的最近形心内（矫正因子 $\log\pi_k$logπk​）
• 因为只计算了到形心的相对距离，所以可以把数据局限于在球空间的形心张成的子空间．
• 子空间可以进一步分解为关于形心分离的最优子空间．这个分解与 Fisher 的分解相同．

降维后的子空间可以看成是数据降维（为了观察）的一个工具．它们是否可以用于分类以及它的基本原理是什么？很明显它们可以用在分类中，正如我们在最初的变形中一样；我们可以简单地把到形心的距离计算限制到选定的子空间中．可以证明这是有着额外限制条件——高斯形心位于 $\rm{IR}^p$IRp 的 $L$L 维子空间中——的高斯分类器．通过极大似然法拟合这样一个模型，然后运用贝叶斯定理构造后验概率，恰巧是上面描述的分类准则．（练习 4.8）

高斯分类器在距离计算时要求矫正因子 $\log\pi_k$logπk​．使用矫正因子的理由可以在图 4.9 中看出．错分类率是基于两个密度计算的重叠部分的面积．如果 $\pi_k$πk​ 是相等的（图中隐含了），则最优的分离点是在投射均值之间．如果 $\pi_k$πk​ 不相等，朝着最小类别移动分类点会改善错误率．正如之前提到的两个类别，可以通过使用 LDA（或者其他任何方法）导出线性规则，然后选择分类点去最小化训练集上的误判率．

作为一个展现降维限制优点的例子，我们回到元音数据．这里有 11 个类别 10 个变量，因此该分类器有 10 个可能的维数．我们可以在每个层次空间计算训练和测试误差；图 4.10 显示了这个结果．图 4.11 显示了基于二维 LDA 的解的分类器的判别边界．

图 4.10. 对于元音数据，训练和测试误差作为判别边界的维数的函数的图象．这种情况下最优的误差率是维数等于 2 的情况．图 4.11 显示了这个空间的判别边界．

Fisher 降秩判别分析 (RDA) 和指示响应矩阵回归之间存在着紧密的联系．事实表明 LDA 意味着 $\hat{\mathbf Y}^T\mathbf Y$Y^TY 的特征值分解后进行回归．在两个类的情形下，存在一个单判别变量，乘上一个标量后能与 $\hat{\mathbf Y}$Y^ 的某一列相等．这些联系将在第 12 章中讨论．一个相关的事实是先将原始的预测变量 $\mathbf X$X 转换为 $\hat{\mathbf Y}$Y^，然后用 $\hat{\mathbf Y}$Y^ 做 LDA 与在原空间中做 LDA 是相同的（练习 4.3）．

图 4.11. 对于元音训练数据，由前两个典则变量张成的二维子空间中的判别边界．注意到在任何高维子空间下，判别边界是高维仿射平面，而且不可以表示成直线．

---

1. Hastie, T., Tibshirani, R. and Buja, A. (1994). Flexible discriminant analysis by optimal scoring, Journal of the American Statistical Association 89: 1255–1270. ↩︎

2. Michie, D., Spiegelhalter, D. and Taylor, C. (eds) (1994). Machine Learning, Neural and Statistical Classification, Ellis Horwood Series in Artificial Intelligence, Ellis Horwood. ↩︎

3. Friedman, J. (1989). Regularized discriminant analysis, Journal of the American Statistical Association 84: 165–175. ↩︎