在协方差矩阵中选择最大特征值和特征向量的意义在于数据分析？

Question

假设有一个矩阵B,其大小为a 500*1000 double(这里,500代表观察1000的数量并代表特征的数量).

sigma是协方差矩阵B,D是一个对角矩阵,其对角元素是特征值sigma.假设A是协方差矩阵的特征向量sigma.

我有以下问题:

1. 我需要选择k = 800对应于具有最大幅度的特征值的第一个特征向量来对所选特征进行排序.最终的矩阵命名Aq.我怎样才能在MATLAB中做到这一点？

2. 这些选定的特征向量是什么意思？

3. 看来最终的矩阵的大小Aq是1000*800 double有一次我计算Aq.时间点/观察信息500已经消失.对于最终的矩阵Aq,什么是价值1000矩阵Aq现在代表什么呢？此外,800矩阵中的值Aq现在代表什么？

Answer 1

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

现在,回到你的问题,你所描述的最终是主成分分析.这背后的机制是计算数据的协方差矩阵,并找到计算结果的特征值和特征向量.众所周知,由于计算大矩阵的特征值和特征向量的数值不稳定性,不推荐这样做.现在最常规的方法是通过奇异值分解.具体地说,V矩阵的列给出协方差矩阵或主成分的特征向量,相关的特征值是矩阵对角线中产生的奇异值的平方根S.

请参阅有关Cross验证的信息,以了解为何首选:

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特征值和给定的输出相关的特征向量eig(800在您的示例),你需要排序的特征值按降序排列,然后重新从产生的特征向量矩阵的列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代表一个特征向量.具体而言,第i ^个列/特征向量的A对应于I ^个本征值中所看到的D.

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

然后,您可以k通过以下方式获取第一个最大的向量和值:

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

问题2

众所周知,协方差矩阵的特征向量等于主成分.具体而言,第一个主成分(即最大的特征向量和相关的最大特征值)为您提供数据中最大可变性的方向.之后的每个主要组成部分都为您提供了降低性质的可变性.值得注意的是,每个主成分彼此正交.

以下是来自维基百科的二维数据的一个很好的例子:

我从维基百科关于主成分分析的文章中提取了上面的图像,我将其链接到上面.这是样本的散点图,其根据以(1,3)大致(0.878, 0.478)方向为标准偏差3 且在正交方向上为1 的标准偏差为中心的双变量高斯分布来分布.标准偏差为3的组件是第一主成分,而正交的组件是第二主成分.所示的矢量是由相应特征值的平方根缩放的协方差矩阵的特征向量,并且移位使得它们的尾部处于平均值.

现在让我们回到你的问题.我们看一下k最大特征值的原因是一种降低维数的方法.从本质上讲,您将执行数据压缩,您可以在其中获取更高维度的数据并将其投影到较低维度的空间.投影中包含的主要组件越多,它就越像原始数据.它实际上开始逐渐减少,但前几个主要组件允许您在大多数情况下忠实地重建数据.

在我过去偶然发现的这个伟大的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主成分分析(特征值顺序)

现在要倒退,你只需要进行逆运算,但是具有特征向量矩阵的特殊属性是,如果你转置它,你得到逆.要恢复原始数据,请撤消上面的操作并将问题添加回问题:

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分析和多维数据重建