sklearn与机器学习系列专题之降维（二）一文弄懂LDA特征筛选&降维

1.PCA算法优缺点

在上一篇推文中，我们详解了PCA算法。这是机器学习中最为常用的降维方法，能降低算法的计算开销，使得数据集更容易处理，且完全无参数的限制。但是，如果用户对观测对象有一定的先验知识，掌握了数据的一些特征，却很难按照预想的方法对处理过程进行干预，可能达不到预期的效果，在非高斯分布的情况下，PCA方法得出的主元可能也并不是最优的。

2.LDA算法简介

这时候，就要线性判别分析降维（Linear discriminant Analysis，LDA）隆重登场啦！线性判别分析(也有叫做Fisher Linear Discriminant)，是一种有监督的线性降维算法，与PCA算法寻找数据集中方差最大的方向作为主成分分量的轴不同，LDA是为了使得降维后的数据点尽可能地容易被区分。

线性判别分析的中心思想是，最大化类间距离，最小化类内距离。

如下图，红色的点代表class1类别的数据，蓝色代表class2的数据，根据PCA算法，数据应该映射到方差最大的方向，即Y轴，但是class1和class2两个不同类别的数据就会完全的混合在一起，很难区分开。所以使用PCA算法进行降维后再进行分类的效果会非常差，这时候就需要我们使用LDA算法，将数据映射到X轴上。

3.枯燥又简洁的理论推导

由于理论推导可以通过网络、书籍等渠道了解到很多，此处以二分类的LDA为例，简要介绍原理。

上一节说道，LDA投影后，希望同一类别数据的投影点尽可能地靠近，而不同类数据的类别中心之间尽可能地大，因此，我们在推导过程中，将这一中心思想量化，如下，分子为不同类别之间中心点距离的差的平方，分母为同类样本之间的离散程度。最终目标是将J(W)最大化。

先将每一类中的散列值展开：

其中，x是投影前的特征，w为投影矩阵。
再对分子展开：

代入原式，目标函数最终可简化为：

其中SB代表类间散布矩阵，Sw代表类内散布矩阵。利用拉格朗日乘子法对上式变换，可得：

观察一下式子，可将它视为线性代数中的特征向量。因此，在线性判别分析中，只需要得到类内和类间散布矩阵，再求解其特征向量，即可得到投影方向，再对数据执行相应的矩阵变化，即可完成全部的降维过程。

4.python实战LDA

LDA在sklearn中有自带的包，导入方法为

from sklearn.discriminant_analysis import LinearDiscriminantAnalysis as LDA

使用自带库时，将上一篇推文代码中的PCA换成LDA即可，此处不再赘述代码过程。为了更清楚地向读者展示LDA算法的原理，以下将以“万能”的鸢尾花数据集为例，手撕LDA算法（选自《跟着迪哥学python》）！

# 自己来定义列名
feature_dict = {i:label for i,label in zip(
                range(4),
                  ('sepal length in cm',
'sepal width in cm',
'petal length in cm',
'petal width in cm', ))}

label_dict = {i:label for i,label in zip(
                range(1,4),
                  ('Setosa',
'Versicolor',
'Virginica'
                 ))}
import pandas as pd
# 数据读取，大家也可以先下载下来直接读取
df = pd.io.parsers.read_csv(
    filepath_or_buffer='https://archive.ics.uci.edu/ml/machine-learning-databases/iris/iris.data',
    header=None,
    sep=',',
    )
# 指定列名
df.columns = [l for i,l in sorted(feature_dict.items())] + ['class label']

df.head()

数据集共有150条数据，每条数据有4个特征，现在需要将四维特征降成二维。观察输出结果可以发现，其特征已经是数值数据，不需要做额外处理，但是需要转换一下标签：

from sklearn.preprocessing import LabelEncoder

X = df[['sepal length in cm','sepal width in cm','petal length in cm','petal width in cm']].values
y = df['class label'].values

# 制作标签{1: 'Setosa', 2: 'Versicolor', 3:'Virginica'}
enc = LabelEncoder()
label_encoder = enc.fit(y)
y = label_encoder.transform(y) + 1

上述代码使用了sklearn工具包中的“LabelEncoder”用于快速完成标签转换，可以发现基本上所有sklearn中的数据处理操作都是分两步走，先fit再transform。

在计算过程中需要基于均值来判断距离，因此先要对数据中各个特征求均值，但是只求4个特征的均值能满足要求吗？不要忘记任务中还有3种花，相当于3个类别，所以也要对每种花分别求其各个特征的均值。

import numpy as np
#设置小数点的位数
np.set_printoptions(precision=4)
#这里会保存所有的均值
mean_vectors = []
# 要计算3个类别
for cl in range(1,4):
# 求当前类别各个特征均值
    mean_vectors.append(np.mean(X[y==cl], axis=0))
    print('均值类别 %s: %s\n' %(cl, mean_vectors[cl-1]))

输出如下：

接下来计算类内散布矩阵

# 原始数据中有4个特征
S_W = np.zeros((4,4))
# 要考虑不同类别，自己算自己的
for cl,mv in zip(range(1,4), mean_vectors):
    class_sc_mat = np.zeros((4,4))
# 选中属于当前类别的数据
for row in X[y == cl]:
# 这里相当于对各个特征分别进行计算，用矩阵的形式
        row, mv = row.reshape(4,1), mv.reshape(4,1)
# 跟公式一样
        class_sc_mat += (row-mv).dot((row-mv).T)
    S_W += class_sc_mat
print('类内散布矩阵:\n', S_W)

输出为：

继续计算类间散布矩阵：

# 全局均值
overall_mean = np.mean(X, axis=0)
# 构建类间散布矩阵
S_B = np.zeros((4,4))
# 对各个类别进行计算
for i,mean_vec in enumerate(mean_vectors):
    #当前类别的样本数
n = X[y==i+1,:].shape[0]
mean_vec = mean_vec.reshape(4,1)
overall_mean = overall_mean.reshape(4,1)
    # 如上述公式进行计算
S_B += n * (mean_vec - overall_mean).dot((mean_vec - overall_mean).T)

print('类间散布矩阵:\n', S_B)

输出为：

得到类内和类间散布矩阵后，还需将它们组合在一起，然后求解矩阵的特征向量：

#求解矩阵特征值，特征向量
eig_vals, eig_vecs = np.linalg.eig(np.linalg.inv(S_W).dot(S_B))
# 拿到每一个特征值和其所对应的特征向量
for i in range(len(eig_vals)):
    eigvec_sc = eig_vecs[:,i].reshape(4,1)
    print('\n特征向量 {}: \n{}'.format(i+1, eigvec_sc.real))
    print('特征值 {:}: {:.2e}'.format(i+1, eig_vals[i].real))

输出特征向量：

输出结果得到4个特征值和其所对应的特征向量。特征向量直接观察起来比较麻烦，因为投影方向在高维上很难理解；特征值还是比较直观的，这里可以认为特征值代表的是其所对应特征向量的重要程度，也就是特征值越大，其所对应的特征向量就越重要，所以接下来可以对特征值按大小进行排序，排在前面的越重要，排在后面的就没那么重要了。

#特征值和特征向量配对
eig_pairs = [(np.abs(eig_vals[i]), eig_vecs[:,i]) for i in range(len(eig_vals))]

# 按特征值大小进行排序
eig_pairs = sorted(eig_pairs, key=lambda k: k[0], reverse=True)

print('特征值排序结果:\n')
for i in eig_pairs:
    print(i[0])

特征值排序结果为：

32.27195779972981
0.27756686384003953
1.1483362279322388e-14
3.422458920849769e-15

print('特征值占总体百分比:\n')
eigv_sum = sum(eig_vals)
for i,j in enumerate(eig_pairs):
print('特征值 {0:}: {1:.2%}'.format(i+1, (j[0]/eigv_sum).real))

特征值占总体百分比：

特征值 1: 99.15%
特征值 2: 0.85%
特征值 3: 0.00%
特征值 4: 0.00%

可以看出，打印出来的结果差异很大，第一个特征值占据总体的99.15%，第二个特征值只占0.85%，第三和第四个特征值，看起来微不足道。这表示对鸢尾花数据进行降维时，可以把特征数据降到二维甚至一维，但没必要降到三维。

既然已经有结论，选择把数据降到二维，只需选择特征值1、特征值2所对应的特征向量即可：

W = np.hstack((eig_pairs[0][1].reshape(4,1), eig_pairs[1][1].reshape(4,1)))
print('矩阵W:\n', W.real)

矩阵W:

[[-0.2049 -0.009 ]
[-0.3871 -0.589 ]
[ 0.5465 0.2543]
[ 0.7138 -0.767 ]]

这也是最终所需的投影方向，只需和原始数据组合，就可以得到降维结果：

# 执行降维操作
X_lda = X.dot(W)
X_lda.shape

现在可以看到数据维度从原始的（150,4）降到（150,2），到此就完成全部的降维工作。接下来对比分析一下降维后结果，为了方便可视化展示，在原始四维数据集中随机选择两维进行绘图展示。

from matplotlib import pyplot as plt
# 可视化展示
def plot_step_lda():

    ax = plt.subplot(111)
    for label,marker,color in zip(
        range(1,4),('^', 's', 'o'),('blue', 'red', 'green')):

        plt.scatter(x=X[:,0].real[y == label],
                y=X[:,1].real[y == label],
                marker=marker,
                color=color,
                alpha=0.5,
                label=label_dict[label]
                )

    plt.xlabel('X[0]')
    plt.ylabel('X[1]')

    leg = plt.legend(loc='upper right', fancybox=True)
    leg.get_frame().set_alpha(0.5)
    plt.title('Original data')

    # 把边边角角隐藏起来
    plt.tick_params(axis="both", which="both", bottom="off", top="off",
            labelbottom="on", left="off", right="off", labelleft="on")

    # 为了看的清晰些，尽量简洁
    ax.spines["top"].set_visible(False)
    ax.spines["right"].set_visible(False)
    ax.spines["bottom"].set_visible(False)
    ax.spines["left"].set_visible(False)

    plt.grid()
    plt.tight_layout
    plt.show()
    plot_step_lda()

从上述输出结果可以发现，如果对原始数据集随机取两维数据，数据集并不能按类别划分开，很多数据都堆叠在一起（尤其是图中方块和圆形数据点）。再来看看降维后的数据点分布，绘图代码保持不变，只需要传入降维后的两维数据即可。

from matplotlib import pyplot as plt
# 可视化展示
def plot_step_lda():

    ax = plt.subplot(111)
    for label,marker,color in zip(
        range(1,4),('^', 's', 'o'),('blue', 'red', 'green')):

        plt.scatter(x=X_lda[:,0].real[y == label],
                y=X_lda[:,1].real[y == label],
                marker=marker,
                color=color,
                alpha=0.5,
                label=label_dict[label]
                )

    plt.xlabel('LD1')
    plt.ylabel('LD2')

    leg = plt.legend(loc='upper right', fancybox=True)
    leg.get_frame().set_alpha(0.5)
    plt.title('LDA on iris')

    # 把边边角角隐藏起来
    plt.tick_params(axis="both", which="both", bottom="off", top="off",
            labelbottom="on", left="off", right="off", labelleft="on")

    # 为了看的清晰些，尽量简洁
    ax.spines["top"].set_visible(False)
    ax.spines["right"].set_visible(False)
    ax.spines["bottom"].set_visible(False)
    ax.spines["left"].set_visible(False)

    plt.grid()
    plt.tight_layout
    plt.show()

plot_step_lda()

可以明显看到，坐标轴变成LD1与LD2，这就是降维后的结果，从数据点的分布来看，混杂在一起的数据不多，划分起来就更容易。这就是经过一步步计算得到的最终降维结果。

在实验过程中，我们也可以发现，当拿到一份规模较大的数据集时，一方面可以通过观察特征值排序结果来决定，另一方面还需要通过实验交叉验证。

5.下篇预告

下一篇依然是以降维为主题，具体哪方面的，小编暂时保密，敬请期待！

欢迎关注公众号“码点联盟”，不定期分享机器学习、深度学习干货，及各种资料大礼包！