Python手撸机器学习系列（三）：线性判别分析（Fisher究极详解，看不懂你锤死我）

1.原理概述

我们的目的是将高维的数据投影到一维直线上并在投影的值中取一个阈值进行分类，如下图所示：（绘画水平有限，将就着看）

在上图，很明显左边的投影更适合分类，因为两种类别（o和x）在投影直线上能轻松地找到一个阈值将其区分开来，而右边的投影方向则不适合当前分类。

所以我们需要求解一个适合的投影方向$w$

在理解fisher的时候，我遇到了很多不理解问题，在经过多本书籍的对比之后终于搞懂了，其大致的思路如下：

1. 问题的初衷在于找到一条线将坐标点向该线上投影，将这条线的方向设为$w$，并用该$w$作为假设带入，最后解出最佳$w$
2. 按照我们假设的$w$，将样本点向该直线中投影，即$w^{T}x$，求出每一类样本点在投影上的均值和方差（或者说是协方差矩阵）
3. 按照类间小，类内大的目标，设立目标函数求解$w$

值得注意的是，我们求得的w是最终投影的平面（在这里为一维直线）方向，而不是感知机或逻辑斯蒂回归中的决策边界，这个问题一度让我以为自己的w求错了！！！

另外，在推导公式中，一定一定记得随时查看当前数据维度，不然极易混淆

2.原理详解

2.1 输入

对于样本点$X=(x_1,x_2,...,x_N{)}^T$，其维度为$N\times p$，即每一个样本有$p$个特征 ，其类别$Y=(y_1,y_2,...,y_N{)}^T$，其中 y i ∈ { + 1 , − 1 } y_i\in\{+1,-1\} yi​∈{+1,−1}

我们将样本按其标签分为两类，数量分别为$N_1$和$N_2$，即$|x_{c1}|=N_1,|x_{c2}|=N_2,N_1+N_2=N$

在本文中，我们用坐标点投影在直线上来直观地表述，即$p=2$

2.2 投影

假设，我是说假设啊，我们现在找到了一个适合的投影方向，则样本$x$在直线$z$上的投影为：$z_i=w^T{x}_i$这个公式的来源：我们假设$w$的模为1（因为方向重要长度不重要，可以缩放）。在两个向量中，你可以很明显地看出$x_i$在$w$上的投影为$\Delta =x_{i}cos\theta$，而$x_i$和$w$的点积为$x_i\cdot w=|x_i||w|cos\theta =|x_i|cos\theta =\Delta$，就是说我们可以用$x_i$和$w$的点积表示其$x_i$在$w$上的投影，写作$w^T{x}_i$的形式。

值得注意的是，在这里投影指的是 ：以原点为起点，$x_i$为终点的向量在直线$z$（方向为$w$）上的投影

然后，所有样本在$z$上的投影的均值为：$\bar{z}=\frac{1}{N}\sum _{i=1}^N{z}_i=\frac{1}{N}\sum _{i=1}^N{w}^T{x}_i$

协方差（在这里可以理解为方差）为$S_z=\frac{1}{N}\sum _{i=1}^N(z_i-\bar{z})(z_i-\bar{z}{)}^T=\frac{1}{N}\sum _{i=1}^N(w^T{x}_i-\bar{z})(w^T{x}_i-\bar{z}{)}^T$

两个类分开写：

类别1：$\overline{z_1}=\frac{1}{N_1}\sum _{i=1}^{N_1}{w}^T{x}_i$ ，$S_1=\frac{1}{N_1}\sum _{i=1}^{N_1}(w^T{x}_i-\overline{z_1})(w^T{x}_i-\overline{z_1}{)}^T$

类别2：$\overline{z_2}=\frac{1}{N_2}\sum _{i=1}^{N_2}{w}^T{x}_i$ ，$S_1=\frac{1}{N_2}\sum _{i=1}^{N_2}(w^T{x}_i-\overline{z_2})(w^T{x}_i-\overline{z_2}{)}^T$

注意，这里由于我们是平面坐标点在一维直线上的投影，所以$\bar{z}、S_z、\overline{z_1}、\overline{z_2}、S_1、S_2$均可以理解为一个数。

还记得我们的目标吗：类内小，类间大

在这里，我们将类间表示为$(\overline{z_1}-\overline{z_2}{)}^2$，即两个类别的样本分别取均值，其均值之差的平方

将类内表示为$S_1+S_2$，即两个类别的样本方差之和

由此我们可以很自然地得到一个最大化的目标函数：
$$J(w)=\frac{(\overline{z_1}-\overline{z_2}{)}^2}{S_1+S_2}$$
极大化这个式子就相当于最小化$S_1+S_2$和最大化$(\overline{z_1}-\overline{z_2}{)}^2$，完美符合我们的目标

对于$J(w)$，我们可以进一步的调整：

对于其分子:
$$\begin{array}{rl}(\overline{z_1}-\overline{z_2}{)}^2& ={\left[\frac{1}{N_1}\sum _{i=1}^{N_1}{w}^T{x}_i-\frac{1}{N_2}\sum _{i=1}^{N_2}{w}^T{x}_i\right]}^2\\ & ={\left[w^T(\frac{1}{N_1}\sum _{i=1}^{N_1}{x}_i-\frac{1}{N_2}\sum _{i=1}^{N_2}{x}_i)\right]}^2\\ & ={\left[w^T(\overline{x_{c1}}-\overline{x_{c2}})\right]}^2\\ & =w^T(\overline{x_{c1}}-\overline{x_{c2}})(\overline{x_{c1}}-\overline{x_{c2}}{)}^{T}w\end{array}$$
对于其分母，我们先取$S1$进行分析：
$$\begin{array}{rl}{S}_1& =\frac{1}{N_1}\sum _{i=1}^{N_1}(w^T{x}_i-\overline{z_1})(w^T{x}_i-\overline{z_1}{)}^T\\ & =\frac{1}{N_1}\sum _{i=1}^{N_1}{w}^T(x_i-\overline{x_{c1}})(x_i-\overline{x_{c1}}{)}^{T}w\\ & =w^T\frac{1}{N_1}\sum _{i=1}^{N_1}(x_i-\overline{x_{c1}})(x_i-\overline{x_{c1}}{)}^{T}w\\ & =w^T{S}_{c1}w\end{array}$$
上式中我们用$S_{c1}$表示$\frac{1}{N_1}\sum _{i=1}^{N_1}(x_i-\overline{x_{c1}})(x_i-\overline{x_{c1}}{)}^T$，很显然，这个式子的意思为类别1中样本的方差

值得注意的是，这里的$S_{c1}$与之前的$S_1$有很大的不同！！！！！！$S_1$表示的是投影的方差，维度1×1，在投影中可以理解为一个数！一个数！一个数！！！！而$S_{c1}$表示类别1样本的协方差矩阵，没有投影之前的协方差矩阵，维度为p×p，是一个矩阵！一个矩阵！一个矩阵！

这个地方一定要搞清楚，不然你就不知道为什么前面是方差是一个数但是后面方差求出来是一个矩阵了，我当时推了几遍才反应过来！

此外，$p\times p$的协方差矩阵为$x_i$的每一个特征的标准差与其他特征标准差的乘积

好了，书接上回：

$S_2$同理，所以其分母等价于：
$$S_1+S_2=w^T{S}_{c1}w+w^T{S}_{c2}w=w^T(S_{c1}+S_{c2})w$$

则$J(w)$可以转化为：
$$J(w)=\frac{w^T(\overline{x_{c1}}-\overline{x_{c2}})(\overline{x_{c1}}-\overline{x_{c2}}{)}^{T}w}{w^T(S_{c1}+S_{c2})w}=\frac{w^T{S}_{b}w}{w^T{S}_{w}w}$$
我们用$S_b$表示$(\overline{x_{c1}}-\overline{x_{c2}})(\overline{x_{c1}}-\overline{x_{c2}}{)}^T$，其被称为类间方差

用$S_w$表示$S_{c1}+S_{c2}$，其被称为类内方差

这里将$J(w)$转化的意义在于，将原来假设的投影$w$的方差和均值转化为了各类别样本本身的方差和均值，这样求解出来的w就可以直接用样本本来的值求解了，你甚至可以理解成把虚的变成实的！

2.3 求解$w$

由上面的式子，可以得到：

$J(w)=\frac{w^T{S}_{b}w}{w^T{S}_{w}w}=w^T{S}_{b}w\cdot (w^T{S}_{w}w{)}^{-1}$

对$w$进行求导并令结果为0：

$$\frac{\partial J(w)}{\partial w}=2S_{b}w\cdot (w^T{S}_{w}w{)}^{-1}+w^T{S}_{b}w\cdot (-1)(w^T{S}_{w}w{)}^{-2}\cdot 2S_{w}w=0$$

两边同时乘上$(w^T{S}_{w}w{)}^2$，注意这里$(w^T{S}_{w}w)$是一个数，维度计算：$1\times p\cdot p\times p\cdot p\times 1=1\times 1$

可得：
$$\begin{gathered} S_{b}w(w^T{S}_{w}w)-w^{T}Sbw\cdot S_{w}w=0 \\ {w}^{T}Sbw\cdot S_{w}w=S_{b}w(w^T{S}_{w}w) \end{gathered}$$
同$(w^T{S}_{w}w)$一样，$(w^T{S}_{b}w)$依然是一个数，所以可得：
$$S_{w}w=\frac{w^T{S}_{w}w}{w^T{S}_{b}w}\cdot S_b\cdot w$$
两边左乘$S_w^{-1}$，可得：
$$w=\frac{w^T{S}_{w}w}{w^T{S}_{b}w}\cdot S_w^{-1}\cdot S_b\cdot w$$
由于$\frac{w^T{S}_{w}w}{w^T{S}_{b}w}$是一个数，并不会影响我们求$w$的方向，所以我们不妨设它为1（这里可以理解为我们投影的直线$z$可以放缩），即可得：
$$w=S_w^{-1}\cdot S_b\cdot w$$
我们继续研究，可以发现$S_b\cdot w=(\overline{x_{c1}}-\overline{x_{c2}})(\overline{x_{c1}}-\overline{x_{c2}}{)}^T\cdot w$，而$(\overline{x_{c1}}-\overline{x_{c2}}{)}^T\cdot w$的维度计算为$1\times p\cdot p\times 1=1\times 1$，即实数，与$w$方向无关。

所以最终我们求得的$w$为：
$$w=S_w^{-1}(\overline{x_{c1}}-\overline{x_{c2}})$$
即两个类协方差矩阵之和的逆矩阵乘上两个类的均值之差

3.代码实现

有了上述结果，我们可以很清楚的得到算法流程：

1. 根据给定的两个类的样本，分别计算其协方差矩阵和样本均值
2. 将协方差矩阵相加得到$S_w$，将样本均值相减得到$(\overline{x_{c1}}-\overline{x_{c2}})$
3. 计算$w=S_w^{-1}(\overline{x_{c1}}-\overline{x_{c2}})$
4. 将新样本根据计算出的$w$投影，并进行分类

以下为python实现代码：

import numpy as np
import matplotlib.pyplot as plt

def fisher(x_1 , y_1 , x_2 ,y_2):

    #类内方差 Sw = S1 + S2
    u_1 = np.mean(x_1 , axis = 0 )
    S1 = (x_1 - u_1).T.dot((x_1-u_1))

    u_2 = np.mean(x_2 , axis = 0 )
    S2 = (x_2 - u_2).T.dot(x_2-u_2)

    Sw = S1 + S2
    # print(Sw)
    w = np.linalg.pinv(Sw).dot(u_1 - u_2) #inv逆矩阵，pinv伪逆矩阵

    return w

def predict(test_data , w , c1 , c2): #根据计算的W进行预测
    u_1 = np.mean(c1, axis=0)
    u_2 = np.mean(c2, axis=0)

    diff_1 = w.dot(u_1.T)
    diff_2 = w.dot(u_2.T)

    diff_cur = w.dot(test_data.T) #根据投影距离远近决定是哪一类

    return [1 if abs(diff_cur[i]-diff_1) < abs(diff_cur[i]-diff_2) else -1 for i in range(len(test_data))]

if __name__ == '__main__':
    X_True = [[1,1],[2,2],[0,4],[3,4]] #正例
    X_False = [[3,3],[4,5],[3,4],[5,4]] #负例
    Y_True = [1] * len(X_True)
    Y_False = [-1] * len(X_False)
    w = fisher(X_True,Y_True,X_False,Y_False) #求投影直线的方向w

    test_point = np.array([[2,5],[4,0]]) #测试数据
    predict_result = predict(test_point,w,X_True,X_False)
    print('测试点预测类别：'+str(predict_result))
    
    #绘图
    plot_x = np.arange(0,6)
    plot_y = -(w[0]*plot_x)/w[1]
    plt.scatter([x[0] for x in X_True] ,[x[1] for x in X_True] ,c = 'r',label = 'class 1')
    plt.scatter([x[0] for x in X_False] ,[x[1] for x in X_False] ,c = 'b',label = 'class -1')
    plt.scatter([x[0] for x in test_point] ,[x[1] for x in test_point] ,c = 'green',label = 'test_point')
    plt.plot(plot_x,plot_y,c = 'pink')
    plt.legend()
    plt.show()

运行结果：

图中的粉色直线即为所求结果，注意这个是投影的直线，而不是决策边界。

4.联系方式

如有问题欢迎评论区指出，有问题也可以联系我：
1759412770@qq.com