使用Python一步步实现PCA算法

使用Python一步步实现PCA算法

标签： PCA Python

---

本文原地址为：
http://sebastianraschka.com/Articles/2014_pca_step_by_step.html
Implementing a Principal Component Analysis (PCA)
– in Python, step by step
Apr 13, 2014
by Sebastian Raschka

此篇为翻译作品，仅作为学习用途。

简介

主成分分析（PCA）的主要目的是通过分析发现数据的模式进行维度的缩减，这个过程的原则是信息损失最小化。

我们希望得到的结果，把初始特征空间映射到一个相对低维度的子空间，同时保证这个低维度空间也能够很好的表达数据的有效信息。在模式分类中，我们希望通过降维操作抽取能够最佳表达数据的特征子集来降低运算时间花费，减少参数估计的误差。

主成分分析（PCA） vs 多元判别式分析（MDA）

PCA和MDA都是线性变换的方法，二者关系密切。在PCA中，我们寻找数据集中最大化方差的成分，在MDA中，我们对类间最大散布的方向更感兴趣。

一句话，通过PCA，我们将整个数据集（不带类别标签）映射到一个子空间中，在MDA中，我们致力于找到一个能够最好区分各类的最佳子集。粗略来讲，PCA是通过寻找方差最大的轴（在一类中，因为PCA把整个数据集当做一类），在MDA中，我们还需要最大化类间散布。

在通常的模式识别问题中，MDA往往在PCA后面。

什么是好的特征子集？

假设我们要把d维空间降到k维（k小于d），那么我们怎么选择k？还有就是我们怎么知道这个特征子集能够很好地代表整个数据集？
接下来，我们会计算数据集的特征向量，将他们保存在散布矩阵（或者通过协方差矩阵计算）中。每个特征向量对应一个特征值，这个值会告诉我们特征向量的长度或者量级。
如果我们观察到所有的特征值尺度相似，这表明我们的数据已经在一个很好的特征空间了，如果有的特征值远远大于其他的，我们会选择这些特征值，因为他们保存了更多的有效信息，相反，如果特征值趋近于0，表名这个维度包含的信息量极少，我们可以考虑扔掉这些维度。

总结一下PCA方法

下面列出实现PCA方法的6个步骤，我们会在随后的章节中一一详解。

1. 取 d $d$-维数据集，不包含类别标签
2. 计算d$d$维矩阵的均值向量（数据集每一维度的均值）
3. 计算数据集的散布矩阵（协方差矩阵）
4. 计算特征向量（ ee1,ee2,...,eed $\pmb e_1, \; \pmb e_2, \; ..., \; \pmb e_d$）以及对应的特征值( λλ1,λλ2,...,λλd $\pmb \lambda_1, \; \pmb \lambda_2, \; ..., \; \pmb \lambda_d$)
5. 对特征向量按照特征值降序排序，选择特征值最大的 k $k$个特征向量，形成一个d×k$d \times k$维的矩阵 WW $\pmb W\;$（每一列代表一个特征向量）
6. 使用这个 d×k $d \times k$维的特征向量矩阵将数据集映射到一个新的子空间上。可以用下面的数学公式描述： yy=WWT×xx $\pmb y = \pmb W^T \times \pmb x$（其中 xx $\pmb x$是一个 d×1 $d \times 1$维的向量代表一个样本，而 yy $\pmb y$是一个转换过的 k×1 $k \times 1$维的数据集）

生成一些3维样本数据

在下面的例子中，我们会随机生成40个符合多元高斯分布的3维样本。
在此，我们假设样本属于两个类，其中一半样本标记为 ω1 $\omega_1$（类别1），另外一个标记为 ω2 $\omega_2$（类别2）。
μ1μ1=⎡⎣⎢000⎤⎦⎥μ2μ2=⎡⎣⎢111⎤⎦⎥ $\pmb{\mu_1} = \begin{bmatrix}0\\0\\0\end{bmatrix} \pmb{\mu_2} = \begin{bmatrix}1\\1\\1\end{bmatrix}$ (样本均值)

Σ1Σ1=⎡⎣⎢100010001⎤⎦⎥Σ2Σ2=⎡⎣⎢100010001⎤⎦⎥ $\pmb{\Sigma_1} = \begin{bmatrix}1\quad 0\quad 0\\0\quad 1\quad0\\0\quad0\quad1\end{bmatrix} \pmb{\Sigma_2} = \begin{bmatrix}1\quad 0\quad 0\\0\quad 1\quad0\\0\quad0\quad1\end{bmatrix}$(协方差矩阵)

为什么我们选取3维样本

多维数据最大的问题在于可视化，比如PCA分析的可视化，我们本可以选取2维样本数据集，但是PCA目标在于降维，至少要将一个维度，因此从3维将到2维比较合适。

import numpy as np

np.random.seed(234234782384239784) # random seed for consistency

# A reader pointed out that Python 2.7 would raise a
# "ValueError: object of too small depth for desired array".
# This can be avoided by choosing a smaller random seed, e.g. 1
# or by completely omitting this line, since I just used the random seed for
# consistency.

mu_vec1 = np.array([0,0,0])
cov_mat1 = np.array([[1,0,0],[0,1,0],[0,0,1]])
class1_sample = np.random.multivariate_normal(mu_vec1, cov_mat1, 20).T
assert class1_sample.shape == (3,20), "The matrix has not the dimensions 3x20"

mu_vec2 = np.array([1,1,1])
cov_mat2 = np.array([[1,0,0],[0,1,0],[0,0,1]])
class2_sample = np.random.multivariate_normal(mu_vec2, cov_mat2, 20).T
assert class2_sample.shape == (3,20), "The matrix has not the dimensions 3x20"

通过上面的代码，我们生成了2个 3×20 $3\times20$的数据集，类别标签分布是 ω1 $\omega_1$， ω2 $\omega_2$。
其中每一列都可以用一个3维向量表示： xx=⎛⎝⎜x1x2x3⎞⎠⎟ $\pmb x = \begin{pmatrix} x_1 \\ x_2 \\ x_3 \end{pmatrix}$
我们的数据集就有如下的形式：

XX=⎛⎝⎜x11x12...x120x21x22...x220x31x32...x320⎞⎠⎟ $\pmb X = \begin{pmatrix} x_{1_1}\; x_{1_2} \; ... \; x_{1_{20}}\\ x_{2_1} \; x_{2_2} \; ... \; x_{2_{20}}\\ x_{3_1} \; x_{3_2} \; ... \; x_{3_{20}}\end{pmatrix}$

为了简要表示我们两个样本的分布情况，我们用3D散布图画一张图。

%pylab inline
from matplotlib import pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from mpl_toolkits.mplot3d import proj3d

fig = plt.figure(figsize=(8,8))
ax = fig.add_subplot(111, projection='3d')
plt.rcParams['legend.fontsize'] = 10   
ax.plot(class1_sample[0,:], class1_sample[1,:], class1_sample[2,:], 'o', markersize=8, color='blue', alpha=0.5, label='class1')
ax.plot(class2_sample[0,:], class2_sample[1,:], class2_sample[2,:], '^', markersize=8, alpha=0.5, color='red', label='class2')

plt.title('Samples for class 1 and class 2')
ax.legend(loc='upper right')

plt.show()

操作步骤

1. 取除了标签外的所有维度

因为PCA分析不需要类别标签，所以我们把两个类别的样本合并起来形成一个 3×40 $3\times40$的一个矩阵

all_samples = np.concatenate((class1_sample, class2_sample), axis=1)
assert all_samples.shape == (3,40), "The matrix has not the dimensions 3x40"

2.计算一个 dd $\pmb d$ 维均值向量

mean_x = np.mean(all_samples[0,:])
mean_y = np.mean(all_samples[1,:])
mean_z = np.mean(all_samples[2,:])

mean_vector = np.array([[mean_x],[mean_y],[mean_z]])

print('Mean Vector:\n', mean_vector)

  Mean Vector:
     [[ 0.50576644]
     [ 0.30186591]
     [ 0.76459177]]

3(a).计算散布矩阵

通过下面的公式计算散布矩阵
S=∑k=1n(xxk−mm)(xxk−mm)T $S = \sum\limits_{k=1}^n (\pmb x_k - \pmb m)\;(\pmb x_k - \pmb m)^T$

其中 mm $\pmb m$是均值向量，
mm=1n∑k=1nxxk $\pmb m = \frac{1}{n} \sum\limits_{k=1}^n \; \pmb x_k$

scatter_matrix = np.zeros((3,3))
for i in range(all_samples.shape[1]):
    scatter_matrix += (all_samples[:,i].reshape(3,1) - mean_vector).dot((all_samples[:,i].reshape(3,1) - mean_vector).T)
print('Scatter Matrix:\n', scatter_matrix)

 Scatter Matrix:
     [[ 48.91593255   7.11744916   7.20810281]
     [  7.11744916  37.92902984   2.7370493 ]
     [  7.20810281   2.7370493   35.6363759 ]]

3(b).计算协方差矩阵（散布矩阵的替代选择）

我们可以也可以通过使用内置的numpy.cov()函数来计算协方差矩阵来替代散布矩阵。计算协方差矩阵和计算散布矩阵的公式非常相似，唯一的区别在于，我们使用 1N−1 $\frac{1}{N-1}$作为尺度因子（此处是： 140−1=139 $\frac{1}{40-1} = \frac{1}{39}$）。因此，二者的特征空间是一样的（一样的特征向量，只是特征值计算的尺度因子常量不同）。
Σi=[σ211σ221σ231σ212σ222σ232σ213σ223σ233] $\Sigma_i = \Bigg[ \begin{array}{cc} \sigma_{11}^2 & \sigma_{12}^2 & \sigma_{13}^2\\ \sigma_{21}^2 & \sigma_{22}^2 & \sigma_{23}^2\\ \sigma_{31}^2 & \sigma_{32}^2 & \sigma_{33}^2\\ \end{array} \Bigg] %]]>$

cov_mat = np.cov([all_samples[0,:],all_samples[1,:],all_samples[2,:]])
print('Covariance Matrix:\n', cov_mat)

Covariance Matrix:
     [[ 1.25425468  0.1824987   0.18482315]
     [ 0.1824987   0.97253923  0.07018075]
     [ 0.18482315  0.07018075  0.91375323]]

4.计算特征向量以及对应的特征值

为了证明通过散布矩阵和协方差矩阵计算出来的特征向量是一样的，我们在代码中加一个断言assert，与此同时，我们也会看到从散布矩阵得到的特征值是被常量因子39整除的。

# eigenvectors and eigenvalues for the from the scatter matrix
eig_val_sc, eig_vec_sc = np.linalg.eig(scatter_matrix)

# eigenvectors and eigenvalues for the from the covariance matrix
eig_val_cov, eig_vec_cov = np.linalg.eig(cov_mat)

for i in range(len(eig_val_sc)):
    eigvec_sc = eig_vec_sc[:,i].reshape(1,3).T
    eigvec_cov = eig_vec_cov[:,i].reshape(1,3).T
    assert eigvec_sc.all() == eigvec_cov.all(), 'Eigenvectors are not identical'

    print('Eigenvector {}: \n{}'.format(i+1, eigvec_sc))
    print('Eigenvalue {} from scatter matrix: {}'.format(i+1, eig_val_sc[i]))
    print('Eigenvalue {} from covariance matrix: {}'.format(i+1, eig_val_cov[i]))
    print('Scaling factor: ', eig_val_sc[i]/eig_val_cov[i])
    print(40 * '-')

Eigenvector 1:
    [[-0.84190486]
     [-0.39978877]
     [-0.36244329]]
    Eigenvalue 1 from scatter matrix: 55.398855957302445
    Eigenvalue 1 from covariance matrix: 1.4204834860846791
    Scaling factor:  39.0
    ----------------------------------------
    Eigenvector 2:
    [[-0.44565232]
     [ 0.13637858]
     [ 0.88475697]]
    Eigenvalue 2 from scatter matrix: 32.42754801292286
    Eigenvalue 2 from covariance matrix: 0.8314755900749456
    Scaling factor:  39.0
    ----------------------------------------
    Eigenvector 3:
    [[ 0.30428639]
     [-0.90640489]
     [ 0.29298458]]
    Eigenvalue 3 from scatter matrix: 34.65493432806495
    Eigenvalue 3 from covariance matrix: 0.8885880596939733
    Scaling factor:  39.0
    ----------------------------------------

检查特征向量与特征值的计算

我们快速检查一下特征向量和特征值的计算是否正确，是否满足下面的公式
ΣΣvv=λvv $\pmb\Sigma\pmb{v} = \lambda\pmb{v}$

其中，
ΣΣ=Covariancematrix  $\pmb\Sigma = Covariance \; matrix\ $
vv=Eigenvector  $\pmb{v} = \; Eigenvector\ $
λ=Eigenvalue $\lambda = \; Eigenvalue$

for i in range(len(eig_val_sc)):
    eigv = eig_vec_sc[:,i].reshape(1,3).T
    np.testing.assert_array_almost_equal(scatter_matrix.dot(eigv), eig_val_sc[i] * eigv,
                                         decimal=6, err_msg='', verbose=True)

将特征向量可视化

在我们进行下一步操作之前，我们以样本均值为中心绘制特征向量。

%pylab inline

from matplotlib import pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from mpl_toolkits.mplot3d import proj3d
from matplotlib.patches import FancyArrowPatch


class Arrow3D(FancyArrowPatch):
    def __init__(self, xs, ys, zs, *args, **kwargs):
        FancyArrowPatch.__init__(self, (0,0), (0,0), *args, **kwargs)
        self._verts3d = xs, ys, zs

    def draw(self, renderer):
        xs3d, ys3d, zs3d = self._verts3d
        xs, ys, zs = proj3d.proj_transform(xs3d, ys3d, zs3d, renderer.M)
        self.set_positions((xs[0],ys[0]),(xs[1],ys[1]))
        FancyArrowPatch.draw(self, renderer)

fig = plt.figure(figsize=(7,7))
ax = fig.add_subplot(111, projection='3d')

ax.plot(all_samples[0,:], all_samples[1,:], all_samples[2,:], 'o', markersize=8, color='green', alpha=0.2)
ax.plot([mean_x], [mean_y], [mean_z], 'o', markersize=10, color='red', alpha=0.5)
for v in eig_vec_sc.T:
    a = Arrow3D([mean_x, v[0]], [mean_y, v[1]], [mean_z, v[2]], mutation_scale=20, lw=3, arrowstyle="-|>", color="r")
    ax.add_artist(a)
ax.set_xlabel('x_values')
ax.set_ylabel('y_values')
ax.set_zlabel('z_values')

plt.title('Eigenvectors')

plt.show()

5.1.对特征向量按照特征值降序排序

现在我们开始将特征空间进行降维，将样本的特征空间通过PCA映射到一个维度相对低一些的子空间，在这个子空间中，特征向量将是这个新的特征空间的轴，然而，特征向量仅仅定义了新轴的方向，因为他们的单位长度都是1，我们可以通过下面的代码确认。

for ev in eig_vec_sc:
    numpy.testing.assert_array_almost_equal(1.0, np.linalg.norm(ev))
    # instead of 'assert' because of rounding errors

所以，为了决定我们要剔除哪些特征向量，我们来看一下特征向量对应的特征值。粗略地说，特征值越小包含的数据分布的信息越小，这些就是我们需要剔除的对象。

最普遍的方法就是对特征值从高到底进行排序，选择前面 kk $\pmb k$个特征值。

# Make a list of (eigenvalue, eigenvector) tuples
eig_pairs = [(np.abs(eig_val_sc[i]), eig_vec_sc[:,i]) for i in range(len(eig_val_sc))]

# Sort the (eigenvalue, eigenvector) tuples from high to low
eig_pairs.sort(key=lambda x: x[0], reverse=True)

# Visually confirm that the list is correctly sorted by decreasing eigenvalues
for i in eig_pairs:
    print(i[0])

   55.3988559573
    34.6549343281
    32.4275480129

5.2.选择 kk $\pmb k$个特征值最大的特征向量

在本文的小例子中，我们将一个3维矩阵降维到2维矩阵，我们通过特征值的计算构建我们的 dd×k $\pmb d\times k$维特征向量 WW $\pmb W$

matrix_w = np.hstack((eig_pairs[0][1].reshape(3,1), eig_pairs[1][1].reshape(3,1)))
print('Matrix W:\n', matrix_w)

    Matrix W:
     [[-0.84190486  0.30428639]
     [-0.39978877 -0.90640489]
     [-0.36244329  0.29298458]]

6.将样本映射到新的子空间

在最后一步中，我们使用一个刚计算出来的 2×3 $2\times 3$的矩阵 WW $\pmb W$将样本映射到新的特征空间中，使用公式 yy=WWT×xx $\pmb y = \pmb W^T \times \pmb x$

transformed = matrix_w.T.dot(all_samples)
assert transformed.shape == (2,40), "The matrix is not 2x40 dimensional."

plt.plot(transformed[0,0:20], transformed[1,0:20], 'o', markersize=7, color='blue', alpha=0.5, label='class1')
plt.plot(transformed[0,20:40], transformed[1,20:40], '^', markersize=7, color='red', alpha=0.5, label='class2')
plt.xlim([-4,4])
plt.ylim([-4,4])
plt.xlabel('x_values')
plt.ylabel('y_values')
plt.legend()
plt.title('Transformed samples with class labels')

plt.show()

使用matplotlib.mlab库的PCA方法

我们一步步完成了PCA的工作，我们可以使用matplotlib库内置的PCA()类来更加便捷的实现PCA算法，原有的参考文档写的比较简单 (http://matplotlib.sourceforge.net/api/mlab_api.html#matplotlib.mlab.PCA) ，我们可以参加一个写的更好的文档：https://www.clear.rice.edu/comp130/12spring/pca/pca_docs.shtml

PCA()类的代码实现在下面的网址中：
https://sourcegraph.com/github.com/matplotlib/matplotlib/symbols/python/lib/matplotlib/mlab/PCA

PCA()类的属性

 Attrs:

    a : a centered unit sigma version of input a

    numrows, numcols: the dimensions of a

    mu : a numdims array of means of a

    sigma : a numdims array of atandard deviation of a

    fracs : the proportion of variance of each of the principal components

    Wt : the weight vector for projecting a numdims point or array into PCA space

    Y : a projected into PCA space

需要指出的是，PCA()类期待的输入类型是np.array()，其中假设numrows>numcols，所以我们必须转置我们的数据集。
matplotlib.mlab.PCA()在变换后保存所有的 dd $\pmb d$维输入（保存在属性PCA.Y中），并且已经是有序的（由于PCA分析在维度的重要性在主成分轴属于降序排序，因此我们看到，frac是单调递减值的。https://www.clear.rice.edu/comp130/12spring/pca/pca_docs.shtml），因此本例子中我们只需要画出前两列即可。

from matplotlib.mlab import PCA as mlabPCA

mlab_pca = mlabPCA(all_samples.T)

print('PC axes in terms of the measurement axes scaled by the standard deviations:\n', mlab_pca.Wt)

plt.plot(mlab_pca.Y[0:20,0],mlab_pca.Y[0:20,1], 'o', markersize=7, color='blue', alpha=0.5, label='class1')
plt.plot(mlab_pca.Y[20:40,0], mlab_pca.Y[20:40,1], '^', markersize=7, color='red', alpha=0.5, label='class2')

plt.xlabel('x_values')
plt.ylabel('y_values')
plt.xlim([-4,4])
plt.ylim([-4,4])
plt.legend()
plt.title('Transformed samples with class labels from matplotlib.mlab.PCA()')

plt.show()

PC axes in terms of the measurement axes scaled by the standard deviations:
     [[ 0.65043619  0.53023618  0.54385876]
     [-0.01692055  0.72595458 -0.68753447]
     [ 0.75937241 -0.43799491 -0.48115902]]

一步步实现PCA和使用matplotlib.mlab.PCA()二者的区别

当我们将转换的数据集绘制到新的二维子空间上时，可以观察到，我们一步步实现的方法和matplotlib.mlab.PCA（）类的画出的散点图看起来不一样。 这是由于matplotlib.mlab.PCA（）类在计算协方差矩阵之前将变量缩放为单位方差。这将/最终可能导致计算特征向量的方差不同，并影响变量对主要成分的贡献。

我们举个例子说明，一个变量以单位英寸（其中另一个变量以厘米为单位）。
然而，对于我们的假设例子，我们假设两个变量具有相同（任意）单位，所以我们跳过了缩放输入数据的步骤。

使用sklearn.decomposition的PCA方法来确认我们的结果

为了确保我们没有算错，我们将使用另一个库，该库在默认情况下不会重新缩放输入数据。
在这里，我们将使用scikit-learn机器学习库中的PCA类。 文档可以通过下面链接找到：
http://scikit-learn.org/stable/modules/generated/sklearn.decomposition.PCA.html.

为了方便起见，我们可以直接指定n_components参数来确定降维的数量。

 n_components : int, None or string

    Number 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

接下来，我们用.fit_transform()方法来进行降维操作。

from sklearn.decomposition import PCA as sklearnPCA

sklearn_pca = sklearnPCA(n_components=2)
sklearn_transf = sklearn_pca.fit_transform(all_samples.T)

plt.plot(sklearn_transf[0:20,0],sklearn_transf[0:20,1], 'o', markersize=7, color='blue', alpha=0.5, label='class1')
plt.plot(sklearn_transf[20:40,0], sklearn_transf[20:40,1], '^', markersize=7, color='red', alpha=0.5, label='class2')

plt.xlabel('x_values')
plt.ylabel('y_values')
plt.xlim([-4,4])
plt.ylim([-4,4])
plt.legend()
plt.title('Transformed samples with class labels from matplotlib.mlab.PCA()')

plt.show()

上面的图看起来像我们一步步得到的结果的一个镜像，这是因为特征向量的符号可以是正的或负的，因为特征向量被缩放到单位长度1，所以我们可以简单地将经变换的数据 ×(−1) $\times(-1)$来反转镜像。

sklearn_transf = sklearn_transf * (-1)

# sklearn.decomposition.PCA
plt.plot(sklearn_transf[0:20,0],sklearn_transf[0:20,1], 'o', markersize=7, color='blue', alpha=0.5, label='class1')
plt.plot(sklearn_transf[20:40,0], sklearn_transf[20:40,1], '^', markersize=7, color='red', alpha=0.5, label='class2')
plt.xlabel('x_values')
plt.ylabel('y_values')
plt.xlim([-4,4])
plt.ylim([-4,4])
plt.legend()
plt.title('Transformed samples via sklearn.decomposition.PCA')
plt.show()

# step by step PCA
plt.plot(transformed[0,0:20], transformed[1,0:20], 'o', markersize=7, color='blue', alpha=0.5, label='class1')
plt.plot(transformed[0,20:40], transformed[1,20:40], '^', markersize=7, color='red', alpha=0.5, label='class2')
plt.xlim([-4,4])
plt.ylim([-4,4])
plt.xlabel('x_values')
plt.ylabel('y_values')
plt.legend()
plt.title('Transformed samples step by step approach')
plt.show()

我们看上面的两个图，沿着主成分轴的数据分布看起来是一模一样的，只有数据中心略有不同，如果我们想要模拟scikit-learn库的PCA类产生的结果，我们可以从样本 XX $\pmb X$中减去均值向量，以便将坐标系原点处的数据居中（作者感谢Alexander Guth的建议）。将

transformed = matrix_w.T.dot(all_samples)

替换为

transformed = matrix_w.T.dot(all_samples - mean_vector)

翻译仓促，如果觉得有迷惑之处，可以参考原文。
原文地址：
http://sebastianraschka.com/Articles/2014_pca_step_by_step.html#generating-some-3-dimensional-sample-data