用几行代码实现sklearn中PCA的PCA.transform、PCA.fit_transform和PCA.inverse_transform功能

前情提要

想在控制器中用C语言实现降维或者降噪。本来以为PCA是由矩阵相乘后得到降维数据，所以苦苦寻找那个中间降维用的矩阵，但是无论怎么看，都不能从sklearn的介绍中找到这个矩阵。

然后，由降维前的数据和降维后的数据反推求到了一个中间矩阵，但发现将降维前的数据乘中间矩阵得到的结果，与sklearn算出来的误差挺大的。

最后，索性去看了PCA的源码，发现这个中间矩阵根本不存在， 这个中间矩阵还是存在的。

但是还是有所收获， 在学习sklearn中的PCA源码后，用几行代码提取出了里头的fit_transform和inverse_transform功能。

sklearn的PCA源码

Python代码

import numpy as np
from scipy import linalg
# svd_flip是为了保证输出的唯一性，也可以注释掉
#from sklearn.utils.extmath import svd_flip

# 降维后的维度是多少
n_components=2
# X是降维前的数据
means = np.mean(X, axis=0)
data = X - means
# 基于SVD分解协方差矩阵实现PCA算法
U, S, V = linalg.svd(data, full_matrices=False)
#U, V = svd_flip(U, V)

# 以下实现PCA.transform的功能
# 降维后的数据
X_transformed = data.dot(V[:n_components].T)

# 以下实现PCA.fit_transform的功能
# X_new = X * V = U * S * V^T * V = U * S
# 降维后的数据
X_new = U[:, :n_components] * S[:n_components]
# X_new 和 X_transformed 在数值上是相等的

# 以下实现PCA.inverse_transform功能
# 升维后的数据
X_inverse = X_new.dot(V[:n_components]) + means

实例测试

在这里我们用鸢尾花数据集看看效果。

from sklearn.datasets import load_iris

X = load_iris()['data']
print(X[:5])

显示结果为：

[[5.1 3.5 1.4 0.2]
[4.9 3. 1.4 0.2]
[4.7 3.2 1.3 0.2]
[4.6 3.1 1.5 0.2]
[5. 3.6 1.4 0.2]]

在带入上面的代码中：

from scipy import linalg
from sklearn.utils.extmath import svd_flip
import numpy as np

# 降维后的维度是多少
n_components=2
# X是降维前的数据
means = np.mean(X, axis=0)
data = X - means
U, S, V = linalg.svd(data, full_matrices=False)
U, V = svd_flip(U, V)

# 以下实现PCA.transform的功能
# 降维后的数据
X_transformed = data.dot(V[:n_components].T)
print('X_transformed:\n',X_transformed[:5],'\n')

# 以下实现PCA.fit_transform的功能
# X_new = X * V = U * S * V^T * V = U * S
# 降维后的数据
X_new = U[:, :n_components] * S[:n_components]
# X_new 和 X_transformed 在数值上是相等的
print('X_new:\n',X_new[:5],'\n')

# 以下实现PCA.inverse_transform功能
# 升维后的数据
X_inverse = X_new.dot(V[:n_components]) + means
print('X_inverse:\n',X_inverse[:5])

显示结果为：

X_transformed:
[[-2.68412563 0.31939725]
[-2.71414169 -0.17700123]
[-2.88899057 -0.14494943]
[-2.74534286 -0.31829898]
[-2.72871654 0.32675451]]

X_new:
[[-2.68412563 0.31939725]
[-2.71414169 -0.17700123]
[-2.88899057 -0.14494943]
[-2.74534286 -0.31829898]
[-2.72871654 0.32675451]]

X_inverse:
[[5.08303897 3.51741393 1.40321372 0.21353169]
[4.7462619 3.15749994 1.46356177 0.24024592]
[4.70411871 3.1956816 1.30821697 0.17518015]
[4.6422117 3.05696697 1.46132981 0.23973218]
[5.07175511 3.52655486 1.36373845 0.19699991]]

从结果可以看出，X_new 和 X_transformed 在数值上是相等的，说明中间矩阵是 V[:n_components].T ，只不过乘之前，要把数据中心化，即减去平均值。