在矩阵中查找比例列

Question

我有一个大矩阵（1,000 行和 50,000 列）。我知道有些列是相关的（排名只有 100），我怀疑有些列甚至是成比例的。我怎样才能找到这样的比例列？ （一种方法是循环corr(M(:,j),M(:,k))），但有什么更有效的方法吗？

Answer 1

如果矩阵的行列式为零，则列是成比例的。

有 50,000 列，或 2 乘以 25,000 列。

最容易求解 2×2 矩阵的行列式。

因此：要找到比例矩阵，最长时间的解决方案是 在电子表格上定义大矩阵。

将行列式公式应用于从左侧第一个正方形开始的正方形。 为每一行和每一列复制它，以便在下一个电子表格中得到答案。

找出行列式为零的列。

这是相当基本的，不是很耗时，应该以结果为导向。 手动或 Excel 电子表格（高效）

Answer 2

如果您通过除以最大值来标准化每列，则比例变为相等。这使问题变得更容易。

现在，要测试是否相等，您可以在列上使用单个（外部）循环；内部循环很容易用bsxfun 向量化。为提高速度，请仅将每列与其右侧的列进行比较。

另外为了节省一些时间，结果矩阵被预先分配到一个近似的大小（你应该设置它）。如果近似大小有误，唯一的惩罚就是速度会慢一点，但代码可以工作。

像往常一样，浮点值之间的相等性测试应该包括一个容差。

结果以 2 列矩阵 (S) 的形式给出，其中每一行包含成比例的两行的索引。

A = [1 5 2 6 3 1
     2 5 4 7 6 1
     3 5 6 8 9 1]; %// example data matrix
tol = 1e-6; %// relative tolerance

A = bsxfun(@rdivide, A, max(A,[],1)); %// normalize A
C = size(A,2);
S = NaN(round(C^1.5),2); %// preallocate result to *approximate* size
used = 0; %// number of rows of S already used
for c = 1:C
    ind = c+find(all(abs(bsxfun(@rdivide, A(:,c), A(:,c+1:end))-1)<tol));
    u = numel(ind); %// number of columns proportional to column c
    S(used+1:used+u,1) = c; %// fill in result
    S(used+1:used+u,2) = ind; %// fill in result
    used = used + u; %// update number of results
end
S = S(1:used,:); %// remove unused rows of S

在这个例子中，结果是

S =
     1     3
     1     5
     2     6
     3     5

意思是第1列与第3列成比例；第 1 列与第 5 列等成比例。

Answer 3

如果我正确理解您的问题，您希望确定矩阵中的那些列线性相关，这意味着一列与另一列成比例或标量倍数。有一个基于QR Decomposition 的非常基本的算法。对于 QR 分解，您可以将任何矩阵分解为两个矩阵的乘积：Q 和 R。换句话说：

A = Q*R

Q 是一个正交矩阵，每一列都是一个单位向量，因此将Q 与其转置相乘即可得到单位矩阵 (Q^{T}*Q = I)。 R 是一个直角三角形或上三角矩阵。 Golub and Van Loan in their 1996 book: Matrix Computations 提出的一个非常有用的理论是，如果 R 的对角线元素的所有值都不为零，则矩阵被认为是满秩。由于计算机上的浮点精度，我们必须设置阈值并检查R 的对角线中大于此容差的任何值。如果是，那么这个对应的列就是一个独立的列。我们可以简单地找到所有对角线的绝对值，然后检查它们是否大于某个容差。

我们可以稍微修改一下，以便我们搜索小于小于的值，这意味着该列不独立。您调用 QR 分解的方式是这样的：

[Q,R] = qr(A, 0);

Q 和R 就是我刚才讲的，你指定矩阵A 作为输入。第二个参数0 代表生成Q 和R 的economy-size 版本，如果这个矩阵是矩形的（就像你的情况一样），这将返回一个方阵，其中尺寸是两种尺寸中最大的。换句话说，如果我有一个像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 分解，那么我可以建议将您的矩阵缩减为它的row-reduced Echelon form，并且您可以确定构成矩阵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 列相互关联，所以返回哪一列并不重要。其中一个是其他两个的基础。

与 QR 方法相比，我还没有测试过使用 rref 的效率。我怀疑rref 进行了高斯行消除，与进行 QR 分解相比，复杂性更差，因为该算法高度优化且高效。因为你的矩阵相当大，我会坚持使用 QR 方法，但还是使用rref 看看你得到了什么！