假设有一个矩阵B，其大小为500*1000 double（这里，500代表观察数，1000代表特征数）。
sigma是B的协方差矩阵，D是一个对角矩阵，其对角元素是sigma的特征值。假设A是协方差矩阵sigma的特征向量。
我有以下问题：
我需要选择最大的特征值对应的第一个k = 800特征向量来对选择的特征进行排序。最终的矩阵命名为Aq。如何在 MATLAB 中做到这一点？
这些被选中的特征向量是什么意思？
一旦我计算出Aq，最终矩阵Aq 的大小似乎是1000*800 double。 500的时间点/观测信息已经消失。对于最终的矩阵Aq，矩阵Aq 中的值1000 现在代表什么？另外，矩阵Aq中的值800现在代表什么？
【问题讨论】：
标签： matlab machine-learning pca data-analysis
我假设您从 eig 函数中确定了特征向量。以后我会向你推荐的是使用eigs 函数。这不仅会为您计算特征值和特征向量，还会为您计算 k largest 特征值及其关联的特征向量。这可以节省计算开销，您不必计算矩阵的所有特征值和相关的特征向量，因为您只需要一个子集。您只需将数据的协方差矩阵提供给eigs，它就会为您返回k 的最大特征值和特征向量。
现在，回到您的问题，您所描述的最终是Principal Component Analysis。这背后的机制是计算数据的协方差矩阵并找到计算结果的特征值和特征向量。众所周知，不推荐这样做，因为计算大型矩阵的特征值和特征向量时数值不稳定。现在最规范的方法是通过Singular Value Decomposition。具体来说，V 矩阵的列给出了协方差矩阵的特征向量或主成分，相关的特征值是矩阵S 对角线中产生的奇异值的平方根。
请参阅 Cross Validated 上的这篇内容丰富的帖子，了解为什么这是首选：
https://stats.stackexchange.com/questions/79043/why-pca-of-data-by-means-of-svd-of-the-data
我还会插入另一个链接，讨论为什么在主成分分析中使用奇异值分解背后的理论：
现在让我们一次一个地回答你的问题。
MATLAB 以 未排序 的方式生成特征值和特征向量的相应排序。如果您希望在给定eig（在您的示例中为 800）的输出中选择最大的 k 特征值和相关的特征向量，您需要按降序对特征值进行排序，然后重新排列从eig 生成的特征向量矩阵的列，然后选择第一个k 值。
我还应该注意，使用eigs 并不能保证排序顺序，因此当涉及到它时，您也必须明确地对它们进行排序。
在 MATLAB 中，执行我们上面描述的操作看起来像这样：
sigma = cov(B);
[A,D] = eig(sigma);
vals = diag(D);
[~,ind] = sort(abs(vals), 'descend');
Asort = A(:,ind);
值得注意的是，您对特征值的绝对值进行排序，因为缩放后的特征值本身也是特征值。这些量表也包括负数。这意味着，如果我们有一个特征值为 -10000 的组件，这很好地表明该组件对您的数据具有重要意义，如果我们仅根据数字本身进行排序，则它位于较低的排名附近.
第一行代码找到了B 的协方差矩阵，尽管你说它已经存储在sigma 中，但是让我们让它可重现。接下来，我们找到协方差矩阵的特征值和相关的特征向量。注意特征向量矩阵A的每一列代表一个特征向量。具体来说，A 的第 ith 列/特征向量对应于 D 中看到的第 ith 特征值。
但是，特征值在对角矩阵中，因此我们使用 diag 命令提取对角线，对它们进行排序并计算它们的顺序，然后重新排列 A 以遵守此顺序。我使用sort 的第二个输出，因为它告诉您未排序结果中的每个值将出现在排序结果中的位置。这是我们需要重新排列特征向量矩阵A 的列的顺序。必须选择'descend' 作为标志，以便最大的特征值和相关的特征向量首先出现，就像我们之前讨论的那样。
然后您可以通过以下方式提取出第一个 k 最大的向量和值：
k = 800;
Aq = Asort(:,1:k);
众所周知，协方差矩阵的特征向量等于主成分。具体来说，第一个主成分（即最大特征向量和相关的最大特征值）为您提供数据中最大可变性的方向。之后的每个主成分都会给您带来递减性质的可变性。还需要注意的是，每个主成分彼此正交。
这是来自维基百科的二维数据的一个很好的例子：
我从维基百科关于主成分分析的文章中提取了上面的图像，我将你链接到上面。这是根据以(1,3) 为中心的二元高斯分布分布的样本散点图，在大致(0.878, 0.478) 方向上的标准差为 3，在正交方向上的标准差为 1。标准差为 3 的分量是第一个主分量，而正交的分量是第二个分量。显示的向量是协方差矩阵的特征向量，按相应特征值的平方根缩放，并移动以使其尾部位于均值处。
现在让我们回到您的问题。我们查看k 最大特征值的原因是执行dimensionality reduction 的一种方式。本质上，您将执行数据压缩，您将获取更高维度的数据并将它们投影到更低维度的空间。您在投影中包含的主要成分越多，它就越类似于原始数据。它实际上在某个点开始逐渐减少，但前几个主要成分允许您忠实地重建大部分数据。
我在过去偶然发现的这篇很棒的 Quora 帖子中发现了一个执行 PCA（或者更确切地说是 SVD）和数据重建的出色视觉示例。
您可以使用此矩阵将高维数据重新投影到低维空间。 1000 的行数仍然存在，这意味着您的数据集中最初有 1000 个要素。 800 是数据的降维量。将此矩阵视为从特征的原始维度 (1000) 到其降维 (800) 的转换。
然后，您可以将此矩阵与重建原始数据结合使用。具体来说，这将为您提供原始数据的近似值，且误差最小。在这种情况下，您不需要使用所有主成分（即仅使用 k 最大向量），您可以使用比以前更少的信息来创建数据的近似值。
如何重建数据非常简单。让我们先用完整的数据来谈谈正向和反向操作。前向操作是获取您的原始数据并重新投影它，但我们将使用所有的组件而不是较低的维度。您首先需要获得原始数据，但意味着减去：
Bm = bsxfun(@minus, B, mean(B,1));
Bm 将生成一个矩阵，其中每个样本的每个特征都被平均减去。 bsxfun 允许减去两个维度不等的矩阵，前提是您可以广播这些维度，以便它们都可以匹配。具体来说，在这种情况下，将计算B 的每一列/特征的平均值，并生成一个与B 一样大的临时复制矩阵。当您使用此复制矩阵减去原始数据时，效果将减去每个数据点及其各自的特征均值，从而分散您的数据，使每个特征的均值为 0。
一旦你这样做了，投影的操作很简单：
Bproject = Bm*Asort;
上面的操作很简单。您正在做的是将每个样本的特征表示为主成分的线性组合。例如，给定去中心化数据的第一个样本或第一行，投影域中的第一个样本的特征是属于整个样本的行向量和作为列向量的第一个主成分的点积。第一个样本在投影域中的第二个特征是整个样本和第二个分量的加权和。您将对所有样本和所有主成分重复此操作。实际上，您正在重新投影数据，使其与主成分相关 - 主成分是将数据从一种表示转换为另一种表示的正交基向量。
可以在此处找到对我刚才谈到的内容的更好描述。看看Amro的回答：
Matlab Principal Component Analysis (eigenvalues order)
现在倒退，您只需进行逆运算，但特征向量矩阵的一个特殊属性是，如果您将其转置，您将得到逆运算。要取回原始数据，您可以撤消上面的操作并将方法添加回问题：
out = bsxfun(@plus, Bproject*Asort.', mean(B, 1));
您想取回原始数据，因此您正在针对我之前执行的操作求解Bm。然而，Asort 的倒数只是这里的转置。执行此操作后发生的情况是您正在取回原始数据，但数据仍然是分散的。要取回原始数据，您必须将每个特征的均值添加回数据矩阵中以获得最终结果。这就是为什么我们在这里使用另一个 bsxfun 调用，以便您可以为每个样本的特征值执行此操作。
通过上面两行代码，你应该可以从原始域和投影域来回切换。现在降维（或原始数据的近似）发挥作用的是反向操作。您首先需要做的是将数据投影到主成分的基础上（即前向操作），但现在要回到我们试图用减少的主成分数量重建数据的原始域，您只需将上述代码中的Asort 替换为Aq，还可以减少您在Bproject 中使用的功能数量。具体来说：
out = bsxfun(@plus, Bproject(:,1:k)*Aq.', mean(B, 1));
执行Bproject(:,1:k) 会选择数据投影域中的k 特征，对应于k 最大特征向量。有趣的是，如果您只想要降低维度的数据表示，您可以使用Bproject(:,1:k)，这就足够了。但是，如果您想继续计算原始数据的近似值，我们需要计算反向步骤。上面的代码只是我们之前对你的数据的全维度所拥有的，但我们使用Aq 以及在Bproject 中选择k 功能。这将为您提供由矩阵中k 最大特征向量/特征值表示的原始数据。
如果你想看一个很棒的例子，我会模仿我链接到你的 Quora 帖子，但使用另一张图片。考虑使用灰度图像执行此操作，其中每一行是一个“样本”，每一列是一个特征。让我们以作为图像处理工具箱一部分的摄影师图像为例：
im = imread('camerman.tif');
imshow(im); %// Using the image processing toolbox
我们得到这张图片：
这是一个 256 x 256 的图像，这意味着我们有 256 个数据点，每个点有 256 个特征。我要做的是将图像转换为double 以精确计算协方差矩阵。现在我要做的是重复上述代码，但每次从 3、11、15、25、45、65 和 125 逐渐增加 k。因此，对于每个 k，我们将引入更多我们应该慢慢开始重建我们的数据。
这里有一些可运行的代码来说明我的观点：
%%%%%%%// Pre-processing stage
clear all;
close all;
%// Read in image - make sure we cast to double
B = double(imread('cameraman.tif'));
%// Calculate covariance matrix
sigma = cov(B);
%// Find eigenvalues and eigenvectors of the covariance matrix
[A,D] = eig(sigma);
vals = diag(D);
%// Sort their eigenvalues
[~,ind] = sort(abs(vals), 'descend');
%// Rearrange eigenvectors
Asort = A(:,ind);
%// Find mean subtracted data
Bm = bsxfun(@minus, B, mean(B,1));
%// Reproject data onto principal components
Bproject = Bm*Asort;
%%%%%%%// Begin reconstruction logic
figure;
counter = 1;
for k = [3 11 15 25 45 65 125 155]
%// Extract out highest k eigenvectors
Aq = Asort(:,1:k);
%// Project back onto original domain
out = bsxfun(@plus, Bproject(:,1:k)*Aq.', mean(B, 1));
%// Place projection onto right slot and show the image
subplot(4, 2, counter);
counter = counter + 1;
imshow(out,[]);
title(['k = ' num2str(k)]);
end
如您所见，大部分代码与我们看到的相同。不同的是我循环了k的所有值，用k最高特征向量投影回原始空间（即计算近似值），然后显示图像。
我们得到了这个漂亮的数字：
如您所见，以k=3 开头并没有真正给我们带来任何好处...我们可以看到一些一般结构，但添加更多内容并没有什么坏处。随着我们开始增加组件的数量，我们开始更清楚地了解原始数据的样子。在k=25，我们实际上可以完美地看到摄影师的样子，而且我们不需要 26 及以上的组件来查看正在发生的事情。这就是我所说的关于数据压缩的内容，您无需处理所有主要组件即可清楚地了解正在发生的事情。
我想通过向您推荐 Chris Taylor 关于主成分分析主题的精彩论述来结束这篇笔记，其中包含代码、图表和一个很好的解释！这是我开始使用 PCA 的地方，但 Quora 帖子巩固了我的知识。
Matlab - PCA analysis and reconstruction of multi dimensional data
【讨论】：
B的协方差矩阵。