matlab - 在协方差矩阵中选择最大特征值和特征向量在数据分析中意味着什么？

Question

假设存在一个矩阵B，其大小为500*1000 double（此处，500表示观测值的数量，1000表示要素的数量）。

sigma是B的协方差矩阵，D是对角线矩阵，其对角线元素是sigma的特征值。假设A是协方差矩阵sigma的特征向量。

我有以下问题：


我需要选择与最大特征值相对应的第一个k = 800特征向量，以对所选特征进行排名。最后的矩阵名为Aq。如何在MATLAB中执行此操作？
这些选定的特征向量是什么意思？
一旦计算出Aq，看来最终矩阵1000*800 double的大小为Aq。 500的时间点/观察信息已消失。对于最终矩阵Aq，矩阵1000中的值Aq现在代表什么？另外，矩阵800中的值Aq现在代表什么？

Answer 1

我假设您从eig函数确定了特征向量。我以后建议您使用eigs函数。这不仅为您计算特征值和特征向量，还将为您计算k最大特征值及其关联的特征向量。这样就可以节省计算开销，因为您只需要一个子集，就不必计算矩阵的所有特征值和相关特征向量。您只需将数据的协方差矩阵提供给eigs，它就会为您返回k最大的特征值和特征向量。



现在，回到您的问题，您所描述的最终是Principal Component Analysis。其背后的机制是计算数据的协方差矩阵，并找到计算结果的特征值和特征向量。众所周知，由于计算大型矩阵的特征值和特征向量时数值不稳定，因此不建议这样做。现在最典型的方法是通过Singular Value Decomposition。具体而言，V矩阵的列为您提供协方差矩阵的特征向量或主成分，并且相关的特征值是矩阵S的对角线中产生的奇异值的平方根。

请参阅有关交叉验证的信息，以了解为什么偏爱此方法：

https://stats.stackexchange.com/questions/79043/why-pca-of-data-by-means-of-svd-of-the-data

我还将在另一个链接中讨论在主成分分析中为何使用奇异值分解的理论：

https://stats.stackexchange.com/questions/134282/relationship-between-svd-and-pca-how-to-use-svd-to-perform-pca



现在让我们一次回答您的问题。

问题1

MATLAB以未排序的方式生成特征值和特征向量的相应顺序。如果要根据给定的k输出（在本例中为800）选择最大的eig特征值和相关特征向量，则需要按降序对特征值进行排序，然后重新排列产生的特征向量矩阵的列从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的第i列/特征向量对应于 D中的第i个特征值。

但是，特征值在对角矩阵中，因此我们使用 diag命令提取对角线，对其进行排序并找出它们的顺序，然后重新排列 A以遵守该顺序。我使用 sort的第二个输出，因为它告诉您未排序结果中每个值将出现在排序结果中的位置。这是我们需要重新排列特征向量矩阵 A的列的顺序。必须选择 'descend'作为标志，以便最大的特征值和相关的特征向量首先出现，就像我们之前讨论过的一样。

然后，您可以通过以下方式提取前 k个最大的向量和值：

k = 800;
Aq = Asort(:,1:k);

问题2

众所周知，协方差矩阵的特征向量等于主成分。具体而言，第一个主成分（即最大特征向量和关联的最大特征值）为您提供数据中最大可变性的方向。此后的每个主要成分都会为您带来递减的可变性。还应注意，每个主成分彼此正交。

这是来自Wikipedia的一个很好的二维数据示例：



我从Wikipedia上关于主成分分析的文章中提取了上面的图片，并将其链接到上面。这是根据以 (1,3)为中心的双变量高斯分布进行分布的样本的散点图，其标准偏差在 (0.878, 0.478)方向上为3，在正交方向上为1。标准偏差为3的分量是第一主分量，而正交的分量是第二分量。所示向量是协方差矩阵的特征向量，该特征向量由相应特征值的平方根缩放，并进行了平移，因此它们的尾部位于平均位置。

现在，让我们回到您的问题。我们看一下 k最大特征值的原因是执行 dimensionality reduction的一种方式。本质上，您将执行数据压缩，其中将获取较高维度的数据并将其投影到较低维度的空间。投影中包含的主要成分越多，它与原始数据的相似度就越高。实际上，它在某个点开始逐渐减小，但是前几个主要组件使您可以在大多数情况下忠实地重建数据。

我过去偶然发现的这篇出色的Quora文章找到了执行PCA（或SVD）和数据重建的绝佳视觉示例。

http://qr.ae/RAEU8a

问题＃3

您将使用此矩阵将较高维度的数据重新投影到较低维度的空间。行数仍为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并在 k中选择 Bproject功能。这将为您提供由矩阵中 k个最大特征向量/特征值表示的原始数据。



如果您想看到一个很棒的示例，我将模仿我链接到您的Quora帖子，但使用另一个图像。考虑使用灰度图像执行此操作，其中每一行都是一个“样本”，每一列都是一个要素。让我们以摄影师图像作为图像处理工具箱的一部分：

im = imread('camerman.tif');
imshow(im); %// Using the image processing toolbox

我们得到这个图像：



这是一个256 x 256的图像，这意味着我们有256个数据点，每个点都有256个特征。我要做的是将图像转换为 double，以精确计算协方差矩阵。现在，我要做的是重复上面的代码，但是每次将 k分别从3、11、15、25、45、65和125递增。因此，对于每个 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