机器学习算法03-线性判别分析LDA

​ 判别分析（Discriminant Analysis）是一种分类方法，它通过一个已知类别的“训练样本”来建立判别准则，并通过预测变量来为未知类别的数据进行分类。线性判别式分析（Linear Discriminant Analysis，简称为LDA）是其中一种，也是模式识别的经典算法。LDA以Bayes判别思想为基础，当分类只有两种且总体服从多元正态分布条件下，Bayes判别与Fisher判别、距离判别是等价的。

1. LDA和PCA比较

LDA与PCA主成分分析都是常用的降维方法，二者的区别在于：

1. 出发思想不同。PCA主要是从特征的协方差角度，去找到比较好的投影方式，即选择样本点投影具有最大方差的方向；而LDA则更多的是考虑了分类标签信息，寻求投影后不同类别之间数据点距离更大化以及同一类别数据点距离最小化，即选择分类性能最好的方向。
2. 学习模式不同。PCA属于无监督式学习，因此大多场景下只作为数据处理过程的一部分，需要与其他算法结合使用，例如将PCA与聚类、判别分析、回归分析等组合使用；LDA是一种监督式学习方法，本身除了可以降维外，还可以进行预测应用，因此既可以组合其他模型一起使用，也可以独立使用。
3. 降维后可用维度数量不同。LDA降维后最多可生成N-1维子空间（分类标签数-1），因此LDA与原始维度数量无关，只有数据标签分类数量有关；而PCA最多有n维度可用，即最大可以选择全部可用维度。

---

从直接可视化的角度，以二维数据降维为例，PCA和LDA的区别如下图:

上图左侧是PCA的降维思想，它所作的只是将整组数据整体映射到最方便表示这组数据的坐标轴上，映射时没有利用任何数据内部的分类信息。因此，虽然PCA后的数据在表示上更加方便（降低了维数并能最大限度的保持原有信息），但在分类上也许会变得更加困难；上图右侧是LDA的降维思想，可以看到LDA充分利用了数据的分类信息，将两组数据映射到了另外一个坐标轴上，使得数据更易区分了（在低维上就可以区分，减少了运算量）。

线性判别分析LDA算法由于其简单有效性在多个领域都得到了广泛地应用，是目前机器学习、数据挖掘领域经典且热门的一个算法；但是算法本身仍然存在一些局限性：

• 当样本数量远小于样本的特征维数，样本与样本之间的距离变大使得距离度量失效，使LDA算法中的类内、类间离散度矩阵奇异，不能得到最优的投影方向，在人脸识别领域中表现得尤为突出
• LDA不适合对非高斯分布的样本进行降维
• LDA在样本分类信息依赖方差而不是均值时，效果不好
• LDA可能过度拟合数据

2. LDA思想

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

​ 可能还是有点抽象，我们先看看最简单的情况。假设我们有两类数据 分别为红色和蓝色，如下图所示，这些数据特征是二维的，我们希望将这些数据投影到一维的一条直线，让每一种类别数据的投影点尽可能的接近，而红色和蓝色数据中心之间的距离尽可能的大。

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

3. 瑞利熵

瑞利商是指这样的函数R(A,x):

R(A,x)=xTAxxTx

​ 其中x 为非零向量，而A为 n×n的Hermitian矩阵。所谓的Hermitian矩阵就是满足共轭转置矩阵和自己相等的矩阵，即AH=A。如果我们的矩阵A是实矩阵，则满足AT=A 的矩阵即为Hermitian矩阵。

​ 瑞利商R(A,x) 有一个非常重要的性质，即它的最大值等于矩阵A最大的特征值，而最小值等于矩阵AA的最小的特征值，也就是满足

λmin⪯xHAxxHx⪯λmax

​ 当向量 x 是标准正交基时，即满足 xHx=1 时，瑞利商退化为: R(A,x)=xHAx .

​ 广义瑞利商。广义瑞利商是指这样的函数 R(A,B,x) :

R(A,B,x)=xHAxxHBx

​ 其中 x 为非零向量，而A, B 为n×n的Hermitian矩阵。B为正定矩阵。它的最大值和最小值是什么呢？其实我们只要通过将其通过标准化就可以转化为瑞利商的格式。我们令x'=B−1/2x, 则分母转化为:

xHBx=x′H(B−1/2)HBB−1/2x′=x′HB−1/2BB−1/2x′=x′Hx′

​ 而分子转化为:

xHAx=x′HB−1/2AB−1/2x′

​ 此时我们的R(A,B,x)转化为R(A,B,x'):

R(A,B,x′)=x′HB−1/2AB−1/2x′x′Hx′

​ 利用前面的瑞利商的性质，我们可以很快的知道，R(A,B,x)的最大值为矩阵B−1/2AB−1/2的最大特征值，或者说矩阵B−1A的最大特征值，而最小值为矩阵B−1A的最小特征值。

4. 二类LDA

​ 训练数据：D=(x1,y1),(x2,y2),...,(xn,yn) , D 中有n个数据, 其中xi 是一个d维的向量, xi=(xi1,xi2,...,xid)T , y 在 {0, 1} 中取值, 其中 1表示正例, 0表示负例.

​ 我们定义Nj(j=0,1) 为第j类样本的个数，xj(j=0,1) 为第j类样本的集合，而μj(j=0,1)为第j类样本的均值向量，定义Σj(j=0,1)为第j类样本的协方差矩阵.

​ μj的表达式为:

μj=1Nj∑x∈Xjx(j=0,1)

​ Σj的表达式为:

Σj=∑x∈Xj(x−μj)(x−μj)T(j=0,1)

​ 由于是两类数据，因此我们只需要将数据投影到一条直线上即可。假设我们的投影直线是向量w, 则对任意一个样本本xi, 它在直线w的投影为wTxi, 我们将这个最佳的向量称为 w (d 维)，那么样例x (d 维)到w上的投影可以用下式来计算: y=wTx, 这里得到的y值不是0/1值，而是x投影到直线上的点到原点的距离。

​ 对于我们的两个类别的中心点μ0, μ1, 在直线w的投影为wTμ0和wTμ1。由于LDA需要让不同类别的数据的类别中心之间的距离尽可能的大，也就是我们要最大化||wTμ0−wTμ1||22,同时我们希望同一种类别数据的投影点尽可能的接近，也就是要同类样本投影点的协方差wTΣ0w和wTΣ1w尽可能的小，即最小化wTΣ0w+wTΣ1w。综上所述，我们的优化目标为:

argmaxwJ(w)=||wTμ0−wTμ1||22wTΣ0w+wTΣ1w=wT(μ0−μ1)(μ0−μ1)TwwT(Σ0+Σ1)w

​ 定义”类内散度矩阵” (within-class scatter matrix)

Sw=Σ0+Σ=∑x∈X0(x−μ0)(x−u0)T+∑x∈X1(x−μ1)(x−μ1)T

​ 定义”类间散度矩阵” (between-class scatter matrix)

Sb=(μ0−μ1)(μ0−μ1)T=∑j=01Nj(μj−μ)(μj−μ)T

​ 这样我们的优化目标重写为:

argmaxwJ(w)=wTSbwwTSww

​ J(w) 的解和 的长度无关, 只与其方向有关, 令 :

argminwJ=−wTSbw

s.t.wTSww=1

​ 由拉格朗日乘子法, 可得:

⟹⟹⟹⟹L(w)=−wTSbw−λ(wTSww−1)∂L(w)∂w=−2∗Sbw−2λ∗Sww=0Sbw=λSwwSbw=λ(u0−u1)Sww=(u0−u1)

​ 利用广义瑞利商的性质，我们知道我们的J(w)最大值为矩阵S−1wSb的最大特征值，而对应的w为S_w^{−1}S_bS−1wSb的最大特征值对应的特征向量.

​ 注意到对于二类的时候，S_bwSbw的方向恒为μ0−μ1μ0−μ1, 因为S_bw=(\mu_0-\mu_1)(\mu_0-\mu_1)^Tw = \lambda(u_0-u_1)Sbw=(μ0−μ1)(μ0−μ1)Tw=λ(u0−u1), 可以得到:

w=S_w^{−1}(μ_0−μ_1)w=S−1w(μ0−μ1)

​ 也就是说我们只要求出原始二类样本的均值和方差就可以确定最佳的投影方向w了.

二维样本的投影的原图和结果图:

=>

5. 多类LDA

​ 数据集D=(x1,y1),(x2,y2),...,(xm,ym) , 样本xi为n 维向量，yi∈C1,C2,...,Ck。我们定义Nj (j=1,2...,k) 为第j类样本的个数，Xj(j=1,2...k)为第j类样本的集合，而μj(j=1,2...k)为第j类样本的均值向量，定义Σj(j=1,2...k)为第j 类样本的协方差矩阵。在二类LDA里面定义的公式可以很容易的类推到多类LDA。

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

​ 此时我们的优化目标应该可以变成为:

J(W)=WTSbWWTSwW

​ 其中μ 为所有样本均值向量.

Sb=∑j=1kNj(μj−μ)(μj−μ)T

Sw=∑j=1kSwj=∑j=1k∑x∈Xj(x−μj)(x−μj)T

​ 总体散度矩阵(total scatter matrix)

St=∑i=1m(xi−μ)(xi−μ)T=Sb+Sw

​ 显然多分类LDA可以有多种实现方法: 使用 Sb,Sw,St 中的任意两个即可实现优化.

​ 但是有一个问题，就是WTSbW和WTSwW 都是矩阵，不是标量，无法作为一个标量函数来优化！也就是说，我们无法直接用二类LDA的优化方法，怎么办呢？一般来说，我们可以用其他的一些替代优化目标来实现。

​ 常见的一个LDA多类优化目标函数定义为：

argmaxWJ(W)=∏diagWTSbW∏diagWTSwW

​ 其中∏diagA 为A的主对角线元素的乘积，W为m×d 的矩阵。

​ J(W) 的优化过程可以转化为：

J(W)=∏i=1dwTiSbwi∏i=1dwTiSwwi=∏i=1dwTiSbwiwTiSwwi

​ 上式最右边是广义瑞利商, 最大值是矩阵S−1wSb 的最大特征值, 最大的d个值的乘积就是矩阵S−1wSb的最大的d个特征值的乘积, 此时对应的矩阵W为这最大的d 个特征值对应的特征向量张成的矩阵。

​ 由于W 是一个利用了样本的类别得到的投影矩阵，因此它的降维到的维度d最大值为k−1 。为什么最大维度不是类别数k呢？因为Sb中每个μj−μ 的秩为1，因此协方差矩阵相加后最大的秩为kk(矩阵的秩小于等于各个相加矩阵的秩的和)，但是由于如果我们知道前k-1个μ_jμj 后，最后一个μ_kμk可以由前k-1个μ_jμj线性表示，因此S_bSb的秩最大为k-1，即特征向量最多有k-1个。

6. LDA算法流程

​ LDA降维的流程总结

​ 输入：D={(x_1,y_1),(x_2,y_2),...,(x_m,y_m)}D=(x1,y1),(x2,y2),...,(xm,ym) , x_ixi为n维向量，y_i∈{C_1,C_2,...,C_k}yi∈C1,C2,...,Ck，降维到的维度d。

​ 输出：降维后的样本集D'

​ 1) 计算类内散度矩阵Sw

​ 2) 计算类间散度矩阵Sb

​ 3) 计算矩阵S−1wSb

​ 4）计算S−1wSb的最大的d个特征值和对应的d个特征向量(w1,w2,...wd),得到投影矩阵W

​ 5) 对样本集中的每一个样本特征xi ,转化为新的样本zi=WTxi

​ 6) 得到输出样本集D'=(z1,y1),(z2,y2),...,(zm,ym)

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

---

参考文章:

1. http://www.cnblogs.com/pinard/p/6244265.html
2. http://www.cnblogs.com/jerrylead/archive/2011/04/21/2024384.html