线性判别分析（LDA）的原理

前言

LDA在模式识别领域（比如人脸识别，舰艇识别等图形图像识别领域）中有非常广泛的应用，因此我们有必要了解一下它的算法原理。除非特别声明，本文中的LDA均指的是线性判别分析（Linear Discriminant Analysis），它与自然语言处理领域中的LDA隐含狄利克雷分布（Latent Dirichlet Allocation）是有本质上的区别的，后者是一种处理文档的主题模型。

一、LDA的思想

LDA是一种监督学习的降维技术，也就是说它的数据集的每个样本是有类别输出的。这点和PCA不同。PCA是不考虑样本类别输出的无监督降维技术。LDA的思想可以用一句话概括，就是“投影后类内方差最小，类间方差最大”。什么意思呢？ 我们要将数据在低维度空间上进行投影，投影后希望同种类别数据的投影点尽可能的接近，而不同类别的数据的类别中心之间的距离尽可能的大。

这里假设我们有两类数据分别用红色和蓝色表示，如下图所示，这些数据特征是二维的，我们希望将这些数据投影到一维的一条直线上，让同种类别数据的投影点尽可能的接近，而且让红色和蓝色数据中心之间的距离尽可能的大：

上图中提供了两种投影方式，那么哪一种能更好的满足我们的标准呢？从直观上可以看出，右图要比左图的投影效果好，因为右图的红色数据和蓝色数据各自较为集中，且类别之间的距离明显。左图则在边界处两类数据比较混杂。以上就是LDA的主要思想了，当然在实际应用中，我们的数据一般是多个类别的，我们的原始数据一般也是超过二维的，投影后的数据也一般不是直线，而是一个低维的超平面。

在我们将上面直观的内容转化为可以度量的问题之前，我们先了解些必要的数学基础知识，这些在后面讲解具体LDA原理时会用到。

二、瑞利商（Rayleigh Quotient）和广义瑞利商（Generalized Rayleigh Quotient）

瑞利商是指这样的函数 $R(A,x)$ ：
$$R(A,x)=\frac{x^{H}Ax}{x^{H}x}$$
其中 $x$ 为非零向量，而 $A$ 为 $n\times n$ 的 $Hermite$ （埃米尔特）矩阵。所谓的 $Hermite$ 矩阵就是满足共轭转置矩阵和自己相等的矩阵，即 $A^H=A$ 。如果我们的矩阵 $A$ 是实矩阵，则满足 $A^H=A^T$ 的矩阵即为 $Hermite$ 矩阵。

瑞利商 $R(A,x)$ 有一个非常重要的性质，即它的最大值等于矩阵 $A$ 最大的特征值，而最小值等于矩阵 $A$ 的最小的特征值，也就是满足：
$${\lambda }_{min}\le \frac{x^{H}Ax}{x^{H}x}\le {\lambda }_{max}$$
具体的证明大家可以参考这篇文章：瑞利商及其极值的计算。
当向量 $x$ 是标准正交基时，即满足 $x^{H}x=1$ 时，瑞利商退化为：$R(A,x)=x^{H}Ax$ ，这个形式在谱聚类和PCA中都有出现。

以上就是瑞利商的内容，现在我们再看看广义瑞利商。广义瑞利商是指这样的函数 $R(A,B,x)$ ：
$$R(A,B,x)=\frac{x^{H}Ax}{x^{H}Bx}$$
其中 $x$ 为非零向量，而 $A,B$ 为 $n\times n$ 的 $Hermite$ 矩阵。$B$ 为正定矩阵。它的最大值和最小值是什么呢？其实我们只要通过将其标准化就可以转化为瑞利商的格式。我们令 $x=B^{-1/2}{x}^{\prime }$ ,则分母转化为：
$$x^{H}Bx=x^{\prime H}(B^{-1/2}{)}^{H}B{B}^{-1/2}{x}^{\prime }=x^{\prime H}{B}^{-1/2}B{B}^{-1/2}{x}^{\prime }=x^{\prime H}{x}^{\prime }$$
而分子转化为：
$$x^{H}Ax=x^{\prime H}{B}^{-1/2}A{B}^{-1/2}{x}^{\prime }$$
此时我们的 $R(A,B,x)$ 转化为 $R(A,B,x\prime )$：
$$R(A,B,x^{\prime })=\frac{x^{\prime H}{B}^{-1/2}A{B}^{-1/2}{x}^{\prime }}{x^{\prime H}{x}^{\prime }}$$
利用前面的瑞利商的性质，我们可以很快的知道，$R(A,B,x^{\prime })$ 的最大值为矩阵 $B^{-1/2}A{B}^{-1/2}$ 的最大特征值，或者说矩阵 $B^{-1}A$ 的最大特征值，而最小值为矩阵 $B^{-1}A$ 的最小特征值。这里使用了一些技巧，即对矩阵进行标准化 $\frac{A_{ij}}{\sqrt{B_i\cdot B_j}}$。

三、二类LDA的原理

现在我们回到LDA的原理上，我们在第一小节介绍到了LDA希望投影后同种类别数据的投影点尽可能的接近，而不同类别的数据的类别中心之间的距离尽可能的大，但这只是一个感官的度量。现在我们首先从比较简单的二类LDA入手，严谨的分析LDA的原理。

假设我们的数据集 D = { ( x 1 , y 1 ) , ( x 2 , y 2 ) , . . . , ( ( x m , y m ) ) } D=\{(x_1, y_1),(x_2, y_2), ... , ((x_m,y_m))\} D={(x1​,y1​),(x2​,y2​),...,((xm​,ym​))} ，其中任意样本 $x_i$ 为 $n$ 维向量， y i ∈ { 0 , 1 } y_i∈\{0,1\} yi​∈{0,1} 。我们定义 $N_j(j=0,1)$ 为第 $j$ 类样本的个数，$X_j(j=0,1)$ 为第 $j$ 类样本的集合，而 ${\mu }_j(j=0,1)$ 为第 $j$ 类样本的均值向量，定义 ${\Sigma }_j(j=0,1)$ 为第 $j$ 类样本的协方差矩阵（严格说是缺少分母部分的协方差矩阵）。

${\mu }_j$ 的表达式为：
$${\mu }_j=\frac{1}{N_j}\sum _{x\in X_j}x(j=0,1)$$
${\Sigma }_j$ 的表达式为：
$${\Sigma }_j=\sum _{x\in X_j}(x-{\mu }_j)(x-{\mu }_j{)}^T(j=0,1)$$
由于是两类数据，因此我们只需要将数据投影到一条直线上即可。假设我们的投影直线是向量 $w$ ，则对任意一个样本 $x_i$，它在直线 $w$ 的投影为 $w^T{x}_i$，对于我们的两个类别的中心点 ${\mu }_0,{\mu }_1$ ,在直线 $w$ 上的投影为 $w^T{\mu }_0$ 和 $w^T{\mu }_1$ 。由于LDA需要让不同类别的数据的类别中心之间的距离尽可能的大，也就是我们要最大化 $||w^T{\mu }_0-w^T{\mu }_1|{|}_2^2$ ,同时我们希望同种类别数据的投影点尽可能的接近，也就是要让同类样本投影点的协方差 $w^T{\Sigma }_{0}w$ 和 $w^T{\Sigma }_{1}w$ 尽可能的小，即最小化 $w^T{\Sigma }_{0}w+w^T{\Sigma }_{1}w$ 。综上所述，我们的优化目标为：
$$\underset{w}{\underbrace{\arg \max }}J(w)=\frac{||w^T{\mu }_0-w^T{\mu }_1|{|}_2^2}{w^T{\Sigma }_{0}w+w^T{\Sigma }_{1}w}=\frac{w^T({\mu }_0-{\mu }_1)({\mu }_0-{\mu }_1{)}^{T}w}{w^T({\Sigma }_0+{\Sigma }_1)w}$$
我们一般定义类内散度矩阵 $S_w$ 为：
$$S_w={\Sigma }_0+{\Sigma }_1=\sum _{x\in X_0}(x-{\mu }_0)(x-{\mu }_0{)}^T+\sum _{x\in X_1}(x-{\mu }_1)(x-{\mu }_1{)}^T$$
同时定义类间散度矩阵 $S_b$为：
$$S_b=({\mu }_0-{\mu }_1)({\mu }_0-{\mu }_1{)}^T$$
这样我们的优化目标重写为：
$$\underset{w}{\underbrace{\arg \max }}J(w)=\frac{w^T{S}_{b}w}{w^T{S}_{w}w}$$
仔细观察上式，这不就是我们的广义瑞利商嘛！那么这就简单了，利用我们第二小节讲到的广义瑞利商的性质，我们知道我们的 $J(w^{\prime })$ 最大值为矩阵 $S_w^{-1/2}{S}_b{S}_w^{-1/2}$ 的最大特征值，而对应的 $w^{\prime }$ 其实为 $S_w^{-1/2}{S}_b{S}_w^{-1/2}$ 的最大特征值对应的特征向量！而 $S_w^{-1}{S}_b$ 的特征值和 $S_w^{-1/2}{S}_b{S}_w^{-1/2}$ 的特征值相同， $S_w^{-1}{S}_b$ 的特征向量 $w$ 和 $S_w^{-1/2}{S}_b{S}_w^{-1/2}$ 的特征向量 $w^{\prime }$ 满足 $w=S_w^{-1/2}{w}^{\prime }$ 的关系！

注意到对于二类的时候，$S_{b}w$ 的方向恒为 ${\mu }_0-{\mu }_1$ ，因为：
$$S_{b}w=({\mu }_0-{\mu }_1)({\mu }_0-{\mu }_1{)}^{T}w$$
$$\iff S_{b}w=({\mu }_0-{\mu }_1)(w^T({\mu }_0-{\mu }_1){)}^T=({\mu }_0-{\mu }_1)(w^T{\mu }_0-w^T{\mu }_1{)}^T$$
我们在之前注意到了 $w^T{\mu }_0-w^T{\mu }_1$ 其实是一个标量，即上式可以表示为：
$$S_{b}w=C({\mu }_0-{\mu }_1)\text{其中C是一个常数}$$
不妨令 $S_{b}w=\lambda ({\mu }_0-{\mu }_1)$ ，将其带入：$(S_w^{-1}{S}_b)w=\lambda w$ ，可以得到 $w=S_w^{-1}({\mu }_0-{\mu }_1)$ ， 也就是说我们只要求出原始二类样本的均值和方差就可以确定最佳的投影方向 $w$ 了。

四、多类LDA原理

有了二类LDA的基础，我们再来看看多类别LDA的原理。

假设我们的数据集 D = { ( x 1 , y 1 ) , ( x 2 , y 2 ) , . . . , ( ( x m , y m ) ) } D=\{(x_1, y_1), (x_2, y_2), ..., ((x_m, y_m))\} D={(x1​,y1​),(x2​,y2​),...,((xm​,ym​))}，其中任意样本 $x_i$ 为 $n$ 维向量， y i ∈ { C 1 , C 2 , . . . , C k } y_i \in \{C_1, C_2, ..., C_k\} yi​∈{C1​,C2​,...,Ck​} 。我们定义 $N_j(j=1,2,...,k)$ 为第j类样本的个数，$X_j(j=1,2,...,k)$ 为第 $j$ 类样本的集合，而 ${\mu }_j(j=1,2,...,k)$ 为第 $j$ 类样本的均值向量，定义 ${\Sigma }_j(j=1,2,...,k)$ 为第 $j$ 类样本的协方差矩阵。在二类LDA里面定义的公式我们可以类推到多类LDA。

由于我们是多类向低维空间投影，则此时投影到的低维空间就不是一条直线了，而是一个超平面。假设我们投影到的低维空间的维度为 $d$ ，对应的基向量为 $(w_1,w_2,...,w_d)$ ，基向量组成的矩阵为 $W$ , 它是一个 $n\times d$ 的矩阵。

此时我们的优化目标可以变成为：
$$J(W)=\frac{W^T{S}_{b}W}{W^T{S}_{w}W}$$
其中 $S_b={\sum }_{j=1}^k{N}_j({\mu }_j-\mu )({\mu }_j-\mu {)}^T$ ，$\mu$为所有样本的均值向量， $S_w={\sum }_{j=1}^k{S}_{wj}={\sum }_{j=1}^k{\sum }_{x\in X_j}(x-{\mu }_j)(x-{\mu }_j{)}^T$ 。
但是这里有一个问题，就是 $W^T{S}_{b}W$ 和 $W^T{S}_{w}W$ 都是矩阵，不是标量，我们无法把它作为一个标量函数来优化！也就是说，我们无法直接用二类LDA的优化方法，怎么办呢？一般来说，我们可以用其他的一些优化目标来进行替换。

常见的一个LDA多类优化目标函数定义为：
$$\underset{W}{\underbrace{\arg \max }}J(W)=\frac{{\prod }_{diag}{W}^T{S}_{b}W}{{\prod }_{diag}{W}^T{S}_{w}W}$$
其中 ${\prod }_{diag}A$ 为 $A$ 的主对角线元素的乘积，$W$ 为 $n\times d$ 的矩阵。
$J(W)$ 的优化过程可以转化为：
$$J(W)=\frac{{\prod }_{i=1}^d{W}_i^T{S}_b{W}_i}{{\prod }_{i=1}^d{W}_i^T{S}_w{W}_i}=\prod _{i=1}^d\frac{W_i^T{S}_b{W}_i}{W_i^T{S}_w{W}_i}$$
仔细观察上式最右边的等式，这不就是广义瑞利商嘛！最大值是矩阵 $S_w^{-1}{S}_b$ 的最大特征值，最大的 $d$ 个这样的广义瑞利商的值的乘积就是矩阵 $S_w^{-1}{S}_b$ 的最大的 $d$ 个特征值的乘积，此时对应的矩阵 $W$ 为这最大的 $d$ 个特征值对应的特征向量张成的矩阵。

由于 $W$ 是一个利用了样本的类别得到的投影矩阵，因此它降维能降到的维度 $d$ 的最大值为 $k-1$ 。为什么最大维度不是类别数 $k$ 呢？因为 $S_b$ 中每个 $({\mu }_j-\mu )({\mu }_j-\mu {)}^T$ 的秩为 $1$ ，因此协方差矩阵相加后最大的秩为 $k$ ( $Rank(A+B)\le Rank(A)+Rank(B)$ )，但是由于如果我们知道前 $k-1$ 个 ${\mu }_j$ 后，最后一个 ${\mu }_k$ 可以由前 $k-1$ 个 ${\mu }_j$ 线性表示（因为这里我们给定样本之后，所有样本的均值 $\mu$ 是已知的），因此 $S_b$ 的秩最大为 $k-1$ ，即特征向量最多有 $k-1$ 个。

五、LDA算法实现的流程

在第三小节和第四小节我们给大家介绍了LDA的原理，现在我们对LDA降维的流程做一个总结。

输入：数据集 D = { ( x 1 , y 1 ) , ( x 2 , y 2 ) , . . . , ( ( x m , y m ) ) } D=\{(x_1, y_1), (x_2, y_2), ..., ((x_m, y_m))\} D={(x1​,y1​),(x2​,y2​),...,((xm​,ym​))}，其中任意样本 $x_i$ 为 $n$ 维向量， y i ∈ { C 1 , C 2 , . . . , C k } y_i \in \{C_1, C_2, ..., C_k\} yi​∈{C1​,C2​,...,Ck​} ，降维到的维度 $d$ 。

输出：降维后的样本集$D^{\prime }$ 。

• ① 计算类内散度矩阵 $S_w$；

• ② 计算类间散度矩阵 $S_b$ ；

• ③ 计算矩阵 $S_w^{-1}{S}_b$ ；

• ④ 计算 $S_w^{-1}{S}_b$ 的最大的 $d$ 个特征值和对应的 $d$ 个特征向量 $(w_1,w_2,...,w_d)$ ，得到投影矩阵 $W$ ；

• ⑤ 对样本集中的每一个样本特征 $x_i$ ，转化成为新的样本 $z_i=W^T{x}_i$ ；

• ⑥ 得到输出样本集 D ′ = { ( z 1 , y 1 ) , ( z 2 , y 2 ) , . . . , ( ( z m , y m ) ) } D′=\{(z_1, y_1), (z_2, y_2), ..., ((z_m, y_m))\} D′={(z1​,y1​),(z2​,y2​),...,((zm​,ym​))} 。

以上就是使用LDA进行降维的算法流程。实际上LDA除了可以用于降维以外，还可以用于分类。一个常见的LDA分类基本思想是假设各个类别的样本数据符合高斯分布，这样利用LDA进行投影后，可以利用极大似然估计计算各个类别投影数据的均值和方差，进而得到该类别高斯分布的概率密度函数。当一个新的样本到来后，我们可以将它投影，然后将投影后的样本特征分别带入各个类别的高斯分布概率密度函数，计算它属于这个类别的概率，最大的概率对应的类别即为该样本的预测类别。

六、LDA与PCA之间的异同

LDA用于降维，和PCA有很多相同，也有很多不同的地方，因此值得我们好好的比较一下两者的异同点。

首先我们看看相同点：

• ① 两者均可以对数据进行降维。

• ② 两者在降维时均使用了矩阵特征分解的思想。

• ③ 两者都假设数据符合高斯分布。

不同点：

• ① LDA是有监督的降维方法，而PCA是无监督的降维方法

• ② LDA降维最多降到类别数k-1的维数，而PCA没有这个限制。

• ③ LDA除了可以用于降维，还可以用于分类。

-④ LDA选择分类性能最好的投影方向，而PCA选择样本点具有最大方差的投影方向。

在某些数据分布下LDA比PCA的降维效果更优，如下图所示：

当然，某些某些数据分布下PCA比LDA的降维效果更优，如下图所示：

七、总结

LDA算法既可以用来降维，又可以用来分类，但是目前来说，主要还是用于降维。在我们进行图像识别图像识别相关的数据分析时，LDA是一个有力的工具。下面总结下LDA算法的优缺点。

LDA算法的主要优点有：

• ① 在降维过程中可以使用类别的先验知识经验，而像PCA这样的无监督学习则无法使用类别先验知识。

• ② LDA在样本分类信息依赖均值而不是方差的时候，优于PCA之类的降维算法。

LDA算法的主要缺点有：

• ① LDA不适合对非高斯分布样本进行降维，PCA也有这个问题。

• ② LDA降维最多降到类别数k-1的维数，如果我们降维的维度大于k-1，则不能使用LDA。当然目前有一些LDA的进化版算法可以绕过这个问题。

• ③ LDA在样本分类信息依赖方差而不是均值的时候，降维效果不是很好。

• ④ LDA可能过度拟合数据。

八、参考文献

① 线性判别分析LDA原理总结 - 刘建平
② 线性代数之矩阵的属性和运算