PCA算法实现

前言

PCA算法是数据降维中最常用的算法之一，利用PCA算法实现的数据降维能够有效减少算法运行时间和算法对硬件的消耗。本篇文章将使用python实现PCA算法，并将其应用于图像处理。

使用PCA算法实现降维

数据可视化

在算法实现之前，首先加载初始数据，并对初始数据进行可视化。这将有利于我们更好的了解PCA算法是如何将2D数据降维至1D数据的。在实际问题中，遇到的数据可能远远超过三维，为了能够实现数据可视化，常常将其降维至三维以下。二维原始数据的可视化实现代码及可视化图像如下所示：

plt.ion()
np.set_printoptions(formatter={'float': '{: 0.6f}'.format})
print('Visualizing example dataset for PCA.')
data = scio.loadmat('ex7data1.mat')
X = data['X']
plt.figure()
plt.scatter(X[:, 0], X[:, 1], facecolors='none', edgecolors='b', s=20)
plt.axis('equal')
plt.axis([0.5, 6.5, 2, 8])
input('Program paused. Press ENTER to continue')

算法实现

数据标准化
为了减少算法实现的复杂度，需要对训练样本进行标准化处理其标准化处理的公式为：
$x_{norm} = \frac{x-\mu}{\sigma}$
其中 $\mu$ 代表 $x$ 的均值， $\sigma$ 代表 $x$ 的标准差，其实现代码如下所示：

def feature_normalize(X):
    mu = np.mean(X, 0)
    sigma = np.std(X, 0, ddof=1)
    X_norm = (X - mu) / sigma
    return X_norm, mu, sigma

算法实现
PCA算法的实现可以调用库函数中的算法实现，其中 $U$ 就是主要成分，而 $S$ 是一个对角矩阵。具体实现代码如下所示：

def pca(X):
    
    U = np.zeros(n)
    S = np.zeros(n)
    #协方差矩阵
    sigma = np.dot(X.T, X) / m
    #full_matrices的取值是为0或者1，默认值为1，这时u的大小为(M,M)，v的大小为(N,N) 。
    #否则u的大小为(M,K)，v的大小为(K,N) ，K=min(M,N)。
    #compute_uv的取值是为0或者1，默认值为1，表示计算u,s,v。为0的时候只计算s。
    U, S, V = scipy.linalg.svd(sigma, full_matrices=True, compute_uv=True) 
    return U, S

如下代码所示，通过调用以上算法，我们已经找到了矩阵的特征向量并且绘制出该数据的主成分部分

X_norm, mu, sigma = feature_normalize(X)
# Run PCA
U,S = pca(X_norm)
draw_line(mu, mu + 1.5 * S[0] * U[:, 0])
draw_line(mu, mu + 1.5 * S[1] * U[:, 1])
# 特征向量
print('Top eigenvector: \nU[:, 0] = {}'.format(U[:, 0]))

PCA算法实现降维

以上，利用python实现了PCA算法，得到了其数据集的特征向量，现在可以使用PCA算法将高维数据投影到更低维的向量空间实现降维。以下代码通过PCA算法将数据集投影到特征向量所在的空间实现降维

数据投影
将 $n$ 维矩阵 $X$ 降维至 $k$ 维矩阵的计算公式如下所示，
$Z = U_k^{T}*X$ ，其中 $U_k$ 表示矩阵 $U$ 的前 $k$ 列向量。具体实现
代码如下所示：

def project_data(X, U, K): 
    Z = np.zeros((X.shape[0], K))
    Ureduce = U[:, np.arange(K)]
    Z = np.dot(X, Ureduce)
    return Z

压缩重现
实现数据降维之后，通过将其投影至原数据所在的高维空间，可以实现数据的压缩重现，数据压缩重现的公式如下所示：
$X_{approx} = U_k*Z$ ，其实现代码如下所示：

def recover_data(Z, U, K):
       X_approx = np.zeros((Z.shape[0], U.shape[0]))   
    Ureduce = U[:, np.arange(K)]
    X_approx = np.dot(Z, Ureduce.T)
    return X_approx

数据可视化
通过图像可视化展示数据降维和压缩重现，可以清楚的展示数据降维和数据压缩重现之间的关系，蓝色圆圈代表原来的数据，红色圆圈代表降维后的数据，实现代码和效果如下所示：

plt.figure()
plt.scatter(X_norm[:, 0], X_norm[:, 1], facecolors='none', edgecolors='b', s=20)
plt.axis('equal')
plt.axis([-4, 3, -4, 3])
K = 1
Z = project_data(X_norm, U, K)
print('Projection of the first example: {}'.format(Z[0]))
print('(this value should be about 1.481274)')
X_rec = recover_data(Z, U, K)
print('Approximation of the first example: {}'.format(X_rec[0]))
print('(this value should be about [-1.047419 -1.047419])')
plt.scatter(X_rec[:, 0], X_rec[:, 1], facecolors='none', edgecolors='r', s=20)
for i in range(X_norm.shape[0]):
    draw_line(X_norm[i], X_rec[i])

PCA算法的应用

PCA算法在图像处理上的应用

图像显示
以上，通过python实现了PCA算法，并利用PCA算法实现了数据降维和数据恢复。在这部分内容中会将PCA算法应用于人脸的图像数据集，人脸图像的矩阵共都包含着1024幅图像，32×32的显示在一个二维图像中，一共有5000个数据集，显示其中的前一百个数据如下所示：

def display_data(x):
    (m, n) = x.shape # m =5000,n=1024

    # 每幅图像都是32×32的灰度级
    example_width = np.round(np.sqrt(n)).astype(int) #width = 32
    example_height = (n / example_width).astype(int) #width = 32

    #计算显示的大小 dsiplay_rows = 70 display_cols = 72
    display_rows = np.floor(np.sqrt(m)).astype(int)
    display_cols = np.ceil(m / display_rows).astype(int)

    # 每幅图像之间的间隔
    pad = 1

    #设置显示矩阵大小 display_array = 2311*2311
    display_array = - np.ones((pad + display_rows * (example_height + pad),
                              pad + display_rows * (example_height + pad)))

    #用每一幅图像逐个替换显示矩阵的每个初始元素
    curr_ex = 0
    for j in range(display_rows):
        for i in range(display_cols):
            if curr_ex > m:
                break

           
            max_val = np.max(np.abs(x[curr_ex]))
            display_array[pad + j * (example_height + pad) + np.arange(example_height),
                          pad + i * (example_width + pad) + np.arange(example_width)[:, np.newaxis]] = \
                          x[curr_ex].reshape((example_height, example_width)) / max_val
            curr_ex += 1

        if curr_ex > m:
            break

    plt.figure()
    plt.imshow(display_array, cmap='gray', extent=[-1, 1, -1, 1])
    plt.axis('off')
data = scio.loadmat('ex7faces.mat')
X = data['X']
display_data(X[0:100])

运行结果如下所示：

PCA算法处理图像
在图像数据中使用PCA算法降维数据，首先要进行数据标准化，运行完PCA算法之后，注意到每一个矩阵 $U$ 主成份都是长度为 $n$ 的向量( $n = 1024$ )，可以将 $U$ 中向量重新设置为32×32的图像矩阵，实现图像数据的可视化，具体实现代码如下所示：

print('Running PCA on face dataset.\n(this might take a minute or two ...)')

# 标准化处理
X_norm, mu, sigma = feature_normalize(X)
U, S = pca(X_norm)
display_data(U[:, 0:36].T)

input('Program paused. Press ENTER to continue')