线性判别分析LDA

线性判别分析LDA

1. LDA的思想

  线性判别分析（Linear Discriminant Analysis，简称LDA）是一种经典的线性学习方法，在二分类问题上最早由[Fisher, 1936]提出，亦称“Fisher判别分析”。LDA和PCA（主成分分析）有一定类似之处。

  LDA的思想比较简单：给定训练集，设法将数据集投影到一条直线上，使得同类数据的投影点尽可能接近，异类数据的投影点尽可能远离，用一句话概括 投影后类内方差最小，类间方差最大。所以从这里可以看出，LDA是一种监督学习的方法。

  我们以下图为例，假设我们有两类数据，分别为红色和蓝色，这些数据特征是二维的，我们希望将这些数据投影到一维的一条直线，让每一种类别数据的投影点尽可能的接近，而红色和蓝色数据中心之间的距离尽可能的大。从直观上来看，我们会认为右图的投影方式更加合理，因为在当前方式下，每类数据彼此分离不重合，并且每一类的数据比较集中；相反在左图的投影方式下，效果可能不如人意。
  在更高维的特征空间中，我们投影后的不一定是直线，很可能是维数较低的超平面。

2. 瑞利商(Rayleigh quotient)与广义瑞利商

  我们先来看一下瑞利商函数的定义：
$$R(A,x)=\frac{x^{H}Ax}{x^{H}x}$$
其中 $x$ 是非零向量，而A为n阶厄米特矩阵。所谓厄米特矩阵，就是A矩阵的共轭转置等于它本身，即 $A^H=A$。而这个瑞利商函数 $R(A,x)$ 有一个非常重要的性质，就是它的最大值等于A的最大特征值，最小值等于A的最小特征值，即：
$${\lambda }_{min}\le R(A,x)\le {\lambda }_{max}$$

---

这个性质的证明如下:

  由于矩阵A是厄米特矩阵，所以一定存在一个酉矩阵$Q$，使得：
$$Q^{H}AQ=diag({\lambda }_1,{\lambda }_2,...,{\lambda }_n)=\Lambda$$
也就是说A可正交相似对角化。那么我们令$x=Qy$，上述式子变为：
$$\begin{array}{rl}R(A,x)& =\frac{y^H{Q}^{H}AQy}{y^H{Q}^{H}Qy}\\ & =\frac{y^H\Lambda y}{y^{H}y}\\ & =\frac{{\lambda }_1{y}_1^2+{\lambda }_2{y}_2^2+...+{\lambda }_n{y}_n^2}{y_1^2+y_2^2+...+y_n^2}\end{array}$$
  对此我们可以进行一个放缩：
$$\begin{array}{r}{\lambda }_{min}\le \frac{{\lambda }_1{y}_1^2+{\lambda }_2{y}_2^2+...+{\lambda }_n{y}_n^2}{y_1^2+y_2^2+...+y_n^2}\le {\lambda }_{max}\end{array}$$
  不妨假设${\lambda }_n$为最大的特征值，${\lambda }_1$为最小的特征值；那么当$y=(0,0,...,1)$时，$R(A,x)$可以取到最大值；当$y=(1,0,...,0)$ 时，$R(A,x)$ 取最小值。

  证毕。

这里我们还应该思考一下，就是当 $R(A,x)$ 取最大值时，$x$ 向量的取值是什么？

我们已经知道当 $y=(0,0,...,1)$时，$R(A,x)$取最大值，根据 $x=Qy$可以得到，$x$为 Q的最后一列；同理可以得到，当 $R(A,x)$ 等于某个特征值 ${\lambda }_i$时， $y=(0,0,...1,..,0,0)$，而 $x$ 为 Q的第 i 列；Q是什么？Q是 $A$ 的特征向量组成的矩阵。

知道这一点是很重要的，对于后面的LDA原理分析奠定基础。

---

当 $x$ 是标准正交基的时候，满足$x^{H}x=1$，这时$R(A,x)=x^{H}Ax$ 。

  下面我们来介绍一下广义瑞利商函数：
$$R(A,B,x)=\frac{x^{H}Ax}{x^{H}Bx}$$
其中 $x$ 是非零向量，而A、B为n阶厄米特矩阵。B为正定矩阵。那么它的最大值和最小值是什么呢？

  令 $x=B^{-1/2}{x}^{\prime }$，则分母变为：
$$\begin{array}{rl}{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 }\end{array}$$
  而分子变为：
$$\begin{array}{r}{x}^{H}Ax=x^{\prime H}{B}^{-1/2}A{B}^{-1/2}{x}^{\prime }\end{array}$$
  此时广义瑞利商函数$R(A,B,x)$变为$R(A,B,x^{\prime })$：
$$R(A,B,x)=\frac{x^{\prime H}{B}^{-1/2}A{B}^{-1/2}{x}^{\prime }}{x^{\prime H}{x}^{\prime }}$$

  利用前面得到瑞利商的性质，我们可以知道$R(A,B,x)$的最大值为矩阵$B^{-1/2}A{B}^{-1/2}$的最大特征值，最小值为矩阵$B^{-1/2}A{B}^{-1/2}$的最小特征值，而矩阵$B^{-1/2}A{B}^{-1/2}$和$B^{-1}A$相似——相似矩阵拥有相同的特征值，故广义瑞利商函数的最大值为矩阵$B^{-1}A$的最大特征值，最小值为$B^{-1}A$的最小特征值。

  同样的，我们更加关注广义瑞利商取到最大值时$x$ 的取值。根据前面瑞利商的性质，$R(A,B,x)$ 取到最大值时，$x^{\prime }$ 为矩阵 $B^{-1/2}A{B}^{1/2}$ 的特征向量，根据特征值特征向量的定义：
$$B^{-1/2}A{B}^{-1/2}\ast x^{\prime }=\lambda \ast x^{\prime }$$
将 $x=B^{-1/2}{x}^{\prime }$ 代入得：
$$\begin{array}{rl}{B}^{-1/2}A{B}^{-1/2}\ast B^{1/2}\ast x& =\lambda \ast B^{1/2}x\\ B^{-1}A\ast x& =\lambda \ast x\end{array}$$
综上可以看出，广义瑞利商的最大值为矩阵$B^{-1}A$的最大特征值，而 $x$的取值为其特征值对应的特征向量。

这里矩阵 B 肯定是半正定的，但不一定是正定的，即不一定可逆，常用的处理手段是：

$B+\rho I$，给矩阵 B 加上一个微小的扰动，这样矩阵从半正定变为正定 (因为 0特征值变为正数了)，之后就可以进行求逆的操作，而对原问题的影响比较小。

还有人可能会矩阵 $B$ 取 1/2次幂有些好奇，这里解释一下

由于矩阵 $B$ 是正定矩阵，故 $B$ 一定相似于主对角元素都为正数的对角阵，也就是说存在可逆阵 $P$，使得 $P^{-1}BP=\Lambda =diag({\lambda }_1,{\lambda }_2,...,{\lambda }_n)$ 是对角阵。令 $C=P^{-1}diag({\sqrt{\lambda }}_1,{\sqrt{\lambda }}_2,...,{\sqrt{\lambda }}_n)P$，可以看到 $C\ast C=B$，故 $C=B^{1/2}$

3. 二类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\in\{0,1\} yi​∈{0,1}，代表数据点的类别。接着定义$N_j(j=0,1)$为第j类样本的个数，$X_j(j=0,1)$代表第 j 类样本的集合，${\mu }_j(j=0,1)$代表第 j 类样本坐标的均值向量，${\sum }_j(j=0,1)$为第 j 类样本坐标的协方差矩阵。

  即：
$$\begin{gathered} {\mu }_j=\frac{1}{N_j}\sum _{x\in X_j}x(j=0,1) \\ \sum _j=\sum _{x\in X_j}(x-{\mu }_j)(x-{\mu }_j{)}^T(j=0,1) \end{gathered}$$
  因为我们要将两类数据投影到一条直线上，设直线向量为 $w$，故对于任意一个样本$x_i$，其在直线上的投影为$w^T{x}_i$，同理对于两类数据的中心点${\mu }_0,{\mu }_1$，它们在直线上的投影为$w^T{\mu }_0$和$w^T{\mu }_1$。而投影之后，类内数据点的方差变为：
$$\begin{array}{rl}\sum _{wj}& =\sum _j(w^{T}x-w^T{\mu }_j{)}^2\\ & =\sum _j(w^T(x-{\mu }_j){)}^2\\ & =\sum _j{w}^T(x-{\mu }_j)(x-{\mu }_j{)}^{T}w\\ & =w^T\sum _j(x-{\mu }_j)(x-{\mu }_j{)}^{T}w\\ & =w^T\sum _{j}w\end{array}$$
  LDA的目的是为了投影后类内方差最小，类间方差最大，所以我们要使$||w^T{\mu }_0-w^T{\mu }_1|{|}_2^2$最大，而最小化$w^T{\sum }_{0}w+w^T{\sum }_{1}w$。综上，我们需要优化：
$$J(w)=\frac{||w^T{\mu }_0-w^T{\mu }_1|{|}_2^2}{w^T{\sum }_{0}w+w^T{\sum }_{1}w}$$
使上述函数最大，也就达到了我们的目的。

  我们定义类内散度矩阵$S_w$：
$$S_w=\sum _0+\sum _1=\sum _{x\in X_0}(x-{\mu }_0)(x-{\mu }_0)+\sum _{x\in X_1}(x-{\mu }_1)(x-{\mu }_1)$$
  定义类间散度矩阵$S_b$：
$$S_b=({\mu }_0-{\mu }_1)({\mu }_0-{\mu }_1{)}^T$$
  这样上述优化目标$J(w)$可以写为：
$$J(w)=\frac{w^T{S}_{b}w}{w^T{S}_{w}w}$$
  而这不就是第二节我们讲过的广义瑞利商吗？利用前面瑞利商的性质，我们可以得到$J(w)$的最大值为矩阵$S_w^{-1}{S}_b$的最大值，最小值为矩阵$S_w^{-1}{S}_b$的最小值。但是我们关注的不仅仅是最值，而是这里的$w$的取值；根据前面广义瑞利商的性质，$w$ 为矩阵 $S_w^{-1}{S}_b$ 特征值对应的特征向量。

---

其他解法：

  可能有些人在西瓜书等地方看到过这样的解法:

  我们可以看到 $J(w)$ 的取值和 $w$ 的长度无关，只和其方向有关，因为可以通过约分手段约去长度，所以我们可以令 $w^T{S}_{w}w=1$。

  (可能有些人对这里会有些怀疑，可以回到普通瑞利商，可以令 $x^{H}x=1$，这一点是毫无疑问的；至于广义瑞利商，同样可以将其化为普通瑞利商，所以令 $w^T{S}_{w}w=1$ 并无大碍。)

  之后这个上述的优化目标可以重新表述为：
$$\begin{gathered} mi{n}_w-w^T{S}_{b}w \\ s.t.w^T{S}_{w}w=1 \end{gathered}$$
根据拉格朗日乘子法，上式等价于：（拉格朗日乘子没有正负限制，所以 $\lambda$ 前面可以是加号或者减号）
$$minL(w,\lambda )=-w^T{S}_{b}w+\lambda (w^T{S}_{w}w-1)$$
求导得：
$$\begin{array}{rl}{S}_{b}w& =\lambda S_{w}w\\ S_w^{-1}{S}_{b}w& =\lambda w\end{array}$$
到这里已经可以看到，这里 $w$ 就是矩阵 $S_w^{-1}{S}_b$的特征向量，代入原式得，最大值为最大特征值。

  多说些，观察一下 $S_{b}w=({\mu }_0-{\mu }_1)({\mu }_0-{\mu }_1{)}^{T}w$，会发现其方向为 ${\mu }_0-{\mu }_1$，设 $S_{b}w=k({\mu }_0-{\mu }_1)$，再将其带入上式得：$w=\frac{k}{\lambda }{S}_w^{-1}({\mu }_0-{\mu }_1)$，而由于 $w$ 只和方向有关，和长度无关，所以这里可以将其长度约去，最终得：$w=S_w^{-1}({\mu }_0-{\mu }_1)$；同样这也是西瓜书上 “令$S_{b}w=\lambda ({\mu }_0-{\mu }_1)$” 的原因。

---

4. 多类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 ∈ { 1 , 2 , . . . , k } y_i\in\{1,2,...,k\} yi​∈{1,2,...,k}，代表数据点的类别。接着定义$N_j(j=1,2,...,k)$为第j类样本的个数，$X_j(j=1,2,...,k)$代表第 j 类样本的集合，${\mu }_j(j=1,2,...,k)$代表第 j 类样本坐标的均值向量，${\sum }_j(j=1,2,...,k)$为第 j 类样本坐标的协方差矩阵。之后我们就可以将二类LDA推广到多类LDA。

  在多类中，我们投影到的低维空间很可能不是一条直线，而是一个超平面。我们假设超平面的维度为 $d$，对应的基向量为 $w_1,w_2,...,w_d$，构成矩阵 $W$。值得注意的是，$W$是一个 $n\ast d$ 的矩阵，这也是显而易见的，因为 $ W ^Tx_i$，即和原空间数据点做内积，将其映射到 $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$，$S_w={\sum }_{j=1}^k{\sum }_{x\in X_j}(x-{\mu }_j)(x-{\mu }_j{)}^T$代表类内散度矩阵。（这里的类间散度矩阵和二类LDA是有一些不同的，这里采用各类中心和所有数据点中心的距离来衡量的）

  和二类LDA不同的是，此时的 $J$ 不是一个标量，（在二类LDA中 $J$ 是函数），而是矩阵，这样的形式一般不好作为优化目标，所以我们可以采用矩阵对角线元素的乘积来进行代替，即：
$$\begin{array}{r}J(W)=\frac{\prod _{diag}{W}^T{S}_{b}W}{\prod _{diag}{W}_T{S}_{w}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}\end{array}$$
  然后仔细观察上面最右边的式子，会发现这不就是前面提到的广义瑞利商吗？（多个瑞利商的乘积）而广义瑞利商的最大值为矩阵$S_w^{-1}{S}_b$的最大特征值，故此时 $J(W)$的最大值就是矩阵$S_w^{-1}{S}_b$的前d个最大特征值的乘积，同时W就是特征值对应的特征向量组成的矩阵。

  虽说矩阵$W$的维度是 $n\ast d$ 的，代表将 $n$ 维空间数据点投影到d维空间；但是观察矩阵 $S_b$，其是由$k$个秩为1的矩阵相加得来的，并且 $k$个矩阵线性相关，即第 $k$ 个矩阵可以由前 $k-1$ 个矩阵线性表出，所以 $S_b$ 的秩小于等于 $k-1$ ，这是因为 $rank(A+B)\le rank(A)+rank(B)$；而矩阵 $S_w^{-1}{S}_b$ 的秩小于等于$k-1$，这是因为 r a n k ( A B ) ≤ m i n { r a n k ( A ) , r a n k ( B ) } rank(AB)\le min\{rank(A), rank(B)\} rank(AB)≤min{rank(A),rank(B)}；所以$S_w^{-1}{S}_b$最多含有 $k-1$ 个非零特征值，故 $d\le k-1$。

5. LDA算法流程

  经过前面对LDA原理的介绍，我们可以给出LDA的算法流程：

​ 输入： 1. 数据集 $D$；2.待降维的空间维数 $d$

​ 输出： 1.降维后的数据集 $D^{\prime }$

1. 计算类内散度矩阵 $S_w={\sum }_{j=1}^k{\sum }_{x\in X_j}(x-{\mu }_j)(x-{\mu }_j{)}^T$；
2. 计算类间散度矩阵 $S_b={\sum }_{j=1}^k{N}_j({\mu }_j-\mu )({\mu }_j-\mu {)}^T$；
3. 计算矩阵 $S_b^{-1}{S}_w$；
4. 计算 $S_b^{-1}{S}_w$ 的 d 个最大特征值及其对应的特征向量，组成 $W$；
5. 对原数据集 $D$ 中的数据点 $x_i$，进行投影，得到 $z_i=W^T{x}_i$:
6. 输出投影之后的数据集 $D^{\prime }$；

举例：

  数据集 $D$:
$$D=\left[\begin{array}{ccc}1& 2& 0\\ 3& 1& 0\\ -2& -2& 1\\ -3& -1& 1\end{array}\right]$$
  1. 类内散度矩阵 $S_w=\left[\begin{array}{cc}2.5& -1.5\\ -1.5& 1\end{array}\right]$

  2. 类间散度矩阵 $S_b=\left[\begin{array}{cc}20.25& 13.5\\ 13.5& 9.0\end{array}\right]$

  3. $S_b^{-1}{S}_w=\left[\begin{array}{cc}162& 108\\ 256.5& 171\end{array}\right]$

  4. ${\lambda }_1=0,{\lambda }_2=333$；$w_1=[-0.5547002,0.83205029],w_2=[-0.53399299,-0.8454889]$

   选择最大特征值对应的特征向量 $w=[-0.53399299,-0.8454889]$

  5. 计算投影之后的点集 $D^{\prime }=\left[\begin{array}{c}-2.21\\ -2.43\\ 2.74\\ 2.43\end{array}\right]$

6. LDA算法小结

  LDA降维的过程基本如上所述，同时也可以将LDA用于分类。不严谨地简述一下：将训练集投影到超平面或者直线之后，对待分类数据点继续投影得到新数据点，利用统计学手段判断其最可能属于某个类别。LDA虽说考虑了不同类别之间的关系，但是还是缺失了原数据集的某些信息，将其用于分类从直观上来说也许并不是很好的手段。

  下面再说一下LDA的优缺点：

  优点：

  (1) 计算速度快；

  (2) 利用数据点的类别信息，充分考虑先验知识；

  缺点：

  (1) 当数据不是高斯分布时候，效果不好，PCA也是；

  (2) 降维之后的维数最多为 $k-1$（k为类别数）；

7. 参考资料

1. 机器学习-LDA(线性判别降维算法)
2. 线性判别分析LDA原理总结——刘建平