用 Numpy 实现 PCA
Posted
技术标签:
【中文标题】用 Numpy 实现 PCA【英文标题】:Implementing PCA with Numpy 【发布时间】:2020-02-28 04:30:46 【问题描述】:我想用类似于sklearn
中的类来实现 PCA。
我的寻找具有 k 个主成分的 PCA 的算法如下:
计算样本均值并转换数据集,使其以原点为中心。 计算新翻译集的协方差矩阵。 找出特征值和特征向量,按降序排列。 将数据集投影到由前 k 个特征向量跨越的向量空间。import numpy as np
class MyPCA:
def __init__(self, n_components):
self.n_components = n_components
def fit_transform(self, X):
"""
Assumes observations in X are passed as rows of a numpy array.
"""
# Translate the dataset so it's centered around 0
translated_X = X - np.mean(X, axis=0)
# Calculate the eigenvalues and eigenvectors of the covariance matrix
e_values, e_vectors = np.linalg.eigh(np.cov(translated_X.T))
# Sort eigenvalues and their eigenvectors in descending order
e_ind_order = np.flip(e_values.argsort())
e_values = e_values[e_ind_order]
e_vectors = e_vectors[e_ind_order]
# Save the first n_components eigenvectors as principal components
principal_components = np.take(e_vectors, np.arange(self.n_components), axis=0)
return np.matmul(translated_X, principal_components.T)
但是,当在 Iris 数据集上运行时,此实现产生的结果与 sklearn 的结果大不相同,并且结果并未表明数据中存在三个不同的组:
from sklearn import datasets
from sklearn.decomposition import PCA
import matplotlib.pyplot as plt
def plot_pca_results(pca_class, dataset, plot_title):
X = dataset.data
y = dataset.target
y_names = dataset.target_names
pca = pca_class(n_components=1)
B = pca.fit_transform(X)
B = np.concatenate([B, np.zeros_like(B)], 1)
scatter = plt.scatter(B[:, 0], B[:, 1], c=y)
scatter_objects, _ = scatter.legend_elements()
plt.title(plot_title)
plt.legend(scatter_objects, y_names, loc="lower left", title="Classes")
plt.show()
dataset = datasets.load_iris()
plot_pca_results(MyPCA, dataset, "Iris - my PCA")
plot_pca_results(PCA, dataset, "Iris - Sklearn")
造成这种差异的原因可能是什么?我的方法在哪里,或者我的计算不正确?
【问题讨论】:
要对 PCA 的数据进行归一化,不仅需要 mean=0,还需要 stdev=1。类似translated_X = (X - np.mean(X, axis=0))/np.std(X, axis=0)
@krubo 即使在我应用了这个归一化/z-score 之后,结果仍然是关闭的。
【参考方案1】:
比较两种方法
问题在于没有标准化数据和提取特征向量(主轴)。该函数比较两种方法。
import numpy as np
from sklearn import datasets
from sklearn.decomposition import PCA
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
def pca_comparison(X, n_components, labels):
"""X: Standardized dataset, observations on rows
n_components: dimensionality of the reduced space
labels: targets, for visualization
"""
# numpy
# -----
# calculate eigen values
X_cov = np.cov(X.T)
e_values, e_vectors = np.linalg.eigh(X_cov)
# Sort eigenvalues and their eigenvectors in descending order
e_ind_order = np.flip(e_values.argsort())
e_values = e_values[e_ind_order]
e_vectors = e_vectors[:, e_ind_order] # note that we have to re-order the columns, not rows
# now we can project the dataset on to the eigen vectors (principal axes)
prin_comp_evd = X @ e_vectors
# sklearn
# -------
pca = PCA(n_components=n_components)
prin_comp_sklearn = pca.fit_transform(X)
# plotting
if n_components == 3:
fig = plt.figure(figsize=(10, 5))
ax = fig.add_subplot(121, projection='3d')
ax.scatter(prin_comp_sklearn[:, 0],
prin_comp_sklearn[:, 1],
prin_comp_sklearn[:, 1],
c=labels)
ax.set_title("sklearn plot")
ax = fig.add_subplot(122, projection='3d')
ax.scatter(prin_comp_evd[:, 0],
prin_comp_evd[:, 1],
prin_comp_evd[:, 2],
c=labels)
ax.set_title("PCA using EVD plot")
fig.suptitle(f"Plots for reducing to n_components-D")
plt.show()
elif n_components == 2:
fig, ax = plt.subplots(1, 2, figsize=(10, 5))
ax[0].scatter(prin_comp_sklearn[:, 0],
prin_comp_sklearn[:, 1],
c=labels)
ax[0].set_title("sklearn plot")
ax[1].scatter(prin_comp_evd[:, 0],
prin_comp_evd[:, 1],
c=labels)
ax[1].set_title("PCA using EVD plot")
fig.suptitle(f"Plots for reducing to n_components-D")
plt.show()
elif n_components == 1:
fig, ax = plt.subplots(1, 2, figsize=(10, 5))
ax[0].scatter(prin_comp_sklearn[:, 0],
np.zeros_like(prin_comp_sklearn[:, 0]),
c=labels)
ax[0].set_title("sklearn plot")
ax[1].scatter(prin_comp_evd[:, 0],
np.zeros_like(prin_comp_evd[:, 0]),
c=labels)
ax[1].set_title("PCA using EVD plot")
fig.suptitle(f"Plots for reducing to n_components-D")
plt.show()
return prin_comp_sklearn, prin_comp_evd[:, :n_components]
加载数据集、预处理和运行实验:
dataset = datasets.load_iris()
X = dataset.data
mean = np.mean(X, axis=0)
# this was missing in your implementation
std = np.std(X, axis=0)
X_std = (X - mean) / std
for n in [3, 2, 1]:
pca_comparison(X_std, n, dataset.target)
结果
3D 图有点杂乱,但是如果您查看 2D 和 1D 的情况,您会发现如果我们将第一个主成分乘以 -1,这些图是相同的; scikit-learn PCA 实现在底层使用奇异值分解,这将提供非唯一的解决方案 (see here)。
测试:
使用来自here 的flip_signs()
函数
def flip_signs(A, B):
"""
utility function for resolving the sign ambiguity in SVD
http://stats.stackexchange.com/q/34396/115202
"""
signs = np.sign(A) * np.sign(B)
return A, B * signs
for n in [3, 2, 1]:
sklearn_pca, evd = pca_comparison(X_std, n, dataset.target)
assert np.allclose(*flip_signs(sklearn_pca, evd))
实施中的问题:
-
如果我们看一下iris数据集中的数据尺度,数据是不同尺度的。这表明我们应该标准化数据 (Read here for more)
引用上述答案的一部分:
由@ttnphns 续
什么时候更愿意对相关性(即 z 标准化变量)而不是协方差(即中心变量)进行 PCA(或因子分析或其他类似类型的分析)?
当变量是不同的测量单位时。很清楚
...
-
为了使用
e_values, e_vectors = np.linalg.eigh(X_cov)
获得主轴,您应该提取e_vectors
(documentation) 的列。您正在提取行。
【讨论】:
以上是关于用 Numpy 实现 PCA的主要内容,如果未能解决你的问题,请参考以下文章
PCA详解-并用scikit-learn实现PCA压缩红酒数据集