前言
线性代数是高等数学里面非常重要的一门课程,可惜在学校的时候是一种自底向上的学习方式,并不知道学出来有什么用,以致彻底沦为应试教育。后来在工作中接触了机器学习,才真正看到了“数学之美”,才发现线性代数是多么的优雅多么的有用。
今天我们来看看下线性代数中非常重要的一个知识点奇异值分解SVD Singular value decomposition。SVD在数据科学当中非常有用,其常见的应用包括: - 自然语言处理中的Latent Semantic Analysis - 推荐系统中的Collaborative Filtering - 降维常用套路Principal Component Analysis
LSA已经在前文中有所讲解,CF的话后面在推荐系统的专题中来写,今天主要聊聊PCA,以及SVD在PCA中的重要作用。同样延续我们“手撕”的传统,使用numpy来理解其中的原理。
PCA
Principal component analysis即主成分分析,是机器学习中一种非常常用的降维方式。其发源也是源自于早期的计算机处理能力有限,当数据样本的维度很高时,预先去除掉数据中的一些冗余信息和噪声(降维),使得数据变得更加简单高效,节省时间和成本。
在深度学习时代,更强调的是原始数据的直接输入,再通过神经网络来做降维工作,最典型是场景就是计算机视觉,直接输入原始图片像素信息,通过CNN卷积层、MaxPooling层来进行降维。因此PCA逐渐开始淡出人们的视线,通常是作为一种数据可视化的手段(二维图表无法展示多维的数据样本)。
其实,在深度学习目前尚未全面攻克的结构化数据领域,PCA仍然有较多的用,其数据处理的思路依然值得我们去学习揣摩。
PCA正常解法
PCA算法的本质,其实就是找到一些投影方向,使得数据在这些投影方向上的方差最大,且这些投影方向是相互正交的。找到新的正交基后,计算原始数据在这些正交基上投影的方差,方差越大,就说明对应正交基上包含了更多的信息量。
关于原始数据的方差,最好的一个工具就是协方差矩阵了。协方差矩阵的特征值越大,对应的方差也就越大,在对应的特征向量上投影的信息量就越大。因此,我们如果将小特征值对应方向的数据删除,就可以达到降维的目的。因此,在数学上,我们可以把问题转化为求原始数据的协方差矩阵,然后计算协方差矩阵的特征值与特征向量。
对于广大程序员来说,学习机器学习最重要的一个坎还是数学。很多实际的代码其实是公式推导后的结果的代码实现,如果没有理清公式推导的过程,那么最后肯定是一头雾水。所以,克服心中的恐惧,翻出压在箱底的《线性代数》,我们上。
首先,求原始数据X的协方差矩阵C,将原始矩阵中心化后,做如下操作
接着,由于协方差矩阵C是方阵,就可以通过特征分解的方式来求C的特征值和特征向量。
最后,选择最大的k个特征值进行保留,求X的k阶PCA(X右乘k阶特征向量)
用SVD来解PCA
根据上面的推导,我们已经可以对矩阵X做PCA了。同学们可能要问了,这跟SVD有什么关系呢?
工程化思维强的同学应该已经想到了,这种纯数学的解法,在实际工程实践中有以下问题: - 在数据量很大时,把原始矩阵进行转置求协方差矩阵,然后再进行特征值分解是一个非常慢的过程。 - 稳定性问题。可以看到X转置乘以X,如果矩阵有非常小的数,很容易在平方中丢失。
工业界中,还是“唯快不破,唯稳不破”。我们知道,奇异值分解相对特征分解,有个很大的优势就是不要求原始矩阵是方阵。这非常符合现实生活中的数据。因此,有大神想到,是否可以用svd来解PCA?推导如下:
我们根据协方差矩阵的公式,把X按照奇异值分解展开,注意后面应用到了一个酋矩阵(unitary)的特性:
看到最后的结果,是否跟上面的
很像?没错。协方差矩阵C的特征值和X的奇异值有以下关系
而C的特征向量即为X的SVD分解后的V向量, 则参考PCA正常解法,X的k阶PCA即为_X右乘k阶V向量_。因此这种方式求PCA,只需要把原始矩阵做一次SVD分解即可,不用转置,不用求协方差矩阵。事实上,在Scikit Learn
等机器学习框架中,就是用的SVD来做PCA。
用numpy来验证
numpy原始解法求PCA
接下来,我们用numpy来验证这种思路。首先是PCA的标准解法: 随机模拟一个原始数据矩阵,5个样本,3个特征:
import numpy as np
X = np.random.rand(5,3)
'''
array([[0.86568791, 0.73022945, 0.17982869],
[0.07201287, 0.99358411, 0.84389196],
[0.61267696, 0.08867997, 0.11770573],
[0.16898969, 0.3093472 , 0.9010064 ],
[0.43840269, 0.97250927, 0.64897872]])
'''
将矩阵中心化,即减去均值:
X_new = X - np.mean(X, axis=0)
'''
array([[ 0.43413389, 0.11135945, -0.35845361],
[-0.35954115, 0.37471411, 0.30560966],
[ 0.18112294, -0.53019003, -0.42057657],
[-0.26256433, -0.3095228 , 0.3627241 ],
[ 0.00684866, 0.35363927, 0.11069642]])
'''
# 确保结果正确,即转换后均值为0
np.allclose(X_new.mean(axis=0), np.zeros(X_new.shape[1]))
求X_new的协方差矩阵C
C = np.dot(X_new.T, X_new) / (X_new.shape[0] - 1)
'''
C
array([[ 0.10488363, -0.02467955, -0.10903811],
[-0.02467955, 0.16369454, 0.05611495],
[-0.10903811, 0.05611495, 0.13564834]])
'''
求C的特征值和特征向量,这里用的是numpy的特征分解函数
eig_vals, eig_vecs = np.linalg.eig(C)
'''
eig_vals
array([0.26474535, 0.00779743, 0.13168373])
eig_vecs
array([[-0.53801107, 0.72610049, -0.42816139],
[ 0.50584138, -0.12820944, -0.85304562],
[ 0.67429117, 0.67552974, 0.29831358]])
'''
求X的PCA结果,就是X右乘k阶特征向量。这里k还是取的3。
X_pca = np.dot(X_new, eig_vecs)
'''
array([[-0.41894072, 0.05880142, -0.38780564],
[ 0.58905292, -0.10265648, -0.07453908],
[-0.64922927, -0.08462316, 0.24926273],
[ 0.22927473, 0.09406657, 0.4846625 ],
[ 0.24984234, 0.03441165, -0.27158052]])
'''
numpy的SVD求PCA
首先,直接求X_new的SVD,同样使用numpy的函数
# 注意这里的Vh其实是公式中的VT
U, Sigma, Vh = np.linalg.svd(X_new, full_matrices=False, compute_uv=True)
'''
U
array([[-0.40710685, 0.53434046, 0.33295236, 0.59843699, -0.28241844],
[ 0.57241386, 0.10270414, -0.58127362, 0.5471169 , -0.15677466],
[-0.63089041, -0.34344824, -0.47916326, 0.32579597, 0.38507162],
[ 0.22279838, -0.66779531, 0.53263483, 0.46842625, 0.03587876],
[ 0.24278501, 0.37419895, 0.19484969, 0.13026931, 0.86376738]])
Sigma
array([1.02906823, 0.72576506, 0.17660612])
Vh
array([[-0.53801107, 0.50584138, 0.67429117],
[ 0.42816139, 0.85304562, -0.29831358],
[ 0.72610049, -0.12820944, 0.67552974]])
'''
我们来根据上面的公式,确认下eig_vals和S的关系,注意在numpy的实现中,特征值和奇异值的排序是不同的
np.allclose(eig_vals, np.square(S) / (X_new.shape[0] - 1))
'''
eig_vals
array([0.26474535, 0.00779743, 0.13168373])
np.square(S) / (X_new.shape[0] - 1)
array([0.26474535, 0.13168373, 0.00779743])
'''
从结果看出,确实跟公式是一致的。 接下来用SVD求PCA就简单了,直接右乘V即可。
# 注意Vh是公式中的VT,因此V=Vh.T
X_pca_svd = np.dot(X_new, Vh.T)
# X_pca_svd = np.dot(U, np.diag(Sigma))
'''
X_pca_svd
array([[-0.41894072, 0.38780564, 0.05880142],
[ 0.58905292, 0.07453908, -0.10265648],
[-0.64922927, -0.24926273, -0.08462316],
[ 0.22927473, -0.4846625 , 0.09406657],
[ 0.24984234, 0.27158052, 0.03441165]])
'''
求出结果后,正当我们信心满满的对比一下X_pca
和X_pca_svd
,以为大功告成打完收工时,却发现二者是不一致的。WTF?
结果分析
仔细研究下X_pca
和X_pca_svd
的结果,可以看出,排除特征值和奇异值的排序导致的列向量顺序不同外,部分列向量的绝对值相同但正负不同。
问题出在哪里?我们搬出Scikit Learn
,再来算一次PCA:
from sklearn.decomposition import PCA
pca = PCA(3)
pca.fit_transform(X) # sklearn自动处理去均值化
'''
array([[ 0.41894072, -0.38780564, -0.05880142],
[-0.58905292, -0.07453908, 0.10265648],
[ 0.64922927, 0.24926273, 0.08462316],
[-0.22927473, 0.4846625 , -0.09406657],
[-0.24984234, -0.27158052, -0.03441165]])
'''
嗯,也是绝对值相同但正负不同。都说Scikit Learn
的PCA就是SVD做的,难道是骗人的? 好在代码不会骗人,我们直接翻出源码。
通过研究Scikit Learn
的源码svd_flip@scikit-learn/extmath.py找到了答案: SVD奇异值分解的结果是唯一的,但是分解出来的U矩阵和V矩阵的正负可以不是唯一,只要保证它们乘起来是一致的就行。因此,sklearn为了保证svd分解结果的一致性,它们的方案是:保证U矩阵的每一行(u_i
)中,绝对值最大的元素一定是正数,否则将u_i
转成-u_i
,并将相应的v_i
转成-v_i
已保证结果的一致。
这又是数学与工程的问题了。在数学上,几种结果都是正确的。但是在工程上,有个很重要的特性叫幂等性(Idempotence)。
Methods can also have the property of “idempotence” in that (aside from error or expiration issues) the side-effects of N > 0 identical requests is the same as for a single request.
这是源自于HTTP规范中的一个概念,可以引申至各种分布式服务的设计当中,即:高质量的服务,一次请求和多次请求,其副作用(结果)应当是一致的。Scikit Learn
正是通过svd_flip
这个函数,把一个数学上并不幂等的操作,转化成了幂等的服务,其设计之讲究可见一斑。
总结
本文通过公式推导和numpy代码实战,展示了PCA的正常解法,以及工业界常用的SVD解法,并最后引申至数学和实现的一些探讨。“part-science, part-art”,这就是我最喜爱的机器学习之道