Python中的主成分分析：分析错误

Question

我正在用 Python 实现一个 主成分分析（PCA）来进行人脸识别，但没有使用 numpy 或 OpenCV 中已经定义好的 PCA 方法。结果简直糟糕透了。

我查看了 OpenCV 的 文档，了解了算法的描述，但有些地方还是不太明白。

• 如果 X = {x1, x2, ..., xn}，那么我觉得这不是用来计算协方差矩阵的 X 矩阵吧？而是需要先减去均值等，正如前两步所解释的那样。
• 特征向量必须按照特征值从大到小的顺序排列。排序的时候是看特征值的绝对值吗？不过顺序不是主要问题，因为我会把所有特征向量画出来，所以可以相应地重新排列。我觉得我可能在分析上犯了错误。

我实现了以下内容：

MU = X.mean(axis=0)
for i in range(n):
    X[i,:] -= MU

S = (np.dot(X, X.T) / float(n))
#Hermitian (or symmetric) matrix.
eigenvalues, eigenvectors = numpy.linalg.eigh(S)
eigenvectors = numpy.dot(X.T, eigenvectors)
for i in range(n):
    eigenvectors[:,i] = normalize_vector(eigenvectors[:,i])

注意样本值是存储在行中，而不是列中。所以 X 的形状是 nxd，其中 n 是样本数量，d 是样本的维度。

上面的图片是参考图。第一张是均值，接下来的三张是最大的特征向量。下面的图片是我的结果。第一张是均值，接下来的都是特征向量，但它们似乎和参考结果不匹配。

使用 cv2.PCACompute(X, 6) 仍然能得到更好的结果。

Answer 1

第一个问题的答案可以在这里找到。

对于第二个问题，答案是你需要根据特征值的大小来排序特征向量。为此，你可以使用Python中的 np.argsort()，然后将这个顺序反转（argsort是从小到大排序的）：

indexes = np.argsort(eigenvalues)[::-1]
eigval = eigenvalues[indexes]
eigvec = eigenvectors[:,indexes]

---

检查你的代码后，我发现了几个问题：

1. 现在，你获取了所有的特征向量，忽略了 nb_components 参数，你应该只取你需要的特征向量。可以这样做：

   eigenvectors = eigenvectors[:,indexes][0:nb_components]

2. 在进行向量归一化时（在pca函数内），你使用了一个从0到 n 的循环，但如果你只需要，比如说，3个特征向量，你只会有3列。为了解决这个问题，应该从0循环到 nb_components。

除此之外，你的代码运行得很好。我尝试只取3个主成分，最后得到了6/6的结果。在我看来，显示特征向量的差异只是一个从浮点数转换为无符号8位整数（uint8）以使用imshow时的表示问题。

关于负特征值，这只是 eigh 的问题。特征值表示某个方向的方差，我们关心的是绝对值，但如果改变了符号，我们也需要改变“方向”（特征向量）。你可以通过将负特征值和它们对应的特征向量都乘以 -1.0 来解决这个问题（见这个链接）：

s = np.where(eigenvalues < 0)
eigenvalues[s] = eigenvalues[s] * -1.0
eigenvectors[:,s] = eigenvectors[:,s] * -1.0

你也可以使用 numpy.linalg.svd（文档）来解决这个问题，但这应该比 numpy.linalg.eigh 慢。

总之，这是我根据你的代码整理出来的代码（为了简洁，我删除了所有注释）：

def pca(X, nb_components=0):
    [n,d] = X.shape
    if (nb_components <= 0) or (nb_components>6):
        nb_components = n

    MU = X.mean(axis=0)
    for i in range(n):
        X[i,:] -= MU

    S = (np.dot(X, X.T) / float(n))

    eigenvalues, eigenvectors = np.linalg.eigh(S)

    s = np.where(eigenvalues < 0)
    eigenvalues[s] = eigenvalues[s] * -1.0
    eigenvectors[:,s] = eigenvectors[:,s] * -1.0

    indexes = np.argsort(eigenvalues)[::-1]
    eigenvalues = eigenvalues[indexes]
    eigenvectors = eigenvectors[:,indexes][:,0:nb_components]

    eigenvectors = np.dot(X.T, eigenvectors)

    for i in range(nb_components):
        eigenvectors[:,i] = normalize_vector(eigenvectors[:,i])

    return (eigenvalues, eigenvectors, MU)