在本章中,我们将会学习到三种特征提取的方法,它们都可以将原始数据集变换到一个维度更低的新的特征子空间。

无监督数据降维技术之主成分分析

主成分分析(PCA)是一种广泛应用于不同领域的无监督线性数据转换技术,突出作用是降维。PCA的目标是在高维数据中找到最大方差的方向,并且将数据映射到一个维度不大于原始数据的新的子空间上。

如果使用PCA技术,我们需要构建一个$ d * k $维的转换矩阵$ W $,从而将原来的d维特征向量转换为k维特征向量(k<d)。PCA算法的步骤如下:

  1. 对原始d维数据做标准化处理
  2. 构造样本的协方差矩阵
  3. 计算协方差矩阵的特征值和相应的特征向量
  4. 选择前k个最大特征对应的特征向量(k为新的特征空间维度)
  5. 通过前k个特征向量构建映射矩阵$ W $
  6. 将原始的d维特征$ x $通过$ W $转换为新的k维特征$ x’ $

总体方差和贡献方差

这一小节完成PCA的前四个步骤。

首先,使用前面用到的葡萄酒数据集:

1
2
import pandas as pd
df_wine = pd.read_csv('https://archive.ics.uci.edu/ml/machine-learning-databases/wine/wine.data', header=None)

接着,将数据集划分为训练集和测试集,同时使用StandardScaler来将其标准化:

1
2
3
4
5
6
7
8
from sklearn.model_selection import train_test_split
from sklearn.preprocessing import StandardScaler

X, y = df_wine.iloc[:, 1:].values, df_wine.iloc[:, 0].values
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.3, random_state=0)
sc = StandardScaler()
X_train_std = sc.fit_transform(X_train)
X_test_std = sc.fit_transform(X_test)

接下来构造协方差矩阵,同时求解协方差矩阵的特征值和特征向量:

1
2
3
4
5
6
7
import numpy as np
cov_mat = np.cov(X_train_std.T)
eigen_vals, eigen_vecs = np.linalg.eig(cov_mat)
print(eigen_vals)
>> [4.8923083 2.46635032 1.42809973 1.01233462 0.84906459 0.60181514
>> 0.52251546 0.08414846 0.33051429 0.29595018 0.16831254 0.21432212
>> 0.2399553 ]

通过使用np.linalg.elg函数,可以得到一个包含有13个特征值的向量(eigen_vals)和一个13 * 13的特征矩阵(eigen_vecs),其中,特征向量以列的方式存在于特征矩阵中。

由于我们需要将数据压缩到一个新的特征子空间上实现降维,我们只需要选择那些包含最多信息的特征向量组成的子集。在此衡量函数是特征值$ \lambda_j $的方差贡献率:
$$
\frac{\lambda_j}{\sum_{i=1}^{d}j}
$$
接下来看一下不同特征值对应的方差贡献率:

1
2
3
4
tot = sum(eigen_vals)
var_exp = [(i / tot) for i in sorted(eigen_vals, reverse=True)]
print(var_exp)
>> [0.3732964772349068, 0.18818926106599568, 0.10896790724757796, 0.07724389477124863, 0.0647859460182618, 0.045920138114781475, 0.03986935597634714, 0.025219142607261574, 0.022581806817679666, 0.01830924471952691, 0.016353362655051454, 0.01284270583749274, 0.006420756933868311]

可以知道,第一主成分占方差总和的$ 40% $左右。

特征转换

接下来继续执行PCA方法的最后三个步骤。

首先,按照特征值的降序排列特征对:

1
2
eigen_pairs = [(np.abs(eigen_vals[i]), eigen_vecs[:, i]) for i in range(len(eigen_vals))]
eigen_pairs.sort(reverse=True)

接下来,我们只选择两个对应的最大的特征向量:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
w = np.hstack((eigen_pairs[0][1][:, np.newaxis],
eigen_pairs[1][1][:, np.newaxis]))
print(w)
>> [[ 0.14669811 0.50417079]
>> [-0.24224554 0.24216889]
>> [-0.02993442 0.28698484]
>> [-0.25519002 -0.06468718]
>> [ 0.12079772 0.22995385]
>> [ 0.38934455 0.09363991]
>> [ 0.42326486 0.01088622]
>> [-0.30634956 0.01870216]
>> [ 0.30572219 0.03040352]
>> [-0.09869191 0.54527081]
>> [ 0.30032535 -0.27924322]
>> [ 0.36821154 -0.174365 ]
>> [ 0.29259713 0.36315461]]

从而我们现在得到了一个13*2的映射矩阵$ W $。接下来转换原始的数据集:

1
X_train_pca = X_train_std.dot(w)

最后,新的数据集被保存在124*2的矩阵中,接下来对其进行可视化:

1
2
3
4
5
6
7
8
9
import matplotlib.pyplot as plt
colors = ['r', 'b', 'g']
markers = ['s', 'x', 'o']
for l, c, m in zip(np.unique(y_train), colors, markers):
plt.scatter(X_train_pca[y_train==l, 0], X_train_pca[y_train==l, 1], c=c, label=l, marker=m)
plt.xlabel('PC 1')
plt.ylabel('PC 2')
plt.legend()
plt.show()

得到的图像如下:

img

从上图可以很直观的看到,线性分类器能够对其有很好的划分。

使用scikit-learn进行主成分分析

我们先使用PCA对葡萄酒数据做预处理,然后再使用逻辑斯蒂回归模型对转换后的数据进行分类,最后绘制出散点图:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
from sklearn.decomposition import PCA

pca = PCA(n_components=2)
X_train_pca = pca.fit_transform(X_train_std)
X_test_pca = pca.transform(X_test_std)

import matplotlib.pyplot as plt
colors = ['r', 'b', 'g']
markers = ['s', 'x', 'o']
for l, c, m in zip(np.unique(y_train), colors, markers):
plt.scatter(X_train_pca[y_train==l, 0], X_train_pca[y_train==l, 1], c=c, label=l, marker=m)
plt.xlabel('PC 1')
plt.ylabel('PC 2')
plt.legend()
plt.show()

得到的图像如下:

img

比较该图和上一节中的图像,可以发现上图实际上就是我们自己完成的PCA图沿着PC1轴翻转的结果。出现此差异的原因在于特征分析方法:特征向量为正或者为负。

接下来使用逻辑斯蒂回归模型进行训练,并且得到训练结果:

1
2
3
4
5
6
7
8
from sklearn.linear_model import LogisticRegression
from sklearn.metrics import accuracy_score

lr = LogisticRegression()
lr.fit(X_train_pca, y_train)
y_pred = lr.predict(X_test_pca)
accuracy_score(y_test, y_pred)
>> 0.9814814814814815

可以发现的逻辑斯蒂回归模型的拟合率很优良。

通过线性判别分析压缩无监督数据

线性判别分析(LDA)是一种可作为特征抽取的技术,它可以提高数据分析过程中的计算效率,同时,对于不适用于正则化的模型,它可以降低因维度灾难带来的过拟合。

LDA方法的步骤如下:

  1. 对d为数据集进行标准化处理
  2. 对于每一类别,计算d维的均值向量
  3. 构造类间的散布矩阵$ S_{B} $以及类内的散布举证$ S_{W} $
  4. 计算矩阵$ s_{W}^{-1}S_{B} $的特征值及对应的特征向量
  5. 选取前k个特征值对应的特征向量,构造一个d*k维的转换矩阵$ W $
  6. 使用转换矩阵$ W $将样本映射到新的特征子空间中

计算散布矩阵

葡萄酒数据我们已经经过标准化处理,接下来求解均值向量$ m_i $:

1
2
3
4
5
6
7
8
9
10
11
np.set_printoptions(precision=4)
mean_vecs = []
for label in range(1, 4):
mean_vecs.append(np.mean(X_train_std[y_train==label], axis=0))
print(mean_vecs)
>> [array([ 0.9259, -0.3091, 0.2592, -0.7989, 0.3039, 0.9608, 1.0515,
>> -0.6306, 0.5354, 0.2209, 0.4855, 0.798 , 1.2017]),
>> array([-0.8727, -0.3854, -0.4437, 0.2481, -0.2409, -0.1059, 0.0187,
>> -0.0164, 0.1095, -0.8796, 0.4392, 0.2776, -0.7016]),
>> array([ 0.1637, 0.8929, 0.3249, 0.5658, -0.01 , -0.9499, -1.228 ,
>> 0.7436, -0.7652, 0.979 , -1.1698, -1.3007, -0.3912])]

通过均值向量,我们计算一下类内散布矩阵$ S_W $:
$$
S_W = \sum_{i=1}^cS_i
$$
这可以通过累加各类别i的散步矩阵$ S_i $来计算:
$$
S_i = \sum_{x \in D_i}^{c}(x-m_i)(x-m_i)^T
$$

1
2
3
4
5
6
7
8
9
10
d = 13
S_W = np.zeros((d, d))
for label, mv in zip(range(1, 4), mean_vecs):
class_scater = np.zeros((d, d))
for row in X[y==label]:
row, mv = row.reshape(d, 1), mv.reshape(d, 1)
class_scater += (row - mv).dot((row - mv).T)
S_W += class_scater
S_W.shape
>> (13, 13)

此前,我们假定对散步矩阵计算时,曾假设训练集的类标是均匀分布的,但是,通过以下程序,我们发现其不遵守这个假设:

1
2
np.bincount(y_train)
>> array([ 0, 40, 49, 35], dtype=int64)

因此,在我们通过累加方式计算散布矩阵$ S_{W} $前,需要对各类别的散步矩阵$ S_i $做缩放处理。但采用此种方式时,此时散布矩阵和协方差矩阵计算方式相同。协方差矩阵可以看作是归一化的散布矩阵:
$$
\frac{1}{N_i}S_{W} = \frac{1}{N_i}\sum_{x \in D_i}^c(x-m_i)(x-m_i)^T
$$

1
2
3
4
5
6
7
d = 13
S_W = np.zeros((d, d))
for label, mv in zip(range(1, 4), mean_vecs):
class_scater = np.cov(X_train_std[y_train==label].T)
S_W += class_scater
S_W.shape
>> (13, 13)

接下来计算类间散布矩阵$ S_B $:
$$
S_B = \sum_{i=1}^cN_i(m_i-m)(m_i-m)^T
$$
其中,m是全局均值,他在计算时用到了所有类别中的全部样本:

1
2
3
4
5
6
7
8
9
10
mean_overall = np.mean(X_train_std, axis=0)
d = 13
S_B = np.zeros((d, d))
for i, mean_vec in enumerate(mean_vecs):
n = X[y==i+1, :].shape[0]
mean_vec = mean_vec.reshape(d, 1)
mean_overall = mean_overall.reshape(d, 1)
S_B += n * (mean_vec - mean_overall).dot((mean_vec - mean_overall).T)
S_B.shape
>> (13, 13)

在新特征子空间上选取线性判别算法

LDA余下的步骤和PCA的步骤相似:

1
2
3
4
5
eigen_vals, eigen_vecs = np.linalg.eig(np.linalg.inv(S_W).dot(S_B))
eigen_pairs = [(np.abs(eigen_vals[i]), eigen_vecs[:, i]) for i in range(len(eigen_vals))]
eigen_pairs = sorted(eigen_pairs, key=lambda k: k[0], reverse=True)
for eigen_pair in eigen_pairs:
print(eigen_pair[0])

得到的结果如下:

1
2
3
4
5
6
7
8
9
10
11
12
13
643.0153843460517
225.08698185416256
8.002675183788468e-14
5.757534614184537e-14
3.5105079604736804e-14
3.4638958368304884e-14
2.587811510007498e-14
2.587811510007498e-14
2.4449817310582036e-14
1.6532199129716054e-14
8.331225171347768e-15
2.3238388797036527e-15
6.522430076120113e-16

从上述输出来看,我们只得到了两个非零特征值(实际得到的3-13个特征值并未严格为0,这是由numpy的浮点数运算导致的),说明只有前面两个特征值对应的特征几乎包含了葡萄酒训练数据集中的全部有用信息。

接下来构造转换矩阵:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
w = np.hstack((eigen_pairs[0][1][:, np.newaxis].real, eigen_pairs[1][1][:, np.newaxis].real))
w
>> array([[-0.0707, 0.3778],
>> [ 0.0359, 0.2223],
>> [-0.0263, 0.3813],
>> [ 0.1875, -0.2955],
>> [-0.0033, -0.0143],
>> [ 0.2328, -0.0151],
>> [-0.7719, -0.2149],
>> [-0.0803, -0.0726],
>> [ 0.0896, -0.1767],
>> [ 0.1815, 0.2909],
>> [-0.0631, -0.2376],
>> [-0.3794, -0.0867],
>> [-0.3355, 0.586 ]])

将样本映射到新的特征空间

通过上一节中构建的转换矩阵$ W $,我们来对原始数据进行转换:

1
2
3
4
5
6
7
8
9
X_train_lda = X_train_std.dot(w)
colors = ['r', 'b', 'g']
markers = ['s', 'x', 'o']
for l, c, m in zip(np.unique(y_train), colors, markers):
plt.scatter(X_train_lda[y_train==l, 0], X_train_lda[y_train==l, 1], c=c, label=l, marker=m)
plt.xlabel('LD 1')
plt.ylabel('LD 2')
plt.legend()
plt.show()

得到的图像如下:

img

通过图像可知,三个葡萄酒类在新的特征子空间上是线性可分的。

使用scikit-laern进行LDA分析

接下来,看一下scikit-laern中对LDA类的实现:

1
2
3
4
5
6
7
8
9
10
11
from sklearn.discriminant_analysis import LinearDiscriminantAnalysis as LDA
lda = LDA(n_components=2)
X_train_lda = lda.fit_transform(X_train_std, y_train)
colors = ['r', 'b', 'g']
markers = ['s', 'x', 'o']
for l, c, m in zip(np.unique(y_train), colors, markers):
plt.scatter(X_train_lda[y_train==l, 0], X_train_lda[y_train==l, 1], c=c, label=l, marker=m)
plt.xlabel('LD 1')
plt.ylabel('LD 2')
plt.legend()
plt.show()

得到的图像如下:

img

此时看一下逻辑斯蒂回归模型的预测准确度:

1
2
3
4
5
6
lr = LogisticRegression()
lr = lr.fit(X_train_lda, y_train)
X_test_lda = lda.fit_transform(X_test_std, y_test)
y_pred = lr.predict(X_test_lda)
accuracy_score(y_pred, y_test)
>> 1.0

可以看到,逻辑斯蒂回归模型在测试数据集上对样本分类可谓完美。

使用核主成分分析进行非线性映射

许多机器学习算法都假定输入数据是线性可分的,但是在现实世界中,大多数的数据是线性不可分的,针对此类问题,使用PCA或者LDA等降维技术,将其转化为线性问题并不是最好的方法。在本节中,我们将了解一下利用核技巧的PCA,或者称其为核PCA,这和第三章中我们介绍的核支持向量机的概念有一定的联系。使用核PCA,我们将学习如何将非线性可分的数据转换到一个适合对其进行线性分类的新的低维子空间中。

核函数

通过核PCA,我们能够得到已经映射到各成分的样本,而不像标准PCA那样去构建一个转换矩阵。简单地说,可以将核函数理解为:通过两个向量点积来度量向量间相似度的函数。最常用的核函数有:

  • 多项式核:
    $$
    k(x^i, x^j) = (x^{iT}x^j + \theta)^p
    $$
    其中,阈值$ \theta $和幂的值$ p $需要自行定义。

  • 双曲正切(sigmoid)核:
    $$
    k(x^i, x^j) = thah(\eta x^{iT}x^j+\theta)
    $$

  • 径向基核函数(RBF)或者称为高斯核函数:
    $$
    k(x^i, x^j) = exp\left(-\frac{||x^i-x^j||^2}{2\sigma^2}\right)
    = exp(-\gamma||x^i - x^j||^2)
    $$

基于RBF核的PCA可以通过如下三个步骤实现:

  1. 为了计算核矩阵$ k $,我们需要做如下计算:
    $$
    k(x^i,x^j) = = exp(-\gamma||x^i - x^j||^2)
    $$
    我们需要计算任意两个样本对之间的值:
    $$
    K = \begin{bmatrix}
    k(x^1,x^1) & k(x^1, x^2) & \cdots & k(x^1, x^n)\
    k(x^2,x^1) & k(x^2, x^2) & \cdots & k(x^2, x^n)\
    \vdots & \vdots & \ddots & \vdots\
    k(x^n,x^1) & k(x^n, x^2) & \cdots & k(x^n, x^n)\
    \end{bmatrix}
    $$

  2. 通过如下公式,使得核矩阵$ k $更为聚集:
    $$
    K’ = K-l_nK-Kl_n+l_nKl_n
    $$
    其中, $ l_n $是一个n*n的矩阵,其所有的值都是$ \frac{1}{n} $。

  3. 将聚集后的核矩阵的特征值按照降序排列,选择前k个特征值对应的特征向量。和标准PCA不同,这里的特征向量不是主成分轴,而是将样本映射到这些轴上。

使用Python实现主成分分析

接下来,借助SciPy和NumPy的函数,我们手动实现一个核PCA:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
from scipy.spatial.distance import pdist, squareform
from scipy import exp
from scipy.linalg import eigh
import numpy as np

def rbf_kernel_pca(X, gamma, n_components):
sq_dists = pdist(X, 'sqeuclidean')
mat_sq_dists = squareform(sq_dists)
K = exp(-gamma * mat_sq_dists)
N = K.shape[0]
one_n = np.ones((N, N)) / N
K = K - one_n.dot(K) - k.dot(one_n) + one_n.dot(K).dot(one_n)
eigvals, eigvecs = eigh(K)
X_pc = np.column_stack((eigvecs[:, -i] for i in range(1, n_components + 1)))
return X_pc

接下来查看几个实例。

  1. 实例一:分离半月形数据

    首先创建一个包含100个样本点的二维数据集,以两个半月形状表示:

    1
    2
    3
    4
    5
    from sklearn.datasets import make_moons
    X, y = make_moons(n_samples=100, random_state=123)
    plt.scatter(X[y==0, 0], X[y==0, 1], color='r', marker='^', alpha=0.5)
    plt.scatter(X[y==1, 0], X[y==1, 1], color='b', marker='o', alpha=0.5)
    plt.show()

    得到的图像如下:

    img

    显然,这两个半月形不是线性可分的,我们的目标是通过核PCA将这两个半月形数据展开,使得数据集成为适用于某一线性分类器的输入数据。

    首先,我们看一下经过标准PCA处理的数据集的图像:

    1
    2
    3
    4
    5
    6
    from sklearn.decomposition import PCA
    scikit_pca = PCA(n_components=2)
    X_spca = scikit_pca.fit_transform(X)
    plt.scatter(X_spca[y==0, 0], X_spca[y==0, 1], color='r', marker='^', alpha=0.5)
    plt.scatter(X_spca[y==1, 0], X_spca[y==1, 1], color='b', marker='o', alpha=0.5)
    plt.show()

    得到的图像如下:

    img

    可以发现,经过标准化PCA处理后,线性分类器未必能很好地发挥作用。

    接下来尝试一下核PCA函数:

    1
    2
    3
    4
    X_kpca = rbf_kernel_pca(X, gamma=15, n_components=2)
    plt.scatter(X_kpca[y==0, 0], X_kpca[y==0, 1], color='r', marker='^', alpha=0.5)
    plt.scatter(X_kpca[y==1, 0], X_kpca[y==1, 1], color='b', marker='o', alpha=0.5)
    plt.show()

    得到的图像如下:

    img

    可以看到,此时两个类别是线性可分的。

  2. 示例二:分离同心圆

    接下俩看一下非线性相关的另外一个例子:同心圆:

    1
    2
    3
    4
    5
    from sklearn.datasets import make_circles
    X, y = make_circles(n_samples=1000, random_state=123, noise=0.1, factor=0.2)
    plt.scatter(X[y==0, 0], X[y==0, 1], color='r', marker='^', alpha=0.5)
    plt.scatter(X[y==1, 0], X[y==1, 1], color='b', marker='o', alpha=0.5)
    plt.show()

    得到的图像如下:

    img

    接下来使用核PCA,观察数据集分布:

    1
    2
    3
    4
    X_kpca = rbf_kernel_pca(X, gamma=15, n_components=2)
    plt.scatter(X_kpca[y==0, 0], X_kpca[y==0, 1], color='r', marker='^', alpha=0.5)
    plt.scatter(X_kpca[y==1, 0], X_kpca[y==1, 1], color='b', marker='o', alpha=0.5)
    plt.show()

    得到的图像如下:

    img

    可以发现,此时两个类别的数据是线性可分的。

映射新的数据点

在标准PCA方法中,我们通过转换矩阵和输入样本之间的点积来对数据进行映射。但是在核PCA中,该如何转换型的数据点呢?实际上,如果我们希望将新的样本$ x’ $映射到此主成分轴,需要进行如下计算:
$$
\phi(x’)^Tv
$$
幸运的是,我们可以使用核技巧,这样就无需精确计算映射$ \phi(x’)^Tv $。通过以下公式计算:
$$
\phi(x’)^Tv = \sum_ia^ik(x’, x^i)^T
$$
其中,核矩阵K的特征向量$ a $和特征值$ \lambda $关系如下:
$$
Ka = \lambda a
$$
通过如下程序实现映射:

1
2
3
4
def project_x(x_new, X, gamma, alphas, lambdas):
pair_dist = np.array([np.sum(x_new - row)**2 for row in X])
k = np.exp(-gamma * pair_dist)
return k.dot(alphas / lambdas)

其中,alphas是前k个特征向量,lambdas是前k个对应的特征值:

1
2
alphas = np.column_stack((eigvecs[:, i] for i in range(1, n_components+1)))
lambdas = [eigvals[-i] for i in range(1, n_components+1)]

将上述两条语句加到rbf_kernel_pca函数末端并且返回他们的值即可。

scikit-learn中的核主成分分析

使用scikit-learn中的API实现核PCA如下:

1
2
3
4
5
6
7
from sklearn.decomposition import KernelPCA
X, y = make_moons(n_samples=100, random_state=123)
scikit_kpca = KernelPCA(n_components=2, kernel='rbf', gamma=15)
X_skpca = scikit_kpca.fit_transform(X)
plt.scatter(X_skpca[y==0, 0], X_skpca[y==0, 1], color='r', marker='^', alpha=0.5)
plt.scatter(X_skpca[y==1, 0], X_skpca[y==1, 1], color='b', marker='o', alpha=0.5)
plt.show()

得到的图像如下:

img

从上图来看,scikit-learn中KernelPCA得到的结果核我们手动实现的结果相一致。