从matlab中的一组向量中提取最大的线性独立向量集

Question

我是使用 matlab 的新手。我找不到一种明显的方法来从给定的向量集中提取线性独立向量的最大子集。

所以给定一个集合 V = [v1 v2 -- vn] where dim(vi) >> n (i=1,2,3,....) 我需要找到一组随机的 r 个线性独立向量“vi”（其中 r 最大），即从 V 中删除所有 (n-r) 线性相关的“vi”向量。抱歉，如果这是模棱两可但找不到更好的方式来说明它。

Answer 1

我将假设您的向量都是n 维，并且我们可以将它们全部连接到一个矩阵中。如果我们进入矩阵和线性代数，您正在寻找的是矩阵的the Column Space。简单地说，列空间被定义为矩阵中的一组列，它们可以在n维空间中唯一地产生另一个向量。或者，它是列向量的所有可能线性组合的集合。因此，如果您想找到最大的线性独立向量集，您所要做的就是确定矩阵的列空间是什么。

因此，给定您的矩阵V，其大小为n x m，其中我们有m 列/向量，每列的大小为n x 1（或n 行），您可以调用rref 或 Row-Reded Echelon Form (RREF) 命令。这会将您的矩阵减少到其row-reduced echelon form。这是为矩阵寻找列空间的开始。你可以这样称呼它：

[R, RB] = rref(V);

R 将包含 V 的 RREF 形式，RB 将包含 构成列空间的 R 的索引或列号。因此，如果你想产生你的线性独立向量，你只需要这样做：

VMax = V(:,RB);

VMax 将仅包含构成列空间的V 的那些列，因此也包含那些是线性独立向量的列。如果您想确定我们有 多少 个独立向量，您只需计算我们有多少 RB 的值：

r = numel(RB);

---

这里有一个简单的例子来说明我的观点。假设我有这个矩阵：

>> V = [1 1 2 0; 2 2 4 9; 3 3 6 7; 4 4 8 3]

V =

     1     1     2     0
     2     2     4     9
     3     3     6     7
     4     4     8     3

第二列/向量只是第一个向量。第三列/向量只是第一个向量加上第二个向量，或者它可以是第一个或第二个向量的两倍。无论哪种方式，这不是线性独立向量，因为它基于第一个向量。第二个向量也是如此。最后一个向量独立于其他三个向量，因为我们无法创建组合或缩放来从其他三个向量生成最后一个向量。如果我们调用rref，将会发生这种情况：

>> [R, RB] = rref(V)

R =

     1     1     2     0
     0     0     0     1
     0     0     0     0
     0     0     0     0


RB =

     1     4

R 包含行缩减梯队形式，而RB 告诉我们哪些列 是线性独立的或构成A 的列空间。如您所见，只有第 1 列和第 4 列是线性独立的，这很有意义。如果您还查看R 的最后两行，我们可以看到它们由全零组成。这暗示了矩阵的rank，它只是非零的总行数（或列数，取决于你在做什么）。这也告诉你有多少向量构成列空间。

要完成您的任务，只需执行以下操作：

>> VMax = V(:,RB)

VMax =

     1     0
     2     9
     3     7
     4     3

如您所见，VMax 的每一列都是来自V 的线性独立向量，这也构成了V 的列空间。

现在，您的下一个任务是每次运行算法时从该列空间中随机选择线性独立向量。请记住，由于有多种方法可以产生列空间，rref 的解决方案只会给你一个这样的列空间。如果我正确解释了您的问题，您希望生成随机列空间并每次选择这些向量的子集。感谢 Luis Mendo（也是来自 knedlsepp 的温和刺激），您可以做的是随机重新排列或置换列，然后在这个置换矩阵上运行 rref。

因此，您可以这样做：

ind = randperm(size(V,2));
Vperm = V(:,ind);
[R,RB] = rref(Vperm);
VMax = Vperm(:,RB);

randperm 将生成一个从 1 到您指定给 randperm 的参数的随机排列数组。在我们的例子中，这是矩阵V 中的总列数。我们使用这个数组随机打乱V 的列，将其存储到Vperm 并运行我们之前完成的代码。通过进行这种随机洗牌，您将输入偏置到rref，以便您强制它选择不同的基向量，但是如果您有多个线性相关的向量，那么在某些情况下我们将选择其中一个线性相关向量来构建我们的基础。

Answer 2

这种寻找由向量跨越的向量空间的基础的简单方法类似于 Luis Mendo 的。然而，与n choose r 相比，它应该具有更好的理论最坏情况复杂度，最大为nrank-computations。缺点是，最好的情况将执行rrank-computations，而不是1。

想法很简单：

• 计算V 的rank r。这是您的集合中向量的最大数量。
• 将chosen 向量集定义为空。
• 对于V 中的每个向量V(:,k)（随机选择）：
  • 如果向量V(:,k)与集合chosen中的向量线性无关：
    • 将其添加到所选向量集。
  • 如果集合chosen 已经有r 元素：
    • 返回 chosen 向量作为解。

生成的代码将如下所示：

function chosen = randBasis(V)
% V is a matrix of column vectors
% chosen are the indices of the randomly selected column vectors,
% so that span(V(:,chosen)) == span(V)

n = size(V,2);
r = rank(V);
chosen = []; 
currentRank = 0;
for k = randperm(n)
    if rank(V(:, [chosen, k])) > currentRank
        currentRank = currentRank+1;
        chosen = [chosen, k];
    end
    if currentRank==r
        break
    end
end

补充说明：

当您处理稀疏矩阵时，这种方法可能很有用，因为检查线性相关性可以用线性系统的解决方案代替秩计算。由于排名计算涉及对svd 的调用，因此相比之下它相当昂贵。尽管 rayryeng​​em> 的 rref-solution 是更优雅的解决方案，但出于某种原因，即使对于非稀疏矩阵，这种幼稚的方法似乎也更快。

使用mldivide 加速版本：

function chosen = randBasis(V)
n = size(V,2);
if issparse(V)
    r = sprank(V)
else
    r = rank(V);
end
chosen = []; 
for k = randperm(n)
    if isLinearlyIndependent(V(:,chosen),V(:,k))
        chosen = [chosen, k];
    end
    if numel(chosen)==r
        break
    end
end
end

function b = isLinearlyIndependent(V, v)
    warning('off','MATLAB:singularMatrix');
    bestapprox = V\v;
    bestapprox(~isfinite(bestapprox)) = 0;
    b = ~(norm(V*bestapprox-v,'inf')<=eps*norm(V,'inf'));
    warning('on','MATLAB:singularMatrix');
end

注意：

这两种迭代方法都存在一个不明显的数值问题： 考虑矩阵：

V = diag(10.^(-40:10:0))
V =
        1e-40            0            0            0            0
            0        1e-30            0            0            0
            0            0        1e-20            0            0
            0            0            0        1e-10            0
            0            0            0            0        1e+00

从数学上讲，您会说这个矩阵的秩为 5。然而，与其他列相比，只有最后两列在数值上是非零的。这就是rank(V) 给出2 的原因。这种迭代方法的问题在于，将线性独立向量添加到我们的集合实际上会导致线性相关集合！ 想象一下这个算法选择了第一个向量。尽管所有其他向量在数学和数值上都是线性独立的，但子集{V(:,1), V(:,5)} 没有数字等级 2，但可以由我们的函数选择。这就是为什么该算法仅适用于具有相似范数的向量的原因。

结论是：

• 数值数学很棘手。
• 想要稳定，就必须付出计算成本的代价。
• rref 是解决这个问题的方法。

Answer 3

我会尝试这些方面的东西。不过，它可能不会很快：

1. 计算向量形成的矩阵的秩。称它为 r。
2. 生成取自1:n 的r 数字的所有组合。
3. 对这些组合进行随机排序。
4. 采用第一个组合并检查该组合指示的向量是否独立（最大等级）。如果不是，请测试下一个组合。通过构造，可以保证某种组合会成功。

代码：

vectors = randi(3,4,3); %// example. Each vector is a row
n = size(vectors,1);
r = rank(vectors); %// step 1
combs = nchoosek(1:n,r); %// step 2
c = size(combs,1);
combs = combs(randperm(c),:); %// step 3
for s = 1:r %// step 4
    pick = combs(s,:);
    if rank(vectors(pick,:))==r
        break %// we're done. Now `pick` contains the selected vectors
    end
end
result = vectors(pick,:); %// each row is a vector

Answer 4

一个例子的部分解决方案......对不起，我没有注意到你的评论 “不同线性独立集”...您的意思是“不同最大线性独立集”吗？抱歉，我只找到了最大线性独立集之一......但是使用[R, RB] = rref(M) 的解决方案给出了相同的答案

+ echo ('on')
+ more ('off')
+ format ('bank')
+ M = [0, 0, 1; 0, 1, 1; 0, 1, 0]'
M =

   0.00   0.00   0.00
   0.00   1.00   1.00
   1.00   1.00   0.00

+ [U, S, V] = svd (M)
U =

   0.00   0.00  -1.00
  -0.71   0.71   0.00
  -0.71  -0.71   0.00

S =

Diagonal Matrix

   1.73      0      0
      0   1.00      0
      0      0  -0.00

V =

  -0.41  -0.71  -0.58
  -0.82   0.00   0.58
  -0.41   0.71  -0.58

+ U (:, 1:2) * S (1:2, 1:2) * V (1:2, 1:2)'
ans =

   0.00   0.00
   0.00   1.00
   1.00   1.00

来自Applications_of_the_SVD的引述

M的非零奇异值对应的左奇异向量跨越M的范围。