tea*_*eef 4 matlab matrix correlation
我有一个大矩阵(1,000行和50,000列).我知道有些列是相关的(排名只有100),我怀疑有些列甚至是成比例的.我怎样才能找到这样的比例列？(一种方式是循环 corr(M(:,j),M(:,k))),但还有什么更有效的吗？
如果我正确理解您的问题,您希望确定矩阵中与线性相关的那些列,这意味着一列是比例或另一列的标量倍数.有一种基于QR分解的非常基本的算法.对于QR分解,您可以使用任何矩阵并将其分解为两个矩阵的乘积:Q和R.换一种说法:
A = Q*R
Q是一个正交矩阵,每列作为一个单位向量,这样乘以Q它的转置就可以得到单位矩阵(Q^{T}*Q = I). R是一个直角三角形或上三角形矩阵.Golub和Van Loan在他们1996年出版的一本书中非常有用的理论:矩阵计算是如果对角线元素的所有值都不为零,则矩阵被认为是满秩R.由于计算机上的浮点精度,我们必须阈值并检查对角线中任何R大于此容差的值.如果是,则该对应列是独立列.我们可以简单地找到所有对角线的绝对值,然后检查它们是否大于一些公差.
我们可以稍微修改它,以便我们搜索小于公差的值,这意味着列不是独立的.您调用QR分解的方式是这样的:
[Q,R] = qr(A, 0);
Q并且R是我刚才谈到的,你指定矩阵A作为输入.第二个参数0代表生成经济尺寸版本的Q和R,如果这个矩阵是矩形的(如你的情况),这将返回一个方形矩阵,其中尺寸是两个尺寸中最大的.换句话说,如果我有一个矩阵5 x 8,生成一个经济大小的矩阵将给你一个输出5 x 8,而不指定0将给你一个8 x 8矩阵.
现在,我们真正需要的是这种调用方式:
[Q,R,E] = qr(A, 0);
在这种情况下,E将是一个置换向量,这样:
A(:,E) = Q*R;
之所以这样,是非常有用的,因为它的命令列Q,并R以这样的方式的重新排序的版本的第一列是最有可能的列,它是独立的,其次是在递减的"实力"为了那些列.因此,E会告诉您每列的线性独立性有多大,并且这些"强度"按递减顺序排列.这种"强度"恰好在R对应于这种重新排序的对角线中捕获.实际上,强度与第一个元素成正比.您应该做的是检查R重新排列版本中的对角线是否大于由公差缩放的第一个系数,并使用这些来确定哪些相应列是线性独立的.
但是,我将翻转它并确定R最后可能的独立列所在的对角线中的点.在这一点之后的任何事情都将被视为线性依赖.这与检查是否有任何对角线小于阈值基本相同,但我们正在使用矩阵的重新排序对我们有利.
在任何情况下,将我在代码中提到的内容放在代码中,这就是你应该做的,假设你的矩阵存储在A:
%// Step #1 - Define tolerance
tol = 1e-10;
%// Step #2 - Do QR Factorization
[Q, R, E] = qr(A,0);
diag_R = abs(diag(R)); %// Extract diagonals of R
%// Step #3 -
%// Find the LAST column in the re-arranged result that
%// satisfies the linearly independent property
r = find(diag_R >= tol*diag_R(1), 1, 'last');
%// Step #4
%// Anything after r means that the columns are
%// linearly dependent, so let's output those columns to the
%// user
idx = sort(E(r+1:end));
请注意,这E将是一个置换向量,我假设您希望对它们进行排序,这就是为什么我们在向量无法变为线性独立的点之后对它们进行排序的原因.让我们测试一下这个理论.假设我有这个矩阵:
A =
1 1 2 0
2 2 4 9
3 3 6 7
4 4 8 3
您可以看到前两列是相同的,第三列是第一列或第二列的倍数.你只需要乘以2乘以得到结果.如果我们运行上面的代码,这就是我得到的:
idx =
1 2
如果你也看看E,这就是我得到的:
E =
4 3 2 1
这意味着第4列是"最佳"线性独立列,后面是第3列.因为我们返回[1,2]为线性相关列,这意味着列1和列2 [1,2,3,4]的列都是其他列的标量倍数.在这种情况下,这将是第3列,因为第1列和第2列是第3列的一半.
希望这可以帮助!
如果您不想进行任何QR分解,那么我可以建议将矩阵缩减为行减少的Echelon形式,并且可以确定构成矩阵列空间的基础向量A.实质上,如果要使用矩阵向量乘法应用此矩阵,则列空间为您提供可以生成输出向量的所有可能线性组合的最小列集.您可以使用该rref命令确定哪些列构成列空间.您将提供第二个输出rref,为您提供一个元素向量,告诉您哪些列是线性独立的,或者形成该矩阵的列空间的基础.因此:
[B,RB] = rref(A);
RB将为您提供列空间的列的位置,并且B将是矩阵的行减少的梯形形式A.因为要查找与线性相关的列,所以需要返回一组不包含这些位置的元素.因此,定义一个线性增加的向量,从1列到尽可能多的列,然后用于RB删除此向量中的这些条目,结果将是您正在寻找的线性相关列.换一种说法:
[B,RB] = rref(A);
idx = 1 : size(A,2);
idx(RB) = [];
通过使用上面的代码,这是我们得到的:
idx =
2 3
请记住,我们确定第2列和第3列是线性相关的,这是有意义的,因为它们都是第1列的倍数.与QR分解方法相比,哪些列线性相关的识别是不同的,因为QR基于列排序关于特定列线性独立的可能性.由于第1列到第3列彼此相关,因此返回哪一列无关紧要.其中一个形成了另外两个的基础.
rref与QR方法相比,我没有测试过使用效率.我怀疑是rref高斯行消除,与进行QR分解相比,复杂性更差,因为该算法是高度优化和高效的.因为你的矩阵相当大,我会坚持使用QR方法,但rref无论如何都要使用,看看你得到了什么!