机器学习入门学习笔记：（2.4）线性判别分析理论推导

LDA

线性判别分析(Linear Discriminant Analysis, 简称LDA)，最早由Fisher提出，也叫“Fisher判别分析”。

线性判别分析的思想：给定样本数据集，设法将样本投影到某一条直线上，使得同类样本的投影点尽可能接近，异类样本的投影点尽可能远；在对新的点进行分类预测时，将其投影到这条直线上，根据投影点的位置来判断样本的类别。

当x是二维时，我们就要寻找一个方向为ω的直线来使得这些样本点的投影分离。

二分类情况

先只考虑二分类情况，给出数据集 {xi,yi}mi=1,yi={0,1}。

设有N0个样本的标签yi=0，其类别为D0；N1个样本的标签yi=1，其类别为D1。

则分别对应类别D0和D1的均值向量μ0和μ1可以求出来：

μ0=1N0∑(x,y)∈D0xi μ1=1N1∑(x,y)∈D1xi

再给出一组新的度量值，称作散度值(scatter)： s0=∑(x,y)∈D0(xi−μ0)2s1=∑(x,y)∈D1(xi−μ1)2

这东西看上去很眼熟吧，不就是方差，少除了样本数量吗？因为这里只需要定量表示样本集合的分散程度，对常数系数不敏感，所以方差中上面的部分就足够了。

给出一组参数ω，假设值为：h(xi)=ωTxi，后面简写为hi=ωTxi。

在线性回归中，我们的目标是使得这个假设值等于样本的标签yi；

而线性判别分析中，假设值hi实质上是样本xi在w上的投影的长度；以二维情况考虑，就是样本xi对应的点(xi1,xi2)到直线ω的投影点到原点的长度。我们的目标是通过区分这个”长度“hi；

再放一次这幅图，看看图不难理解：

接着，根据ω可以求出投影后的样本均值μ0~和μ1~：

μ0~=1N0∑(x,y)∈D0hi=1N0∑(x,y)∈D0ωTxi=ωT1N0∑(x,y)∈D0xi=ωTμ0

同理得：μ1~=ωTμ1

还有，投影后的散度矩阵(scatter)：

s0~=1N0∑(x,y)∈D0(hi−μ0~)2=1N0∑(x,y)∈D0(ωTxi−ωTμ0)2=1N0∑(x,y)∈D0ωT(xi−μ0)2ω=ωT⎡⎣1N0∑(x,y)∈D0(xi−μ0)2⎤⎦ω=ωTs0ωT

同理得：s1~=ωTs1ω

回到最开始说的思想：希望同类样例的投影点尽可能接近，即散度矩阵尽可能小，也即s0~和s1~都要尽可能小，简单表示为s0~+s1~；异类样例的投影点尽可能远，也即均值差值尽可能大，(μ0~−μ1~)2尽可能大；

同时考虑上面两个条件，可以给出一个目标函数：

J=(μ0~−μ1~)2s0~+s1~

代入前面求出的μ0~、μ1~、s0~、s1~:

J=(ωTμ0−ωTμ1)2ωTs0ω+ωTs1ω=ωT(μ0−μ1)2ωωTs0ω+ωTs1ω=ωT(μ0−μ1)T(μ0−μ1)ωωT(s0+s1)ω

这坨东西看得挺复杂的，定义类内散度矩阵(within-class scatter matrix):

Sω=S0+S1=∑(x,y)∈D0(xi−μ0)2+∑(x,y)∈D1(xi−μ1)2=∑(x,y)∈D0(xi−μ0)T(xi−μ0)+∑(x,y)∈D1(xi−μ1)T(xi−μ1)

还有定义类间散度矩阵(between-class scatter matrix)：Sb=(μ0−μ1)2=(μ0−μ1)T(μ0−μ1)

则优化的目标函数变为：J=ωTSbωωTSωω

这个东西就是LDA希望最大化的目标，实质是Sb和Sω的“广义瑞利商”(generalized Rayleigh quotient)。

好的，接下来的任务就是通过J来确定最优的ω。

观察式子可以发现，J的大小不随ω的大小变化，上下都是ω的二次项相互抵消，结果而只与其方向有关。

所以可以令ωTSωω=1，则式子等价于：

minω−ωTSbωs.t.ωTSωω=1

使用拉格朗日乘子法：

设未知数λ，写出拉格朗日函数：

L(ω)=ωTSbω−λ(ωTSωω−1)

对拉格朗日函数求导，且导数为0：

dL(ω)dω=2Sbω−2λSωω=0

矩阵求导中ωTSbω可以简单看作是Sbω2。

得到结果：

Sbω=λSωω

这是一个典型的求矩阵特征值的问题。

从前面的公式：

J=ωTSbωωTSωω

我们知道ω的大小并不影响结果，而是方向才会影响结果。

另外，由这个式子：

Sbω=(μ0−μ1)(μ0−μ1)Tω

观察发现Sbω的方向恒为μ0−μ1，因为任取其中一个(μ0−μ1)T与ω相乘之后是常数，总会剩下一个μ0−μ1。

所以，设一个新的常量λω，使得：

Sbω=(μ0−μ1)λω

回到前面拉格朗日方程求导得到的结果：

Sbω=λSωω

得到：

(μ0−μ1)λω=λSωω

若ω可逆，则有：

ω=λωλS−1ω(μ0−μ1)

注意，由于最后结果只与ω的方向有关，与其大小无关；而这个结果前面的常量λωλ可以舍去：

ω=S−1ω(μ0−μ1)

这个就是最终的结果。

前面已经推导出了Sω、μ0、μ1，代入即可。我们只要有散度矩阵和均值即可求出最优的ω。

这里还有一点，考虑到数值解的稳定性，在实践中通常是对Sω进行奇异值分解，即Sω=UΣVT，这里Σ是一个实对角矩阵，其对角线上的元素是Sω的奇异值，然后再由S−1ω=VΣ−1UT得到Sω。（摘自西瓜书，没学矩阵论，发现好多不懂的）

多分类情况

其实结果跟前面呢二分类差不多，不过是扩展到了多维情况下。

假设存在N个类别，且第i类中的样本数表示为mi。

定义μ为所有样本的均值向量，如图中就是二维下的情况。

μi表示第i类的所有样本的均值向量；Sωi表示第i类的散度矩阵，表示第i类相对于这一类的中心μi的分散程度。

考虑所有类的情况，定义全局散度矩阵：

Sω=∑i=1NSωi

展开类似与二分类的类内散度矩阵：

Sω=∑x∈Xi(x−μi)(x−μi)T

实质是将所有类别的散度矩阵加在一起。

接下来考虑类间散度矩阵Sb，在二分类中，只考虑了两个均值点μ0和μ1的情况；现在在多分类情况下，考虑每个均值点μi与全局的均值点μ之间的距离。

由于每个类别的样本数量不同会对全局均值点μ产生影响：

μ=1N∑x=1N∑i=0N(∑n=0mixn)

注：共有N个类别，且第i类中的样本数为mi。

所以还要引入加权求和，每个类的权值为：mi∑Ni=0mi。由于J对倍数不敏感，所以可以把下面的总和去掉，直接使用mi表示权值。

写出类间散度矩阵Sb：

Sb=∑i=1Nmi(μi−μ)(μi−μ)T

与二分类时的步骤一样，求出投影后的Sb和Sω，步骤就不作赘述了：

Sb~=ωTSbωSω~=ωTSωω

好了，现在可以写出目标函数了：

我们希望：同类样例的投影点尽可能接近，即散度矩阵尽可能小，也即Sω~要尽可能小；异类样例的投影点尽可能远，也即类间距离尽可能大，Sb~尽可能大；

写出与二分类时一样的目标函数：

Sb~Sω~=ωTSbωωTSωω

由于我们得到的分子分母都是散列矩阵，要将矩阵变成实数，需要取行列式。又因为行列式的值实际上是矩阵特征值的积，一个特征值可以表示在该特征向量上的发散程度。因此我们使用行列式来计算（此处我感觉有点牵强，道理不是那么有说服力）。

J=∣∣Sb~∣∣∣∣Sω~∣∣=∣∣ωTSbω∣∣∣∣ωTSωω∣∣

现在又回到了求J的最大值的问题了，跟前面一样的步骤进行求解。

使用拉格朗日乘子法，得到特征方程：

Sbω=λSωω

强调内容接下来就是求矩阵的特征值的问题了，先求出矩阵S−1ωSb的特征值，之后取前N−1个特征向量，构成ω即可。（好吧，这又是个要填的“坑”）

参考资料：

1、《机器学习》周志华

2、线性判别分析（Linear Discriminant Analysis）（一）

写了快一下午了，满脑子数学公式，还有矩阵分析真的有必要补补了。

(~.~)