用numpy进行二维PCA直线拟合

Question

我正在尝试使用numpy实现2D PCA。 代码相当简单：

import numpy as np

n=10
d=10
x=np.linspace(0,10,n)
y=x*d

covmat = np.cov([x,y])
print(covmat)

eig_values, eig_vecs = np.linalg.eig(covmat)
largest_index = np.argmax(eig_values)
largest_eig_vec = eig_vecs[largest_index]

协方差矩阵为：

[[   11.31687243   113.16872428]
 [  113.16872428  1131.6872428 ]]

然后我有一个简单的帮助方法，可以绘制一条围绕给定中心的线（作为一系列点），在给定方向上。这是要由pyplot使用的，因此我正在准备x坐标和y坐标的分别列表。

def plot_line(center, dir, num_steps, step_size):
    line_x = []
    line_y = []
    for i in range(num_steps):
        dist_from_center = step_size * (i - num_steps / 2)
        point_on_line = center + dist_from_center * dir
        line_x.append(point_on_line[0])
        line_y.append(point_on_line[1])
    return (line_x, line_y)

最后是剧情设置：

lines = []
mean_point=np.array([np.mean(x),np.mean(y)])
lines.append(plot_line(mean_point, largest_eig_vec, 200, 0.5))

import matplotlib.pyplot as plt

fig = plt.figure()
ax = fig.add_subplot(111)

ax.scatter(x,y, c="b", marker=".", s=10
           )
for line in lines:
    ax.plot(line[0], line[1], c="r")

ax.scatter(mean_point[0], mean_point[1], c="y", marker="o", s=20)

plt.axes().set_aspect('equal', 'datalim')
plt.show()

很不幸，主成分分析似乎无法工作。 下面是图表：

很抱歉，我不知道出了什么问题。

• 我手动计算了协方差 -> 得到相同的结果。
• 我检查了其他特征值 -> 与红线垂直。
• 我使用方向(1,10)测试plot_line。它与我的点完全对齐：

最终的图表显示，由pca拟合的线是正确的结果，只是在y轴上镜像。

事实上，如果我更改特征向量的x坐标，该线就会完美拟合：

显然这是一个基本问题。不知何故我误解了如何使用pca。

我的错误在哪里？ 在线资源似乎完全描述了我实现的PCA。 我不相信我必须要在y轴上彻底镜像我的线性拟合。一定是其他原因。

- lhk

Answer 1

你的错误在于你提取了特征向量数组的最后一行。但是特征向量是由np.linalg.eig返回的特征向量数组的列而不是行。从文档中可以看到:

[...] 数组 a、w 和 v 满足方程 dot(a[:,:], v[:,i]) = w[i] * v[:,i]，其中 i 为每个元素。

其中 a 是应用了 np.linalg.eig 的数组，w 是特征值的一维数组，v 是特征向量的二维数组。因此，列 v[:, i] 是特征向量。

在这个简单的二维情况下，由于两个特征向量是相互正交的（因为我们从一个对称矩阵开始），并且长度为单位（因为np.linalg.eig按这种方式归一化它们），因此特征向量数组具有以下两种形式之一:

[[ cos(t)  sin(t)]
 [-sin(t)  cos(t)]]

或者

[[ cos(t)  sin(t)]
 [ sin(t) -cos(t)]]

针对某些实数 t，在第一种情况下，读取第一行（例如）而不是第一列将在 [cos(t), -sin(t)] 的位置上给出 [cos(t), sin(t)]。这可以解释您看到的表面反射现象。

请更换该行。

largest_eig_vec = eig_vecs[largest_index]

使用

largest_eig_vec = eig_vecs[:, largest_index]

并且您应该获得预期的结果。

- Mark Dickinson

---

感谢您的快速回答。 - lhk