机器学习算法原理笔记（三）—— LDA线性判别分析

一、简述

线性判别分析(Linear discriminant Analysis，LDA)是一种监督学习的降维技术，与PCA不同的是，PCA是寻找数据集中方差最大的方向作为主成分分量的轴，而LDA是最优化分类的特征子空间。

LDA的思想可以用一句话概括，就是“投影后类内方差最小，类间方差最大”。也就是说，要将数据在低维度上进行投影，投影后希望每一种类别数据的投影点尽可能的接近，而不同类别的数据的类别中心之间的距离尽可能的大：

与PCA之间的对比：

可以看出，如果是PCA的话，为了方差最大化，会选择投影到左边，而LDA则会选择投影到下面。

二、二类LDA

假设我们的数据集 $D=\{(x_1,y_1), (x_2,y_2), ...,((x_m,y_m))\}$ ，其中任意样本xi为n维向量，yi∈{0,1}。我们定义 $N_{j}$ (j=0,1)为第j类样本的个数， $X_{j}$ (j=0,1)为第j类样本的集合，而 $μ_{j}$ (j=0,1)为第j类样本的均值向量，定义 $Σ_{j}$ (j=0,1)为第j类样本的协方差矩阵（严格说是缺少分母部分的协方差矩阵）。

μj的表达式为：

$\mu_j = \frac{1}{N_j}\sum\limits_{x \in X_j}x\;\;(j=0,1)$

Σj的表达式为：

$\Sigma_j = \sum\limits_{x \in X_j}(x-\mu_j)(x-\mu_j)^T\;\;(j=0,1)$

由于是两类数据，因此我们只需要将数据投影到一条直线上即可。假设我们的投影直线是向量w，则对任意一个样本本xi，它在直线w的投影为 $w^Tx_i$ ，我们希望同一种类别数据的投影点尽可能的接近，也就是要同类样本投影点的协方差 $w^T\Sigma_0w$ 和 $w^T\Sigma_1w$ 尽可能的小，即最小化下面的式子：

$w^T\Sigma_0w+w^T\Sigma_1w$

另一方面，我们希望让不同类别的数据尽可能地远离，考虑类的中心（也就是均值） $\mu_0,\mu_1$ ，只要使得 $||w^T\mu_0-w^T\mu_1||_2^2$ 尽可能大即可。

综合上面的考量，我们可以得出下面的优化目标：

$\underbrace{arg\;max}_w\;\;J(w) = \frac{||w^T\mu_0-w^T\mu_1||_2^2}{w^T\Sigma_0w+w^T\Sigma_1w} = \frac{w^T(\mu_0-\mu_1)(\mu_0-\mu_1)^Tw}{w^T(\Sigma_0+\Sigma_1)w}$

为了表示方便，我们做如下定义：

$S_w = \Sigma_0 + \Sigma_1 = \sum\limits_{x \in X_0}(x-\mu_0)(x-\mu_0)^T + \sum\limits_{x \in X_1}(x-\mu_1)(x-\mu_1)^T$

$S_b = (\mu_0-\mu_1)(\mu_0-\mu_1)^T$

优化目标重写为：

$\underbrace{arg\;max}_w\;\;J(w) = \frac{w^TS_bw}{w^TS_ww}$

利用凸优化的知识，我们使用和SVM中同样的技巧，约束 ${w^TS_ww}=1$ （和SVM中一样，上下式子一同缩放不影响最后求 $\underbrace{arg\;max}_w\;\;J(w)$ 的结果），这样使用拉格朗日乘子法就可以得到：

$J(w)=w^TS_bw-\lambda(w^TS_ww-1)$

$\frac{dJ}{dw}=2S_bw-2\lambda S_ww=0$

$S_bw=\lambda S_ww$

如果 $S_w$ 可逆，就又变回了求特征向量的套路：

$S_w^{-1}S_bw=\lambda w$

当然，在实际问题中， $S_w$ 常出现不可逆的问题，不过一般情况下可以通过其他方法解决：

令 $S_w=S_w+\gamma I$ ，其中 $\gamma$ 是一个特别小的数，这样 $S_w$ 一定可逆
先使用PCA对数据进行降维，使得在降维后的数据上 $S_w$ 可逆，再使用LDA

再往前一步：

我们考虑到两类的这种情况下，其实 $s_b w$ 和 $(\mu_1 - \mu_2)$ 是平行的，所以：

$S_b w= (\mu_1 - \mu_2) (\mu_1 - \mu_2)^T w = (\mu_1 - \mu_2) * k$

那么就算在式子两边左乘 $S_w^{-1}$ 也是成立的，那么原先的式子就可以转化为：

$S_w^{-1} S_b w = S_w^{-1} (\mu_1 - \mu_2) * k = w * \lambda$

两边同除以 k 对式子没有影响，所以，就可以得到得到特征向量：

$w = S_w^{-1}{(\mu_1 - \mu_2)}$

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

三、广义瑞利商求解二类LDA

瑞利商 $R(A,x)$ 的定义：

$R(A,x) = \frac{x^HAx}{x^Hx}$

其中，x为非零向量，而A为n×n的Hermitan矩阵(Hermitan矩阵就是满足共轭转置矩阵和自己相等的矩阵，即 $A^H=A$ ，如果矩阵A是实矩阵，则满足 $A^T = A$ 的矩阵即为Hermitan矩阵)

瑞利商 $R(A,x)$ 有一个非常重要的性质，即它的最大值等于矩阵A最大的特征值，而最小值等于矩阵A的最小的特征值，也就是满足(证明可以参考范数不等式):

$\lambda_{min} \leq \frac{x^HAx}{x^Hx} \leq \lambda_{max}$

广义瑞利商 $R(A,B,x)$ 的定义:

$R(A,x) = \frac{x^HAx}{x^HBx}$

其中，x为非零向量，而A,B为n×n的Hermitan矩阵，B为正定矩阵

广义瑞利商是不是很眼熟，我们完全可以套用广义瑞利商来求解之前的最大值问题，但先要得到广义瑞利商的最大最小值。

通过将广义瑞利商通过标准化就可以转化为瑞利商的形式，从而得到广义广义瑞利商的最大最小值：

令 $x=B^{-1/2}x'$ ，则：

分母： $x^HBx = x'^H(B^{-1/2})^HBB^{-1/2}x' = x'^HB^{-1/2}BB^{-1/2}x' = x'^Hx'$

分子： $x^HAx = x'^HB^{-1/2}AB^{-1/2}x'$

最后我们就可以得到：

$R(A,B,x') = \frac{x'^HB^{-1/2}AB^{-1/2}x'}{x'^Hx'}$

$R(A,B,x)$ 的最大值为矩阵 $B^{-1/2}AB^{-1/2}$ 的最大特征值，或者说矩阵 $B^{-1}A$ 的最大特征值，而最小值为矩阵 $B^{-1}A$ 的最小特征值。

四、多类LDA

假设我们的数据集 $D=\{(x_1,y_1), (x_2,y_2), ...,((x_m,y_m))\}$ ,其中任意样本 $x_i$ 为n维向量， $y_i$ ∈{C1,C2,...,Ck}。我们定义 $N_j$ (j=1,2...k)为第j类样本的个数， $X_j$ (j=1,2...k)为第j类样本的集合，而 $μ_j$ (j=1,2...k)为第j类样本的均值向量，定义 $Σ_j$ (j=1,2...k)为第j类样本的协方差矩阵。

假设我们投影到的低维空间的维度为d，对应的基向量为 $(w_1,w_2,...w_d)$ ，基向量组成的矩阵为W，它是一个n×d的矩阵，按我们在二类里的经验，此时我们的优化目标应该可以写成为:

$\frac{W^TS_bW}{W^TS_wW}$

$S_w$ 不难写出，毕竟二类里是两个，那多类里就是多弄几个而已：
$S_w = \sum\limits_{j=1}^{k}S_{wj} = \sum\limits_{j=1}^{k}\sum\limits_{x \in X_j}(x-\mu_j)(x-\mu_j)^T$

$S_b$ 应该怎么考虑呢，直接按之前的套路复杂度太高了，有没有什么可以代替的好办法呢？

我们先考虑整体的散度矩阵，也就是：

$S_t = \sum\limits_{x \in X_j}(x- \overline{x})(x- \overline{x})^T$

其中， $\overline{x}$ 是整体的均值

那么 $S_b$ 可以用 $S_t - S_w$ 来表示（下面的推导来自熊猫学写字的博客，符号有一些差异）

其中：

代入原式：

也就是：
$S_b = \sum\limits_{j=1}^{k}N_j(\mu_j-\mu)(\mu_j-\mu)^T$

现在我们的优化目标完全求解出来了，但是还有一个问题就是现在 $W^TS_bW$ 和 $W^TS_wW$ 都是矩阵，不是标量，无法作为一个标量函数来优化：

$\underbrace{arg\;max}_W\;\;J(W) = \frac{W^TS_bW}{W^TS_wW}$

一种解决方法是使用行列式的值来替换矩阵，解释是说实际上是矩阵特征值的积，一个特征值可以表示在该特征向量上的发散程度，因此我们可以使用行列式来计算。

也就是说我们的优化目标变为：

$\underbrace{arg\;max}_W\;\;J(W) = \frac{|W^TS_bW|}{|W^TS_wW|}$

这样使用之前的拉格朗日乘子法就可以轻松解决，结果也和上面是一样的。

不过需要注意的是由于 $S_b$ 的秩最大为k−1，所以 $S_2^{-1}S_b$ 的特征向量数目不会超过k−1，所以我们投影后的d<=(k−1)。(因为 $S_b$ 中每个 $μ_j − μ$ 的秩为1，因此协方差矩阵相加后最大的秩为k(矩阵的秩小于等于各个相加矩阵的秩的和)，但是由于如果我们知道前k-1个 $μ_j$ 后，最后一个 $μ_k$ 可以由前k-1个 $μ_j$ 线性表示，因此 $S_b$ 的秩最大为k-1，即特征向量最多有k-1个。

另一种解决方法是取矩阵对角线之和，也就是特征值之和：

$\underbrace{arg\;max}_W\;\;J(W) = \frac{\prod\limits_{diag}W^TS_bW}{\prod\limits_{diag}W^TS_wW}$

优化过程可以转化为：

$J(W) = \frac{\prod\limits_{i=1}^dw_i^TS_bw_i}{\prod\limits_{i=1}^dw_i^TS_ww_i} = \prod\limits_{i=1}^d\frac{w_i^TS_bw_i}{w_i^TS_ww_i}$

那么根据广义广义瑞利商即可求得和上面一样的结果。

五、LDA算法流程

输入：数据集 $D=\{(x_1,y_1), (x_2,y_2), ...,((x_m,y_m))\}$ ，其中任意样本 $x_i$ 为n维向 $y_i \in \{C_1,C_2,...,C_k\}$ ，降维到的维度d:

计算 $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^Tx_i$
得到降维后的数据

六、LDA和PCA

LDA的优缺点：

LDA和PCA的比较：

七、sklearn实现

class sklearn.discriminant_analysis.LinearDiscriminantAnalysis(solver='svd', shrinkage=None, priors=None, n_components=None, store_covariance=False, tol=0.0001)

Parameters:

solver: str类型，默认值为‘svd’
svd：使用奇异值分解求解，不用计算协方差矩阵，适用于特征数量很大的情形，无法使用参数收缩（shrinkage）
lsqr：最小平方QR分解，可以结合shrinkage使用
eigen：特征值分解，可以结合shrinkage使用
shrinkage: str or float类型，默认值为None
是否使用参数收缩
None：不使用参数收缩
auto：str，使用Ledoit-Wolf lemma
浮点数：自定义收缩比例
components：int类型，需要保留的特征个数，小于等于n-1

Attributes:

covariances_：每个类的协方差矩阵， shape = [n_features, n_features]
means_：类均值，shape = [n_classes, n_features]
priors_：归一化的先验概率
rotations_：LDA分析得到的主轴，shape [n_features, n_component]
scalings_：数组列表，每个高斯分布的方差σ

import numpy as np
from sklearn.discriminant_analysis import LinearDiscriminantAnalysis

X = np.array([[-1, -1], [-2, -1], [-3, -2], [1, 1], [2, 1], [3, 2]])
y = np.array([1, 1, 1, 2, 2, 2])

clf = LinearDiscriminantAnalysis()
clf.fit(X, y)

LinearDiscriminantAnalysis(n_components=None, priors=None, shrinkage=None,
              solver='svd', store_covariance=False, tol=0.0001)

print(clf.predict([[-0.8, -1]]))

>> [1]

参考：

最后编辑于：2020.01.11 13:14:53

人面猴
序言：七十年代末，一起剥皮案震惊了整个滨河市，随后出现的几起案子，更是在滨河造成了极大的恐慌，老刑警刘岩，带你破解...
沈念sama阅读 198,932评论 5赞 466
死咒
序言：滨河连续发生了三起死亡事件，死亡现场离奇诡异，居然都是意外死亡，警方通过查阅死者的电脑和手机，发现死者居然都...
沈念sama阅读 83,554评论 2赞 375
救了他两次的神仙让他今天三更去死
文/潘晓璐我一进店门，熙熙楼的掌柜王于贵愁眉苦脸地迎上来，“玉大人，你说我怎么就摊上这事。” “怎么了？”我有些...
开封第一讲书人阅读 145,894评论 0赞 328
道士缉凶录：失踪的卖姜人
文/不坏的土叔我叫张陵，是天一观的道长。经常有香客问我，道长，这世上最难降的妖魔是什么？我笑而不...
开封第一讲书人阅读 53,442评论 1赞 268
港岛之恋（遗憾婚礼）
正文为了忘掉前任，我火速办了婚礼，结果婚礼上，老公的妹妹穿的比我还像新娘。我一直安慰自己，他们只是感情好，可当我...
茶点故事阅读 62,347评论 5赞 359
恶毒庶女顶嫁案：这布局不是一般人想出来的
文/花漫我一把揭开白布。她就那样静静地躺着，像睡着了一般。火红的嫁衣衬着肌肤如雪。梳的纹丝不乱的头发上，一...
开封第一讲书人阅读 47,899评论 1赞 275
城市分裂传说
那天，我揣着相机与录音，去河边找鬼。笑死，一个胖子当着我的面吹牛，可吹牛的内容都是我干的。我是一名探鬼主播，决...
沈念sama阅读 37,325评论 3赞 390
双鸳鸯连环套：你想象不到人心有多黑
文/苍兰香墨我猛地睁开眼，长吁一口气：“原来是场噩梦啊……” “哼！你这毒妇竟也来了？” 一声冷哼从身侧响起，我...
开封第一讲书人阅读 35,980评论 0赞 254
万荣杀人案实录
序言：老挝万荣一对情侣失踪，失踪者是张志新（化名）和其女友刘颖，没想到半个月后，有当地人在树林里发现了一具尸体，经...
沈念sama阅读 40,196评论 1赞 294
护林员之死
正文独居荒郊野岭守林人离奇死亡，尸身上长有42处带血的脓包…… 初始之章·张勋以下内容为张勋视角年9月15日...
茶点故事阅读 35,163评论 2赞 317
白月光启示录
正文我和宋清朗相恋三年，在试婚纱的时候发现自己被绿了。大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
茶点故事阅读 37,085评论 1赞 328
活死人
序言：一个原本活蹦乱跳的男人离奇死亡，死状恐怖，灵堂内的尸体忽然破棺而出，到底是诈尸还是另有隐情，我是刑警宁泽，带...
沈念sama阅读 32,826评论 3赞 316
日本核电站爆炸内幕
正文年R本政府宣布，位于F岛的核电站，受9级特大地震影响，放射性物质发生泄漏。R本人自食恶果不足惜，却给世界环境...
茶点故事阅读 38,389评论 3赞 302
男人毒药：我在死后第九天来索命
文/蒙蒙一、第九天我趴在偏房一处隐蔽的房顶上张望。院中可真热闹，春花似锦、人声如沸。这庄子的主人今日做“春日...
开封第一讲书人阅读 29,501评论 0赞 19
一桩弑父案，背后竟有这般阴谋
文/苍兰香墨我抬头看了看天上的太阳。三九已至，却和暖如春，着一层夹袄步出监牢的瞬间，已是汗流浃背。一阵脚步声响...
开封第一讲书人阅读 30,753评论 1赞 255
情欲美人皮
我被黑心中介骗来泰国打工，没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留，地道东北人。一个月前我还...
沈念sama阅读 42,171评论 2赞 344
代替公主和亲
正文我出身青楼，却偏偏与公主长得像，于是被迫代替她去往敌国和亲。传闻我的和亲对象是个残疾皇子，可洞房花烛夜当晚...
茶点故事阅读 41,616评论 2赞 339