矩阵分解: SVD-PCA

矩阵分解

矩阵分解（Decomposition Factorization）是将矩阵拆解为若干个矩阵的相乘的过程。在数值分析中，常常被用来实现一些矩阵运算的快速算法，在机器学习领域有非常重要的作用。有的推荐系统采用SVD算法来实现整套系统中的矩阵分解过程。

SVD算法即为奇异值分解法，相对于矩阵的特征值分解法，它可以对非方阵形式的矩阵进行分解，将一个矩阵A分解为如下形式：

$A=UΣV^T$

其中：

A代表需要被分解的矩阵，设其维度是 $m×n$
U矩阵是被分解为的3个矩阵之一，它是一个 $m×m$ 的方阵，构成这个矩阵的向量是正交的，被称为左奇异向量
Σ是一个 $m×n$ 的向量，它的特点是除了对角线中的元素外，其余元素都为0
V是一个 $n×n$ 的方阵，它的转置也是一个方阵，与U矩阵类似，构成这个矩阵的向量也是正交的，被称为右奇异向量

基于Numpy实现SVD

import numpy as np

matrix = np.array([[1, 2], [3, 4]])

another_matrix = np.dot(matrix, matrix.T)

U, s, V = np.linalg.svd(another_matrix)  # 使用奇异值分解法将矩阵进行分解，得到3个子矩阵U,s,V

# 在s矩阵的基础上，生成S矩阵为
S = np.array([[s[0], 0], [0, s[1]]])

# 利用生成的USV三个矩阵，可以重建回原来的矩阵another_matrix
np.dot(U, np.dot(S, V))

其中得到的结果如下：

s: array([29.86606875,  0.13393125])
S: array([[29.86606875,  0.        ],
       [ 0.        ,  0.13393125]])
       
U：array([[-0.40455358, -0.9145143 ],
       [-0.9145143 ,  0.40455358]])

V：array([[-0.40455358, -0.9145143 ],
       [-0.9145143 ,  0.40455358]])
       
重建：array([[ 5., 11.],
       [11., 25.]])

在上面的代码片段中，s向量表示的是分解后的Σ矩阵中对角线上的元素，所以在这里面引入了一个S矩阵，将s向量中的元素放置在这个矩阵中，用以验证分解后的矩阵重建回原先的矩阵A的过程。

基于SVD实现PCA算法

# 零均值预处理
def zero_centered(data):
    matrix_mean = np.mean(data, axis=0)
    return data - matrix_mean

def pca_eig(data, n):
    new_data = zero_centered(data)
    cov_mat = np.dot(new_data.T, new_data)
    eig_values, eig_vectors = np.linalg.eig(np.mat(cov_mat))
    value_indices = np.argsort(eig_values)  # 特征值从小到大排序
    n_vectors = eig_vectors[:, value_indices[-1: -(n+1): -1]]  # 最大的n个特征值对应的特征向量
    return new_data * n_vectors  # 返回低维空间的数据 x * v

def pca_svd(data, n):
    new_data = zero_centered(data)
    cov_mat = np.dot(new_data.T, new_data)  # 协方差矩阵
    U, s, V = np.linalg.svd(cov_mat)
    pc = np.dot(new_data, U)  # 返回矩阵的第1个列向量即是降维后的结果 
    return pc[:, 0]

def unit_test():
    data = np.array(
        [[2.5, 2.4], [0.5, 0.7], [2.2, 2.9], [1.9, 2.2], [3.1, 3.0], [2.3, 2.7],
        [2, 1.6], [1, 1.1], [1.5, 1.6], [1.1, 0.9]]
    )
    result_eig = pca_eig(data, 1)  # 使用常规的特征矩阵分解，将二维数据降到一维
    print(result_eig)
    result_svd = pca_svd(data, 1)  # 使用奇异值分解法将协方差矩阵分解，得到降维结果
    print(result_svd)
    

if __name__ == '__main__':
    unit_test()

[[-0.82797019]
 [ 1.77758033]
 [-0.99219749]
 [-0.27421042]
 [-1.67580142]
 [-0.9129491 ]
 [ 0.09910944]
 [ 1.14457216]
 [ 0.43804614]
 [ 1.22382056]]
 
[-0.82797019  1.77758033 -0.99219749 -0.27421042 -1.67580142 -0.9129491  0.09910944  1.14457216  0.43804614  1.22382056]

可以看到，数据已经从二维变为一维了，这两个PCA算法的计算结果是相同的。其中pca_eig()函数使用常规的特征值分解方法来求解，读者可以参照前面讲述的PCA算法过程来理解这段代码。pca_svd()函数是使用奇异值分解法来求解的。

下面简要阐述一下PCA算法中奇异值分解的步骤:

第一步，PCA算法中得到样本的协方差矩阵是经过零均值化处理的：

$C=X^T X$

其中，X是经过中心化处理后的样本矩阵，一个矩阵与其转置矩阵相乘的结果是一个对称矩阵，所以C是对称矩阵，将其进行奇异值分解后可以表示为：

$C=UΣV^T$
第二步，将经过中心化的样本矩阵X进行奇异值分解，可以得到：

$X=UΣV^T$

$X^TX \\ = (UΣV^T) ^T (UΣV^T) \\ = VΣ^T U^T UΣV^T \\ = VΣ^2 V^T$

奇异矩阵V中的列对应着PCA算法主成分中的主方向，因此可以得到主成分为：

$XV=UΣV^TV =UΣ$

SVD与PCA等价，所以PCA问题可以转化为SVD问题求解，那转化为SVD问题有什么好处？

有三点：

一般 $X$ 的维度很高， $X^TX$ 的计算量很大
方阵的特征值分解计算效率不高
SVD除了特征值分解这种求解方式外，还有更高效且更准确的迭代求解法，避免了 $X^TX$ 的计算

其实，PCA只与SVD的右奇异向量的压缩效果相同：

如果取 $V$ 的前 $k$ 行作为变换矩阵 $P_{k×n}$ ，则 $Y_{k×m}=P_{k×n}X_{n×m}$ ，起到压缩行即降维的效果
如果取 $U$ 的前 $d$ 行作为变换矩阵 $P_{d×m}$ ，则 $Y_{n×d}=P_{n×m}X_{m×d}$ ，起到压缩列即去除冗余样本的效果

参考

stats.stackexchange.com/questions/1…