五步实现LDA
完成以上几项准备工作后,我们就可以实际运行LDA了。
第一步:计算数据的 d 维均值向量
首先做一个简单的计算:分别求三种鸢尾花数据在不同特征维度上的均值向量 mi:
np.set_printoptions(precision=4)mean_vectors = []for cl in range(1,4): mean_vectors.append(np.mean(X[y==cl], axis=0)) print('Mean Vector class %s: %sn' %(cl, mean_vectors[cl-1]))
得到
Mean Vector class 1: [ 5.006 3.418 1.464 0.244]Mean Vector class 2: [ 5.936 2.77 4.26 1.326]Mean Vector class 3: [ 6.588 2.974 5.552 2.026]第二步:计算散布矩阵
计算两个 4×4 维矩阵:类内散布矩阵和类间散布矩阵。
2.1a 类内散布矩阵 Sw
类内散布矩阵 Sw 由该式算出:
S_W = np.zeros((4,4))for cl,mv in zip(range(1,4), mean_vectors): class_sc_mat = np.zeros((4,4)) # scatter matrix for every class for row in X[y == cl]: row, mv = row.reshape(4,1), mv.reshape(4,1) # make column vectors class_sc_mat += (row-mv).dot((row-mv).T) S_W += class_sc_mat # sum class scatter matricesprint('within-class Scatter Matrix:n', S_W)
得到
within-class Scatter Matrix: [[ 38.9562 13.683 24.614 5.6556] [ 13.683 17.035 8.12 4.9132] [ 24.614 8.12 27.22 6.2536] [ 5.6556 4.9132 6.2536 6.1756]]2.1b 贝塞尔纠偏项的影响
在计算类-协方差矩阵时,可以在类内散布矩阵上添加尺度因数 1N−1,这样计算式就变为
Ni 是类样本大小(此处为50),因为三个类别的样本数是相同的,这里可以舍去 (Ni−1)。
因为本征向量是相同的,只是本征值有一个常数项的尺度变化,所以即便将其忽略不计,最后得到的特征空间也不会改变(这一点在文末还有体现)。
2.2 类间散布矩阵 SB
mm 是全局均值,而 mmi 和 Ni 是每类样本的均值和样本数。
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) # make column vector overall_mean = overall_mean.reshape(4,1) # make column vector S_B += n * (mean_vec - overall_mean).dot((mean_vec - overall_mean).T)print('between-class Scatter Matrix:n', S_B)
得到
between-class Scatter Matrix: [[ 63.2121 -19.534 165.1647 71.3631] [ -19.534 10.9776 -56.0552 -22.4924] [ 165.1647 -56.0552 436.6437 186.9081] [ 71.3631 -22.4924 186.9081 80.6041]]第三步:求解矩阵 S−1WSB 的广义本征值
接下来求解矩阵 S−1WSB 的广义本征值问题,从而得到线性判别“器”:
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('nEigenvector {}: n{}'.format(i+1, eigvec_sc.real)) print('Eigenvalue {:}: {:.2e}'.format(i+1, eig_vals[i].real))
得到
Eigenvector 1:[[-0.2049] [-0.3871] [ 0.5465] [ 0.7138]]Eigenvalue 1: 3.23e+01Eigenvector 2:[[-0.009 ] [-0.589 ] [ 0.2543] [-0.767 ]]Eigenvalue 2: 2.78e-01Eigenvector 3:[[ 0.179 ] [-0.3178] [-0.3658] [ 0.6011]]Eigenvalue 3: -4.02e-17Eigenvector 4:[[ 0.179 ] [-0.3178] [-0.3658] [ 0.6011]]Eigenvalue 4: -4.02e-17
计算结果应当作何解释呢?从第一节线性代数课开始我们就知道,本征向量和本征值表示了一个线性变换的形变程度——本征向量是形变的方向;本征值是本征向量的尺度因数,刻画了形变的幅度。
如果将LDA用于降维,本征向量非常重要,因为它们将会组成新特征子空间的坐标轴。对应的本征值表示了这些新坐标轴的信息量的多少。
再检查一遍计算过程,然后对本征值做进一步讨论。
检查本征向量-本征值的运算
要检查该运算是否正确,看其是否满足等式:
for i in range(len(eig_vals)): eigv = eig_vecs[:,i].reshape(4,1) np.testing.assert_array_almost_equal(np.linalg.inv(S_W).dot(S_B).dot(eigv), eig_vals[i] * eigv, decimal=6, err_msg='', verbose=True)print('ok')
输出:
ok第四步:选择线性判别“器”构成新的特征子空间
4.1 按本征值降序排列本征向量
本征向量只是定义了新坐标轴的方向,它们的大小都是单位长度1。
想要构成更低维的子空间,就得选择丢弃哪个(些)本征向量,所以我们得看看对应本征向量的那些本征值。 大体上说,对于数据的分布情况,本征值最小的那些本征向量所承载的信息最少,所以要舍弃的就是这些本征向量。
将本征向量根据本征值的大小从高到低排序,然后选择前 k 个本征向量:
eig_pairs = [(np.abs(eig_vals[i]), eig_vecs[:,i]) for i in range(len(eig_vals))]# Sort the (eigenvalue, eigenvector) tuples from high to loweig_pairs = sorted(eig_pairs, key=lambda k: k[0], reverse=True)# Visually confirm that the list is correctly sorted by decreasing eigenvaluesprint('Eigenvalues in decreasing order:n')for i in eig_pairs: print(i[0])
得到
Eigenvalues in decreasing order:32.27195779970.277566863845.71450476746e-155.71450476746e-15
可以看到有两个本征值接近0,不是因为它们不含信息量,而是浮点数精度的关系。其实,这后两个本征值应该恰好为0。
在LDA中,线性判别器的数目最多是 c−1,c 是总的类别数,这是因为类内散布矩阵 SB 是 c 个秩为1或0的矩阵的和。
注意到很少有完全共线的情况(所有样本点分布在一条直线上),协方差矩阵秩为1,这导致了只有一个非零本征值和一个对应的本征向量。
现在,把“可释方差”表达为百分数的形式:
print('Variance explained:n')eigv_sum = sum(eig_vals)for i,j in enumerate(eig_pairs): print('eigenvalue {0:}: {1:.2%}'.format(i+1, (j[0]/eigv_sum).real))
得到
Variance explained:eigenvalue 1: 99.15%eigenvalue 2: 0.85%eigenvalue 3: 0.00%eigenvalue 4: 0.00%
第一个本征对是信息量最大的一组,如果我们考虑建立一个一维的向量空间,使用该本征对就不会丢失太多信息。
4.2 选择 k 个最大本征值对应的本征向量
按本征值的大小得到降序排列的本征对之后,现在就可以组建我们的 k×d-维本征向量矩阵 W 了(此时大小为 4×2),这样就从最初的4维特征空间降到了2维的特征空间。
W = np.hstack((eig_pairs[0][1].reshape(4,1), eig_pairs[1][1].reshape(4,1)))print('Matrix W:n', W.real)
得到
Matrix W: [[-0.2049 -0.009 ] [-0.3871 -0.589 ] [ 0.5465 0.2543] [ 0.7138 -0.767 ]]
第五步:将样本变换到新的子空间中
使用上一步刚刚计算出的 4×2-维矩阵 W, 将样本变换到新的特征空间上:
X_lda = X.dot(W)assert X_lda.shape == (150,2), "The matrix is not 150x2 dimensional."
得到
from matplotlib import pyplot as pltdef 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: Iris projection onto the first 2 linear discriminants') # hide axis ticks plt.tick_params(axis="both", which="both", bottom="off", top="off", labelbottom="on", left="off", right="off", labelleft="on") # remove axis spines 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()
上方散点图展示了我们通过 LDA 构建的新的特征子空间。可以看到第一个线性判别器“LD1”把不同类数据区分得不错,第二个线性判别器就不行了。其原因在上面已经做了简单解释。
PCA 和 LDA 的对比
为了与使用线性判别分析得到的特征子空间作比较,我们将使用 scikit-learn 机器学习库中的 PCA 类。
文档看这里:
http://scikit-learn.org/stable/modules/generated/sklearn.decomposition.PCA.html。
方便起见,我们直接在 n_components 参数上设定了希望从输入数据集中得到的成分数量。
n_components : int, None or stringNumber of components to keep. if n_components is not set all components are kept: n_components == min(n_samples, n_features) if n_components == ‘mle’, Minka’s MLE is used to guess the dimension if 0 < n_components < 1, select the number of components such that the amount of variance that needs to be explained is greater than the percentage specified by n_components
但是在观察两种线性变换的结果之前,让我们快速复习一下 PCA 和 LDA 的目标:PCA 在整个数据集中寻找方差最大的坐标轴,而 LDA 则寻找对于类别区分度最佳的坐标轴。
实际的降维步骤中,经常是LDA放在PCA之后使用。
from sklearn.decomposition import PCA as sklearnPCAsklearn_pca = sklearnPCA(n_components=2)X_pca = sklearn_pca.fit_transform(X)def plot_pca(): ax = plt.subplot(111) for label,marker,color in zip( range(1,4),('^', 's', 'o'),('blue', 'red', 'green')): plt.scatter(x=X_pca[:,0][y == label], y=X_pca[:,1][y == label], marker=marker, color=color, alpha=0.5, label=label_dict[label] ) plt.xlabel('PC1') plt.ylabel('PC2') leg = plt.legend(loc='upper right', fancybox=True) leg.get_frame().set_alpha(0.5) plt.title('PCA: Iris projection onto the first 2 principal components') # hide axis ticks plt.tick_params(axis="both", which="both", bottom="off", top="off", labelbottom="on", left="off", right="off", labelleft="on") # remove axis spines 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.tight_layout plt.grid() plt.show()
然后执行
plot_pca()plot_step_lda()
回想我们之前讨论的内容:PCA 获取数据集中方差最大的成分,而LDA给出不同类别间方差最大的坐标轴。
使用 scikit-learn 中的 LDA
我们已经看到,线性判别分析是如何一步步实现的了。其实通过使用 scikit-learn
机器学习库中的 LDA ,我们可以更方便地实现同样的结果。
from sklearn.discriminant_analysis import LinearDiscriminantAnalysis as LDA# LDAsklearn_lda = LDA(n_components=2)X_lda_sklearn = sklearn_lda.fit_transform(X, y)
作图:
def plot_scikit_lda(X, title): ax = plt.subplot(111) for label,marker,color in zip( range(1,4),('^', 's', 'o'),('blue', 'red', 'green')): plt.scatter(x=X[:,0][y == label], y=X[:,1][y == label] * -1, # flip the figure 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(title) # hide axis ticks plt.tick_params(axis="both", which="both", bottom="off", top="off", labelbottom="on", left="off", right="off", labelleft="on") # remove axis spines 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()plot_scikit_lda(X_lda_sklearn, title='Default LDA via scikit-learn')
关于规范化
最近被问到了这个问题,所以我做一下说明:对特征做尺度变换,比如【规范化】,不会改变LDA的最终结果,所以这是可选步骤。的确,散布矩阵会因为特征的缩放而发生变化,而且本征向量也因此而改变。但是关键的部分在于,最终本征值不会变——唯一可见的不同只是成分轴的尺度大小。这在数学上是可以证明的(将来我可能在这儿插入一些公式)。
%matplotlib inlineimport pandas as pdimport matplotlib.pyplot as pltimport pandas as pddf = pd.read_csv('https://archive.ics.uci.edu/ml/machine-learning-databases/iris/iris.data', header=None)df[4] = df[4].map({'Iris-setosa':0, 'Iris-versicolor':1, 'Iris-virginica':2})df.tail()
sepal length in cm
sepal width in cm
petal length in cm
pedal width in cm
class label
145
6.7
3.0
5.2
2.3
Iris-virginica
146
6.3
2.5
5.0
1.9
Iris-virginaica
147
6.5
3.0
5.2
2.0
Iris-virginaica
148
6.2
3.4
5.4
2.3
Iris-virginaica
149
5.9
3.0
5.1
1.8
Iris-virginaica
读取数据集,对 X 中的列做规范化。规范化就是把数据用均值做中心化、用标准差做单位化:
这样所有的列就都是 0 均值(μxstd=0)、标准差为 1 的了(sigma_{x_{std}}=1)。
y, X = df.iloc[:, 4].values, df.iloc[:, 0:4].valuesX_cent = X - X.mean(axis=0)X_std = X_cent / X.std(axis=0)
接下来我简单复制了一下LDA的各个步骤,这些之前我们都讨论过了。为简便都写成了Python函数。
import numpy as npdef comp_mean_vectors(X, y): class_labels = np.unique(y) n_classes = class_labels.shape[0] mean_vectors = [] for cl in class_labels: mean_vectors.append(np.mean(X[y==cl], axis=0)) return mean_vectorsdef scatter_within(X, y): class_labels = np.unique(y) n_classes = class_labels.shape[0] n_features = X.shape[1] mean_vectors = comp_mean_vectors(X, y) S_W = np.zeros((n_features, n_features)) for cl, mv in zip(class_labels, mean_vectors): class_sc_mat = np.zeros((n_features, n_features)) for row in X[y == cl]: row, mv = row.reshape(n_features, 1), mv.reshape(n_features, 1) class_sc_mat += (row-mv).dot((row-mv).T) S_W += class_sc_mat return S_Wdef scatter_between(X, y): overall_mean = np.mean(X, axis=0) n_features = X.shape[1] mean_vectors = comp_mean_vectors(X, y) S_B = np.zeros((n_features, n_features)) for i, mean_vec in enumerate(mean_vectors): n = X[y==i+1,:].shape[0] mean_vec = mean_vec.reshape(n_features, 1) overall_mean = overall_mean.reshape(n_features, 1) S_B += n * (mean_vec - overall_mean).dot((mean_vec - overall_mean).T) return S_Bdef get_components(eig_vals, eig_vecs, n_comp=2): n_features = X.shape[1] 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) W = np.hstack([eig_pairs[i][1].reshape(4, 1) for i in range(0, n_comp)]) return W
先把还没有做尺度变换的数据的本征值、本征向量和投影矩阵打印出来:
S_W, S_B = scatter_within(X, y), scatter_between(X, y)eig_vals, eig_vecs = np.linalg.eig(np.linalg.inv(S_W).dot(S_B))W = get_components(eig_vals, eig_vecs, n_comp=2)print('EigVals: %snnEigVecs: %s' % (eig_vals, eig_vecs))print('nW: %s' % W)
结果:
EigVals: [ 2.0905e+01 +0.0000e+00j 1.4283e-01 +0.0000e+00j -2.8680e-16 +1.9364e-15j -2.8680e-16 -1.9364e-15j]EigVecs: [[ 0.2067+0.j 0.0018+0.j 0.4846-0.4436j 0.4846+0.4436j] [ 0.4159+0.j -0.5626+0.j 0.0599+0.1958j 0.0599-0.1958j] [-0.5616+0.j 0.2232+0.j 0.1194+0.1929j 0.1194-0.1929j] [-0.6848+0.j -0.7960+0.j -0.6892+0.j -0.6892-0.j ]]W: [[ 0.2067+0.j 0.0018+0.j] [ 0.4159+0.j -0.5626+0.j] [-0.5616+0.j 0.2232+0.j] [-0.6848+0.j -0.7960+0.j]]
作图:
X_lda = X.dot(W)for label,marker,color in zip( np.unique(y),('^', 's', 'o'),('blue', 'red', 'green')): plt.scatter(X_lda[y==label, 0], X_lda[y==label, 1], color=color, marker=marker)
再对规范化之后的数据集重复上述操作:
S_W, S_B = scatter_within(X_std, y), scatter_between(X_std, y)eig_vals, eig_vecs = np.linalg.eig(np.linalg.inv(S_W).dot(S_B))W_std = get_components(eig_vals, eig_vecs, n_comp=2)print('EigVals: %snnEigVecs: %s' % (eig_vals, eig_vecs))print('nW: %s' % W_std)
得到
EigVals: [ 2.0905e+01 1.4283e-01 -6.7207e-16 1.1082e-15]EigVecs: [[ 0.1492 -0.0019 0.8194 -0.3704] [ 0.1572 0.3193 -0.1382 -0.0884] [-0.8635 -0.5155 -0.5078 -0.5106] [-0.4554 0.7952 -0.2271 0.7709]]W: [[ 0.1492 -0.0019] [ 0.1572 0.3193] [-0.8635 -0.5155] [-0.4554 0.7952]]
作图:
X_std_lda = X_std.dot(W_std)X_std_lda[:, 0] = -X_std_lda[:, 0]X_std_lda[:, 1] = -X_std_lda[:, 1]for label,marker,color in zip( np.unique(y),('^', 's', 'o'),('blue', 'red', 'green')): plt.scatter(X_std_lda[y==label, 0], X_std_lda[y==label, 1], color=color, marker=marker)
可以看到,不管做不做尺度变化,本征值是一样的(要注意到,W 的秩是2,4维数据集中最小的两个本征值本来就应该是0)。并且,你也看到了,除了成分轴的尺度不太一样,以及图形做了中心对称翻转,最后的投影结果基本没有区别。
注1:文中出现了线性代数术语“eigenvalue”“eigenvector”,中文教材对应有“特征值”“本征值”两种常见译法。为了与“feature”相区分,本文使用“本征”翻译。
注2:文中提到 “k×d-维本征向量矩阵”,原文写作 “k×d-dimensional eigenvector matrix”,指的是“k个d维的本征向量组成的矩阵”。
中文与公式混合排版之后可能引发歧义,故在此做单独解释。
译者 / 李宇琛Asher Li
原文来源
http://sebastianraschka.com/Articles/2014_python_lda.html
谷歌开发大使Joshua Gordon:“所谓深度学习,就是让算法帮你解决问题”返回搜狐,查看更多