线性判别分析(Linear Discriminant Analysis, LDA）(含类内散度矩阵类间散度矩阵全局散度矩阵推导

LDA算法概述：
线性判别式分析(Linear Discriminant Analysis, LDA)，也叫做Fisher线性判别(Fisher Linear Discriminant ,FLD)，是模式识别的经典算法，它是在1996年由Belhumeur引入模式识别和人工智能领域的。线性鉴别分析的基本思想是将高维的模式样本投影到最佳鉴别矢量空间，以达到抽取分类信息和压缩特征空间维数的效果，投影后保证模式样本在新的子空间有最大的类间距离和最小的类内距离，即模式在该空间中有最佳的可分离性。因此，它是一种有效的特征抽取方法。使用这种方法能够使投影后模式样本的类间散布矩阵最大，并且同时类内散布矩阵最小。就是说，它能够保证投影后模式样本在新的空间中有最小的类内距离和最大的类间距离，即模式在该空间中有最佳的可分离性。

说白了就是想让数据降维之后, 同一类的数据靠紧一点, 不同类的离远一点.

为什么要用LDA

既然PCA已经是一种表现很好的数据降维的方法，LDA也是一种将数据降维的方法。那为什么还要有LDA呢？
PCA是一种无监督的数据降维方法，与之不同的是LDA是一种有监督的数据降维方法。我们知道即使在训练样本上，我们提供了类别标签，在使用PCA模型的时候，我们是不利用类别标签的，而LDA在进行数据降维的时候是利用数据的类别标签提供的信息的。
从几何的角度来看，PCA和LDA都是讲数据投影到新的相互正交的坐标轴上。只不过在投影的过程中他们使用的约束是不同的，也可以说目标是不同的。PCA是将数据投影到方差最大的几个相互正交的方向上，以期待保留最多的样本信息。样本的方差越大表示样本的多样性越好，在训练模型的时候，我们当然希望数据的差别越大越好。否则即使样本很多但是他们彼此相似或者相同，提供的样本信息将相同，相当于只有很少的样本提供信息是有用的。样本信息不足将导致模型性能不够理想。这就是PCA降维的目标：将数据投影到方差最大的几个相互正交的方向上。这种约束有时候很有用，比如在下面这个例子：

对于这个样本集我们可以将数据投影到x轴或者y轴，但这都不是最佳的投影方向，因为这两个方向都不能最好地反映数据的分布。很明显还存在最佳的方向可以描述数据的分布趋势，那就是图中红色直线所在的方向。也是数据样本做投影，方差最大的方向。向这个方向做投影，投影后数据的方差最大，数据保留的信息最多。
但是，对于另外的一些不同分布的数据集，PCA的这个投影后方差最大的目标就不太合适了。比如对于下面图片中的数据集：

针对这个数据集，如果同样选择使用PCA，选择方差最大的方向作为投影方向，来对数据进行降维。那么PCA选出的最佳投影方向，将是图中红色直线所示的方向。这样做投影确实方差最大，但是是不是有其他问题。聪明的你一定发现了，这样做投影之后两类数据样本将混合在一起，将不再线性可分，甚至是不可分的。这对我们来说简直就是地狱，本来线性可分的样本被我们亲手变得不再可分。
帅气英俊的你也一定发现了，图中还有一条耀眼的黄色直线，向这条直线做投影即能使数据降维，同时还能保证两类数据仍然是线性可分的。上面的这个数据集如果使用LDA降维，找出的投影方向就是黄色直线所在的方向。
这其实就是LDA的思想，或者说LDA降维的目标：将带有标签的数据降维，投影到低维空间同时满足三个条件：
尽可能多地保留数据样本的信息（即选择最大的特征是对应的特征向量所代表的的方向）。
寻找使样本尽可能好分的最佳投影方向。
投影后使得同类样本尽可能近，不同类样本尽可能远。
其实第二个和第三个条件是基本等价的，我们很容易想到，使样本尽可能好分的投影方向，就是要使投影后使得同类样本尽可能近，不同类样本尽可能远。
上面大致讲解的LDA的基本思想，以及与PCA的不同，下面就来介绍一下LDA模型。

符号说明：

C:一个空间一共有C个类,C:一个空间一共有C个类,C:一个空间一共有C个类,

x∶表示训练样本，使用列向量表示x∶表示训练样本，使用列向量表示x∶表示训练样本，使用列向量表示

xj(i):第i个类的第j个样本,每个xj(i)是一个m行1列的矩阵，x_j^{(i)}:第i个类的第j个样本, 每个x_j^{(i)}是一个m行1列的矩阵，xj(i):第i个类的第j个样本,每个xj(i)是一个m行1列的矩阵，

Mi∶表示第i类训练样本的数目（i=1,2,…,C）M_i ∶表示第i类训练样本的数目（i=1,2,…,C）Mi∶表示第i类训练样本的数目（i=1,2,…,C）

M=∑i=1CMi∶表示训练样本的总数目M=\sum_{i=1}^CM_i ∶表示训练样本的总数目M=∑i=1CMi∶表示训练样本的总数目

μi=1Mi∑j=1MiXj(i)∶表示第i类训练样本的均值（i=1,2,…,C）\mu_i=\frac1{M_i}\sum_{j=1}^{M_i}X_j^{(i)}∶表示第i类训练样本的均值（i=1,2,…,C）μi=Mi1∑j=1MiXj(i)∶表示第i类训练样本的均值（i=1,2,…,C）

μ=1M∑i=1Mxi∶表示全部训练样本的均值\mu=\frac1{M}\sum_{i=1}^M x_i ∶表示全部训练样本的均值μ=M1∑i=1Mxi∶表示全部训练样本的均值

w∶投影直线w∶ 投影直线w∶投影直线

Σi=∑j=1Mi(xj(i)−μi)(xj(i)−μi)T∶表示i类样本的协方差矩阵\Sigma_i =\sum_{j=1}^{M_i}(x_j^{(i)}−\mu_i)(x_j^{(i)}−\mu_i)^T∶表示i类样本的协方差矩阵Σi=∑j=1Mi(xj(i)−μi)(xj(i)−μi)T∶表示i类样本的协方差矩阵

p(i)∶表示第i类的出现概率p(i) ∶表示第i类的出现概率p(i)∶表示第i类的出现概率

p(i,j)∶表示xj(i)的出现概率p(i,j) ∶表示x_j^{(i)}的出现概率p(i,j)∶表示xj(i)的出现概率

二类LDA原理:

由于是两类数据，因此我们只需要将数据投影到一条直线上即可。假设我们的投影直线是向量www,则对任意一个样本xix_ixi,它在直线www的投影为wTxiw^T x_iwTxi,对于我们的两个类别的中心点μ1μ_1μ1,μ2μ_2μ2,在直线www的投影为wTμ1w^T μ_1wTμ1 和wTμ2w^T μ_2wTμ2。由于LDA需要让不同类别的数据的类别中心之间的距离尽可能的大，也就是我们要最大化∣∣wTμ1−wTμ2∣∣22||w^T μ_1-w^T μ_2 ||_2^2∣∣wTμ1−wTμ2∣∣22,同时我们希望同一种类别数据的投影点尽可能的接近，也就是要同类样本投影点的协方差wTΣ1ww^T Σ_1 wwTΣ1w和wTΣ2ww^T Σ_2 wwTΣ2w尽可能的小，即最小化wTΣ1w+wTΣ2ww^T Σ_1 w+w^T Σ_2 wwTΣ1w+wTΣ2w。同时考虑二者，则可得到欲最大化的目标
J=∣∣wTμ1−wTμ2∣∣22wTΣ1w+wTΣ2w=wT(μ1−μ2)(μ1−μ2)TwwT(Σ1+Σ2)w(1)J=\frac{||w^T μ_1-w^T μ_2 ||_2^2}{w^T Σ_1 w+w^T Σ_2 w}=\frac{w^T (μ_1-μ_2 ) (μ_1-μ_2 )^T w}{w^T (Σ_1+Σ_2 )w} \tag1 J=wTΣ1w+wTΣ2w∣∣wTμ1−wTμ2∣∣22=wT(Σ1+Σ2)wwT(μ1−μ2)(μ1−μ2)Tw(1)

类内散度矩阵:

SW=Σ1+Σ2=∑i=1M1(xi(1)−μ1)(xi(1)−μ1)T+∑i=1M2(xi(2)−μ2)(xi(2)−μ2)T(2)S_W=Σ_1+Σ_2= ∑_{i=1}^{M_1}(x_i^{(1)}-μ_1 ) (x_i^{(1)}-μ_1)^T+ ∑_{i=1}^{M_2}(x_i^{(2)}-μ_2 ) (x_i^{(2)}-μ_2 )^T \tag2SW=Σ1+Σ2=i=1∑M1(xi(1)−μ1)(xi(1)−μ1)T+i=1∑M2(xi(2)−μ2)(xi(2)−μ2)T(2)

类间散度矩阵:

SB=(μ1−μ2)(μ1−μ2)T,(3)S_B=(μ_1-μ_2 ) (μ_1-μ_2 )^T,\tag3SB=(μ1−μ2)(μ1−μ2)T,(3)
则欲最大化的目标可重写为
J=wTSBwwTSWw(4)J=\frac{w^T S_B w}{w^T S_W w}\tag4J=wTSWwwTSBw(4)
不失一般性，令wTSWw=1w^T S_W w=1wTSWw=1，则优化目标等价于
max⁡w⁡wTSBws.t.wTSWw=1(5)\max_w⁡ \quad w^T S_B w \qquad s.t. \quad w^T S_W w=1 \tag5wmax⁡wTSBws.t.wTSWw=1(5)
由拉格朗日乘子法，上式等价于
SBw=λSWw,(6)S_B w=λS_W w,\tag6SBw=λSWw,(6)
其中λ是拉格朗日乘子. 注意到SBwS_B wSBw的方向恒为μ1−μ2μ_1-μ_2μ1−μ2, 不妨令
SBw=λ(μ1−μ2),(7)S_B w=λ(μ_1-μ_2 ),\tag7SBw=λ(μ1−μ2),(7)
代入式(6)中即得
w=SW−1(μ1−μ2).(8)w=S_W^{-1} (μ_1-μ_2).\tag8w=SW−1(μ1−μ2).(8)
考虑到数值解的稳定性，在实践中通常是对SWS_WSW进行奇异值分解，即SW=UΣVTS_W=UΣV^TSW=UΣVT，这里ΣΣΣ是一个实对称矩阵，其对角线上的元素是SWS_WSW的奇异值，然后再由SW−1=VΣ−1UTS_W^{-1}=VΣ^{-1}U^TSW−1=VΣ−1UT得到SW−1S_W^{-1}SW−1，或者是令SW=SW+γIS_W=S_W+γISW=SW+γI，γγγ是一个特别小的数，使得SWS_WSW可逆也是一种方法.

多类LDA原理:

在多类问题时，将DDD维的向量xxx投影到d<Dd<Dd<D维的yyy.投影方程为:
y=WTx(9)y=W^T x\tag9y=WTx(9)
这里WWW是一个投影矩阵，每一个列向量wiw_iwi表示一个投影方向.
设数据集合为X={x1(1)，x2(1)，…，xM1(1)，x1(2)，…xMC(C)}X=\lbrace x_1^{(1)}，x_2^{(1)}，…，x_{M_1}^{(1)}，x_1^{(2)}，…x_{M_C}^{(C)}\rbraceX={x1(1)，x2(1)，…，xM1(1)，x1(2)，…xMC(C)}. 在多类的时候，过程与上面一样，不过优化的目标函数需要更改一下.

类内散度矩阵与二类时的定义一样，如下表示:

SW=∑i=1C∑j=1Mip(i,j)(xj(i)−μi)(xj(i)−μi)T(10)S_W=∑_{i=1}^C∑_{j=1}^{M_i}p(i,j)(x_j^{(i)}-μ_i ) (x_j^{(i)}-μ_i )^T \tag{10}SW=i=1∑Cj=1∑Mip(i,j)(xj(i)−μi)(xj(i)−μi)T(10)

类间散度矩阵根据定义以及推导，得

SB=12∑i=1C∑j=1Cp(i)p(j)(μi−μj)(μi−μj)T=12∑i=1C∑j=1Cp(i)p(j)(μiμiT−μiμjT−μjμiT+μjμjT)=12∑i=1Cp(i)μiμiT∑j=1Cp(j)−p(i)μi∑j=1Cp(j)μjT−p(i)(∑j=1Cp(j)μj)μiT−p(i)∑j=1Cp(j)μjμjT=12∑i=1Cp(i)μiμiT−p(i)μiμT−p(i)μμiT−p(i)∑j=1Cp(j)μjμjT=12∑i=1Cp(i)μiμiT−μμT−μμT−∑j=1Cp(j)μjμjT=∑i=1Cp(i)(μiμiT−μμT)=∑i=1Cp(i)(μiμiT−μμT)+2∑i=1Cp(i)μμT−∑i=1Cp(i)μiμT−∑i=1Cp(i)μμiT=∑i=1Cp(i)(μiμiT−μiμT−μμiT+μμT)=∑i=1Cp(i)(μi−μ)(μi−μ)T\begin{aligned} S_B&=\frac12 ∑_{i=1}^C∑_{j=1}^Cp(i)p(j)(μ_i-μ_j ) (μ_i-μ_j )^T \\ &=\frac12 ∑_{i=1}^C∑_{j=1}^Cp(i)p(j)(μ_i μ_i^T-μ_i μ_j^T-μ_j μ_i^T+μ_j μ_j^T ) \\ &=\frac12 ∑_{i=1}^C{p(i) μ_i μ_i^T ∑_{j=1}^Cp(j)-p(i) μ_i ∑_{j=1}^Cp(j)μ_j^T- p(i)(∑_{j=1}^Cp(j)μ_j)μ_i^T- p(i) ∑_{j=1}^Cp(j)μ_j μ_j^T }\\ &=\frac12 ∑_{i=1}^C{p(i) μ_i μ_i^T-p(i) μ_i μ^T-p(i)μμ_i^T-p(i) ∑_{j=1}^Cp(j)μ_j μ_j^T }\\ &=\frac12 {∑_{i=1}^Cp(i) μ_i μ_i^T-μμ^T-μμ^T -∑_{j=1}^Cp(j)μ_j μ_j^T}\\ &=∑_{i=1}^Cp(i) (μ_i μ_i^T-μμ^T) \\ &=∑_{i=1}^Cp(i) (μ_i μ_i^T-μμ^T)+2∑_{i=1}^Cp(i)μμ^T-∑_{i=1}^Cp(i) μ_i μ^T- ∑_{i=1}^Cp(i)μμ_i^T \\ &=∑_{i=1}^Cp(i) (μ_i μ_i^T-μ_i μ^T-μμ_i^T+μμ^T)\\ &=∑_{i=1}^Cp(i)(μ_i-μ) (μ_i-μ)^T \end{aligned}SB=21i=1∑Cj=1∑Cp(i)p(j)(μi−μj)(μi−μj)T=21i=1∑Cj=1∑Cp(i)p(j)(μiμiT−μiμjT−μjμiT+μjμjT)=21i=1∑Cp(i)μiμiTj=1∑Cp(j)−p(i)μij=1∑Cp(j)μjT−p(i)(j=1∑Cp(j)μj)μiT−p(i)j=1∑Cp(j)μjμjT=21i=1∑Cp(i)μiμiT−p(i)μiμT−p(i)μμiT−p(i)j=1∑Cp(j)μjμjT=21i=1∑Cp(i)μiμiT−μμT−μμT−j=1∑Cp(j)μjμjT=i=1∑Cp(i)(μiμiT−μμT)=i=1∑Cp(i)(μiμiT−μμT)+2i=1∑Cp(i)μμT−i=1∑Cp(i)μiμT−i=1∑Cp(i)μμiT=i=1∑Cp(i)(μiμiT−μiμT−μμiT+μμT)=i=1∑Cp(i)(μi−μ)(μi−μ)T
即“类间散度矩阵”可表示为:
SB=∑i=1Cp(i)(μi−μ)(μi−μ)T(11)S_B=∑_{i=1}^Cp(i)(μ_i-μ) (μ_i-μ)^T \tag{11}SB=i=1∑Cp(i)(μi−μ)(μi−μ)T(11)
同时也可以定义“总体散度矩阵”STS_TST，表示为:
ST=∑i=1C∑j=1Mip(i,j)(xj(i)−μ)(xj(i)−μ)T(12)S_T=∑_{i=1}^C∑_{j=1}^{M_i}p(i,j)(x_j^{(i)}-μ) (x_j^{(i)}-μ)^T \tag{12}ST=i=1∑Cj=1∑Mip(i,j)(xj(i)−μ)(xj(i)−μ)T(12)
“类间散度矩阵”也可通过“总体散度矩阵”与“类内散度矩阵”的差值求出，即:
SB=ST−SW=∑i=1C∑j=1Mip(i,j)(xj(i)−μ)(xj(i)−μ)T−∑i=1C∑j=1Mip(i,j)(xj(i)−μi)(xj(i)−μi)T=∑i=1C∑j=1Mip(i,j){(xj(i)−μ)(xj(i)−μ)T−(xj(i)−μi)(xj(i)−μi)T}=∑i=1C{−p(i)μiμT−p(i)μμiT+p(i)μμT+p(i)μiμiT}=∑i=1Cp(i)(μi−μ)(μi−μ)T\begin{aligned} S_B&=S_T-S_W\\ &=∑_{i=1}^C∑_{j=1}^{M_i}p(i,j)(x_j^{(i)}-μ) (x_j^{(i)}-μ)^T -∑_{i=1}^C∑_{j=1}^{M_i}p(i,j)(x_j^{(i)}-μ_i ) (x_j^{(i)}-μ_i )^T \\ &=∑_{i=1}^C∑_{j=1}^{M_i}p(i,j)\lbrace(x_j^{(i)}-μ) (x_j^{(i)}-μ)^T-(x_j^{(i)}-μ_i ) (x_j^{(i)}-μ_i )^T \rbrace \\ &=∑_{i=1}^C\lbrace -p(i)μ_i μ^T- p(i)μμ_i^T+p(i)μμ^T+p(i)μ_i μ_i^T\rbrace \\ &=∑_{i=1}^Cp(i)(μ_i-μ) (μ_i-μ)^T \end{aligned}SB=ST−SW=i=1∑Cj=1∑Mip(i,j)(xj(i)−μ)(xj(i)−μ)T−i=1∑Cj=1∑Mip(i,j)(xj(i)−μi)(xj(i)−μi)T=i=1∑Cj=1∑Mip(i,j){(xj(i)−μ)(xj(i)−μ)T−(xj(i)−μi)(xj(i)−μi)T}=i=1∑C{−p(i)μiμT−p(i)μμiT+p(i)μμT+p(i)μiμiT}=i=1∑Cp(i)(μi−μ)(μi−μ)T
得到SBS_BSB，进一步可以得到目标函数:
J=∣WTSBW∣∣WTSWW∣(13)J=\frac{|W^T S_B W|}{|W^T S_W W|} \tag{13}J=∣WTSWW∣∣WTSBW∣(13)
由于SB，STS_B，S_TSB，ST不是标量，因此在目标函数中使用它们的行列式. 类似于二类推导的式子，可以得出
SW−1SBW=λW(14)S_W^{-1} S_B W=λW\tag{14}SW−1SBW=λW(14)
计算SW−1SBS_W^{-1} S_BSW−1SB的最大的ddd个特征值和对应的ddd个特征向量(w1,w2,…,wd)(w_1,w_2,…,w_d)(w1,w2,…,wd),即可得到投影矩阵W.

总结一下LDA的实现过程:

1)计算每个类样本的均值向量μiμ_iμi和所有样本的均值向量μμμ。
2)计算类内散度矩阵SWS_WSW和类间散度矩阵SBS_BSB。
3)计算SW−1SBS_W^{-1} S_BSW−1SB的特征向量W(w1,w2,…,wD)W(w_1,w_2,…,w_D)W(w1,w2,…,wD)和对应的特征值λ(λ1,λ2,…,λD)λ(λ_1,λ_2,…,λ_D)λ(λ1,λ2,…,λD)。
4)选择ddd个最大特征值对应的矩阵Wm∗dW_{m*d}Wm∗d。
5)对数据集XXX进行降维，得到对应的降维数据集Y=WTXY=W^T XY=WTX。

等价模型

根据LDA的基本思想，即让不同类别的数据的类别中心之间的距离尽可能的大，同时希望同一种类别数据的投影点尽可能的接近，可以得到不同的优化函数，如:
J(W)=∣WTSBW∣∣WTSWW∣J(W)=\frac{|W^T S_B W|}{|W^T S_W W|}J(W)=∣WTSWW∣∣WTSBW∣
J(W)=∣WTSBW∣∣WTSTW∣J(W)=\frac{|W^T S_B W|}{|W^T S_T W|}J(W)=∣WTSTW∣∣WTSBW∣
J(W)=∣WT(SB−SW)W∣∣WTW∣J(W)=\frac{|W^T (S_B-S_W )W|}{|W^T W|} J(W)=∣WTW∣∣WT(SB−SW)W∣
以上优化目标的目的都是一样的，为了最大化类间距离同时最小化类内距离，根据以上不同优化目标函数能够得到不同的LDA模型。

正则化

让模型泛化性能更高，防止模型过拟合。也可为了让SWS_WSW一定可逆而加。
J(W)=∣WTSBW∣∣WT(SW+λI)W∣J(W)=\frac{|W^T S_B W|}{|W^T (S_W+λI)W|}J(W)=∣WT(SW+λI)W∣∣WTSBW∣

比较fisherface和eigenface

LDA用于降维，和PCA有很多相同，也有很多不同的地方。

相同点：

两者均可以对数据进行降维。
两者在降维时均使用了矩阵特征分解的思想。

不同点：

LDA是有监督的降维方法，而PCA是无监督的降维方法
LDA降维最多降到类别数C-1的维数，而PCA没有这个限制。
LDA除了可以用于降维，还可以用于分类。
LDA选择分类性能最好的投影方向，而PCA选择样本点投影具有最大方差的方向。

参考

1.机器学习-线性判别分析.周志华
2.线性判别分析(LDA)