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

LDA

  线性判别分析(Linear Discriminant Analysis, 简称LDA)，最早由Fisher提出，也叫“Fisher判别分析”。
  线性判别分析的思想：给定样本数据集，设法将样本投影到某一条直线上，使得同类样本的投影点尽可能接近，异类样本的投影点尽可能远；在对新的点进行分类预测时，将其投影到这条直线上，根据投影点的位置来判断样本的类别。
  当x是二维时，我们就要寻找一个方向为$\omega$的直线来使得这些样本点的投影分离。

二分类情况

  先只考虑二分类情况，给出数据集 $\left\{xi,yi\right\}_{i=1}^m, y_i=\left\{0,1\right\}$。
  设有$N_0$个样本的标签$y_i=0$，其类别为$D_0$；$N_1$个样本的标签$y_i=1$，其类别为$D_1$。
  则分别对应类别$D_0$和$D_1$的均值向量$\mu_0$和$\mu_1$可以求出来：

$$\mu_0=\frac{1}{N_0}\sum_{(x,y)\in{D_0}}{x_i}$$ $$\mu_1=\frac{1}{N_1}\sum_{(x,y)\in{D_1}}{x_i}$$
  再给出一组新的度量值，称作散度值(scatter)： $$s_0=\sum_{(x,y)\in{D_0}}{(x_i-\mu_0)^2}\\ s_1=\sum_{(x,y)\in{D_1}}{(x_i-\mu_1)^2}$$
  这东西看上去很眼熟吧，不就是方差，少除了样本数量吗？因为这里只需要定量表示样本集合的分散程度，对常数系数不敏感，所以方差中上面的部分就足够了。

  给出一组参数$\omega$，假设值为：$h(x_i)=\omega^T x_i$，后面简写为$h_i=\omega^T x_i$。
  在线性回归中，我们的目标是使得这个假设值等于样本的标签$y_i$；
而线性判别分析中，假设值$h_i$实质上是样本$x_i$在$w$上的投影的长度；以二维情况考虑，就是样本$x_i$对应的点$(x_1^i,x_2^i)$到直线$\omega$的投影点到原点的长度。我们的目标是通过区分这个”长度“$h_i$；
再放一次这幅图，看看图不难理解：

  接着，根据$\omega$可以求出投影后的样本均值$\tilde{\mu_0}$和$\tilde{\mu_1}$：

$$\tilde{\mu_0}=\frac{1}{N_0}\sum_{(x,y)\in{D_0}}{h_i}\\ =\frac{1}{N_0}\sum_{(x,y)\in{D_0}}{\omega^Tx_i}\\ =\omega^T\frac{1}{N_0}\sum_{(x,y)\in{D_0}}{x_i}\\ =\omega^T\mu_0$$
  同理得： $$\tilde{\mu_1}=\omega^T\mu_1$$

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

$$\tilde{s_0}=\frac{1}{N_0}\sum_{(x,y)\in{D_0}}{(h_i-\tilde{\mu_0})^2}\\ =\frac{1}{N_0}\sum_{(x,y)\in{D_0}}{(\omega^Tx_i-\omega^T\mu_0)^2}\\ =\frac{1}{N_0}\sum_{(x,y)\in{D_0}}{\omega^T(x_i-\mu_0)^2\omega}\\ =\omega^T \left[\frac{1}{N_0}\sum_{(x,y)\in{D_0}}{(x_i-\mu_0)^2}\right] \omega\\ =\omega^Ts_0\omega^T$$
  同理得： $$\tilde{s_1}=\omega^Ts_1\omega$$

  回到最开始说的思想：希望同类样例的投影点尽可能接近，即散度矩阵尽可能小，也即$\tilde{s_0}$和$\tilde{s_1}$都要尽可能小，简单表示为$\tilde{s_0}+\tilde{s_1}$；异类样例的投影点尽可能远，也即均值差值尽可能大，$(\tilde{\mu_0}-\tilde{\mu_1})^2$尽可能大；
  同时考虑上面两个条件，可以给出一个目标函数：

$$J=\frac{(\tilde{\mu_0}-\tilde{\mu_1})^2}{\tilde{s_0}+\tilde{s_1}}$$
  代入前面求出的 $\tilde{\mu_0}$、 $\tilde{\mu_1}$、 $\tilde{s_0}$、 $\tilde{s_1}$:
$$J=\frac{(\omega^T\mu_0-\omega^T\mu_1)^2}{\omega^Ts_0\omega+\omega^Ts_1\omega}\\ =\frac{\omega^T(\mu_0-\mu_1)^2\omega}{\omega^Ts_0\omega+\omega^Ts_1\omega}\\ =\frac{\omega^T(\mu_0-\mu_1)^T(\mu_0-\mu_1)\omega}{\omega^T(s_0+s_1)\omega}\\$$

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

$$S_\omega=S_0+S_1\\ =\sum_{(x,y)\in{D_0}}{(x_i-\mu_0)^2}+\sum_{(x,y)\in{D_1}}{(x_i-\mu_1)^2}\\ =\sum_{(x,y)\in{D_0}}{(x_i-\mu_0)^T(x_i-\mu_0)}+\sum_{(x,y)\in{D_1}}{(x_i-\mu_1)^T(x_i-\mu_1)}$$
  还有定义类间散度矩阵(between-class scatter matrix)： $$S_b=(\mu_0-\mu_1)^2\\ =(\mu_0-\mu_1)^T(\mu_0-\mu_1)$$
  则优化的目标函数变为： $$J=\frac{\omega^TS_b\omega}{\omega^TS_\omega\omega}$$
  这个东西就是LDA希望最大化的目标，实质是 $S_b$和 $S_\omega$的“广义瑞利商”(generalized Rayleigh quotient)。

  好的，接下来的任务就是通过$J$来确定最优的$\omega$。
  观察式子可以发现，$J$的大小不随$\omega$的大小变化，上下都是$\omega$的二次项相互抵消，结果而只与其方向有关。
  所以可以令$\omega^TS_\omega\omega=1$，则式子等价于：

$$\min_\omega -\omega^TS_b\omega\\ s.t. \quad \omega^TS_\omega\omega=1$$
  使用拉格朗日乘子法：
  设未知数 $\lambda$，写出拉格朗日函数：
$$L(\omega)=\omega^TS_b\omega-\lambda(\omega^TS_\omega\omega-1)$$
  对拉格朗日函数求导，且导数为0：
$$\frac{\mathrm{d} L(\omega)}{\mathrm{d} \omega}=2S_b\omega-2\lambda S_\omega\omega=0$$
  矩阵求导中 $\omega^TS_b\omega$可以简单看作是 $S_b\omega^2$。
  得到结果：
$$S_b\omega=\lambda S_\omega\omega$$
  这是一个典型的求矩阵特征值的问题。
  从前面的公式：
$$J=\frac{\omega^TS_b\omega}{\omega^TS_\omega\omega}$$
  我们知道 $\omega$的大小并不影响结果，而是方向才会影响结果。
  另外，由这个式子：
$$S_b\omega=(\mu_0-\mu_1)(\mu_0-\mu_1)^T\omega$$
  观察发现 $S_b\omega$的方向恒为 $\mu_0-\mu_1$，因为任取其中一个 $(\mu_0-\mu_1)^T$与 $\omega$相乘之后是常数，总会剩下一个 $\mu_0-\mu_1$。
  所以，设一个新的常量 $\lambda_\omega$，使得：
$$S_b\omega=(\mu_0-\mu_1)\lambda_\omega$$

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

$$S_b\omega=\lambda S_\omega\omega$$

  得到：

$$(\mu_0-\mu_1)\lambda_\omega=\lambda S_\omega\omega$$
  若 $\omega$可逆，则有：
$$\omega=\frac{\lambda_\omega}{\lambda}S_\omega^{-1}(\mu_0-\mu_1)$$
  注意，由于最后结果只与 $\omega$的方向有关，与其大小无关；而这个结果前面的常量 $\frac{\lambda_\omega}{\lambda}$可以舍去：
$$\omega=S_\omega^{-1}(\mu_0-\mu_1)$$
  这个就是最终的结果。
  前面已经推导出了 $S_\omega$、 $\mu_0$、 $\mu_1$，代入即可。我们只要有散度矩阵和均值即可求出最优的 $\omega$。

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

多分类情况

  其实结果跟前面呢二分类差不多，不过是扩展到了多维情况下。
  假设存在$N$个类别，且第$i$类中的样本数表示为$m_i$。
  定义$\mu$为所有样本的均值向量，如图中就是二维下的情况。

  $\mu_i$表示第$i$类的所有样本的均值向量；$S_{\omega_i}$表示第$i$类的散度矩阵，表示第$i$类相对于这一类的中心$\mu_i$的分散程度。
  考虑所有类的情况，定义全局散度矩阵：

$$S_\omega=\sum_{i=1}^NS_{\omega_i}$$
  展开类似与二分类的类内散度矩阵：
$$S_\omega=\sum_{x\in X_i}(x-\mu_i)(x-\mu_i)^T$$
  实质是将所有类别的散度矩阵加在一起。

  接下来考虑类间散度矩阵$S_b$，在二分类中，只考虑了两个均值点$\mu_0$和$\mu_1$的情况；现在在多分类情况下，考虑每个均值点$\mu_i$与全局的均值点$\mu$之间的距离。
  由于每个类别的样本数量不同会对全局均值点$\mu$产生影响：

$$\mu=\frac{1}{N}\sum x\\ =\frac{1}{N}\sum_{i=0}^N(\sum_{n=0}^{m_i}x_n)$$
  注：共有 $N$个类别，且第$i$类中的样本数为 $m_i$。
  所以还要引入加权求和，每个类的权值为： $\frac{m_i}{\sum_{i=0}^N m_i}$。由于 $J$对倍数不敏感，所以可以把下面的总和去掉，直接使用${m_i}$表示权值。
  写出类间散度矩阵 $S_b$：
$$S_b=\sum_{i=1}^N m_i(\mu_i-\mu)(\mu_i-\mu)^T$$
  与二分类时的步骤一样，求出投影后的 $S_b$和 $S_\omega$，步骤就不作赘述了：
$$\tilde{S_b}=\omega^T S_b \omega\\ \tilde{S_\omega} = \omega^T S_\omega \omega$$
  好了，现在可以写出目标函数了：
  我们希望：同类样例的投影点尽可能接近，即散度矩阵尽可能小，也即 $\tilde{S_\omega}$要尽可能小；异类样例的投影点尽可能远，也即类间距离尽可能大， $\tilde{S_b}$尽可能大；
  写出与二分类时一样的目标函数：
$$\frac{\tilde{S_b}}{\tilde{S_\omega}}=\frac{\omega^TS_b\omega}{\omega^TS_\omega\omega}$$
  由于我们得到的分子分母都是散列矩阵，要将矩阵变成实数，需要取行列式。又因为行列式的值实际上是矩阵特征值的积，一个特征值可以表示在该特征向量上的发散程度。因此我们使用行列式来计算（此处我感觉有点牵强，道理不是那么有说服力）。
$$J= \frac{\left|\tilde{S_b}\right|}{\left|\tilde{S_\omega}\right|}= \frac{\left|\omega^TS_b\omega\right|}{\left|\omega^TS_\omega\omega\right|}$$
  现在又回到了求 $J$的最大值的问题了，跟前面一样的步骤进行求解。
  使用拉格朗日乘子法，得到特征方程：
$$S_b\omega=\lambda S_\omega\omega$$
   强调内容接下来就是求矩阵的特征值的问题了，先求出矩阵 $S_\omega^{-1}S_b$的特征值，之后取前 $N-1$个特征向量，构成 $\omega$即可。（好吧，这又是个要填的“坑”）

参考资料：
1、《机器学习》周志华
2、线性判别分析（Linear Discriminant Analysis）（一）

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