从matlab中的一组向量中提取最大的一组线性无关向量
Extracting the largest set of linearly independent vectors from a set of vectors in matlab
我刚开始使用 matlab。我找不到一种明显的方法来从给定的向量集中提取线性无关向量的最大子集。
所以给定一个集合 V = [v1 v2 -- vn] 其中 dim(vi) >> n (i=1,2,3,....) 我需要线性地找到一组随机的 r独立向量 "vi"（其中 r 最大），即从 V 中删除所有 (n-r) 线性相关 "vi" 向量。很抱歉，如果这不明确但找不到更好的方式来说明它。
我会按照这些思路尝试一些东西。不过，它可能不会很快：
- 计算由您的向量构成的矩阵的秩。调用 r.
- 生成取自
1:n. 的 r 个数字的所有组合
- 对这些组合进行随机排序。
- 取第一个组合并检查该组合指示的向量是否独立（最大秩）。如果不是，请测试下一个组合。通过构造，可以保证某些组合会成功。
代码：
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
我假设您的向量都是 n 维的，并且我们可以将它们全部连接成一个矩阵。如果我们进入矩阵和线性代数，您正在寻找的是矩阵的 the Column Space 。简单地说，列 space 被定义为矩阵中的一组列，可以在 n 维 space 中唯一地生成另一个向量。或者，它是列向量所有可能的线性组合的集合。因此，如果你想找到最大的线性无关向量集，你所要做的就是确定矩阵的第 space 列。
因此，给定大小为 n x m 的矩阵 V，我们有 m 列/向量，每列的大小为 n x 1（或n 行），您将调用 rref or the Row-Reduced Echelon Form (RREF) command. This reduces your matrix down to its row-reduced echelon form。这是查找矩阵的 space 列的开始。你可以这样称呼它：
[R, RB] = rref(V);
R 将包含 V 的 RREF 形式，RB 将包含 R 的索引或列号， 构成列 space。因此，如果你想生成线性无关向量，你只需要做：
VMax = V(:,RB);
VMax 将仅包含构成列 space 的 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 告诉我们哪些 列 是线性独立的或构成 space 的列A。如您所见，只有第 1 列和第 4 列是线性独立的，这非常合理。如果您还看一下 R 的最后两行，我们可以看到它们由全零组成。这暗指矩阵的 rank，其中它只是非零的行（或列，具体取决于您正在做的事情）的总数。这也告诉您有多少向量构成列 space.
要完成您的任务，只需执行以下操作：
>> VMax = V(:,RB)
VMax =
1 0
2 9
3 7
4 3
如你所见，VMax的每一列都是来自V的线性独立向量，它也构成了V的第space列。
现在，您的下一个任务是每次 运行 算法时，从该列 space 中随机选择线性无关向量。请记住，因为有多种方法可以生成列 space，rref 的解决方案只会给你一个这样的列 space。如果我正确地解释了你的问题，你希望生成随机列 spaces 并每次选择这些向量的一个子集。感谢 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，这样您就可以迫使它选择不同的基向量，但是如果您有多个线性相关的向量，那么在某些情况下我们会选择一个这些线性相关的向量来构建我们的基础。
一个例子的部分解决方案......抱歉，我没有注意到你的评论
“different 线性独立集”......你是说"different maximal linearly independent sets"？抱歉，我只找到了最大线性独立集之一......但是使用 [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
The left-singular vectors corresponding to the non-zero singular values of M span the range of M.
这种寻找向量所跨越的向量 space 的基础问题的天真方法类似于 Luis Mendo 的。然而，与 n choose r 相比，它应该具有更好的理论最坏情况复杂度，最多 n rank 次计算。缺点是，最好的情况将执行 r rank 计算，而不是 1.
想法很简单：
- 计算
V 的 rank r。这是您的集合中的最大向量数。
- 将
chosen 向量集定义为空。
- 对于
V中的每个向量V(:,k)（随机选择）：
- 如果向量
V(:,k) 与集合 chosen 中的向量线性无关：
- 将其添加到所选向量集中。
- 如果集合
chosen 已经有 r 个元素：
- Return
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 的 rref 解决方案更优雅，但出于某种原因，即使对于非稀疏矩阵，这种朴素的方法似乎也更快。
加速版本使用 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
在数学上你会说这个矩阵有五阶。然而，与其他列相比，在数值上只有最后两列是非零的。这就是 rank(V) 给出 2 的原因。这种迭代方法的问题在于，向我们的集合中添加一个线性无关的向量实际上会导致 线性相关的集合 ！
想象一下这个算法选择第一个向量。尽管所有其他向量在数学上和数值上都是线性独立的，但子集 {V(:,1), V(:,5)} 没有数值等级 2，但可以由我们的函数选择。这就是为什么该算法只适用于相似范数的向量。
结论是：
- 数值数学很棘手。
- 想要稳定，就得付出计算成本的代价
rref 是解决这个问题的方法。
我刚开始使用 matlab。我找不到一种明显的方法来从给定的向量集中提取线性无关向量的最大子集。
所以给定一个集合 V = [v1 v2 -- vn] 其中 dim(vi) >> n (i=1,2,3,....) 我需要线性地找到一组随机的 r独立向量 "vi"（其中 r 最大），即从 V 中删除所有 (n-r) 线性相关 "vi" 向量。很抱歉，如果这不明确但找不到更好的方式来说明它。
我会按照这些思路尝试一些东西。不过，它可能不会很快：
- 计算由您的向量构成的矩阵的秩。调用 r.
- 生成取自
1:n. 的 - 对这些组合进行随机排序。
- 取第一个组合并检查该组合指示的向量是否独立（最大秩）。如果不是，请测试下一个组合。通过构造，可以保证某些组合会成功。
r 个数字的所有组合
代码：
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
我假设您的向量都是 n 维的，并且我们可以将它们全部连接成一个矩阵。如果我们进入矩阵和线性代数，您正在寻找的是矩阵的 the Column Space 。简单地说，列 space 被定义为矩阵中的一组列，可以在 n 维 space 中唯一地生成另一个向量。或者，它是列向量所有可能的线性组合的集合。因此，如果你想找到最大的线性无关向量集，你所要做的就是确定矩阵的第 space 列。
因此，给定大小为 n x m 的矩阵 V，我们有 m 列/向量，每列的大小为 n x 1（或n 行），您将调用 rref or the Row-Reduced Echelon Form (RREF) command. This reduces your matrix down to its row-reduced echelon form。这是查找矩阵的 space 列的开始。你可以这样称呼它：
[R, RB] = rref(V);
R 将包含 V 的 RREF 形式，RB 将包含 R 的索引或列号， 构成列 space。因此，如果你想生成线性无关向量，你只需要做：
VMax = V(:,RB);
VMax 将仅包含构成列 space 的 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 告诉我们哪些 列 是线性独立的或构成 space 的列A。如您所见，只有第 1 列和第 4 列是线性独立的，这非常合理。如果您还看一下 R 的最后两行，我们可以看到它们由全零组成。这暗指矩阵的 rank，其中它只是非零的行（或列，具体取决于您正在做的事情）的总数。这也告诉您有多少向量构成列 space.
要完成您的任务，只需执行以下操作：
>> VMax = V(:,RB)
VMax =
1 0
2 9
3 7
4 3
如你所见，VMax的每一列都是来自V的线性独立向量，它也构成了V的第space列。
现在，您的下一个任务是每次 运行 算法时，从该列 space 中随机选择线性无关向量。请记住，因为有多种方法可以生成列 space，rref 的解决方案只会给你一个这样的列 space。如果我正确地解释了你的问题，你希望生成随机列 spaces 并每次选择这些向量的一个子集。感谢 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，这样您就可以迫使它选择不同的基向量，但是如果您有多个线性相关的向量，那么在某些情况下我们会选择一个这些线性相关的向量来构建我们的基础。
一个例子的部分解决方案......抱歉，我没有注意到你的评论
“different 线性独立集”......你是说"different maximal linearly independent sets"？抱歉，我只找到了最大线性独立集之一......但是使用 [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
The left-singular vectors corresponding to the non-zero singular values of M span the range of M.
这种寻找向量所跨越的向量 space 的基础问题的天真方法类似于 Luis Mendo 的。然而，与 n choose r 相比，它应该具有更好的理论最坏情况复杂度，最多 n rank 次计算。缺点是，最好的情况将执行 r rank 计算，而不是 1.
想法很简单：
- 计算
V的rankr。这是您的集合中的最大向量数。 - 将
chosen向量集定义为空。 - 对于
V中的每个向量V(:,k)（随机选择）：- 如果向量
V(:,k)与集合chosen中的向量线性无关：- 将其添加到所选向量集中。
- 如果集合
chosen已经有r个元素：- Return
chosen个向量作为解决方案。
- Return
- 如果向量
生成的代码如下所示：
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 的 rref 解决方案更优雅，但出于某种原因，即使对于非稀疏矩阵，这种朴素的方法似乎也更快。
加速版本使用 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
在数学上你会说这个矩阵有五阶。然而，与其他列相比，在数值上只有最后两列是非零的。这就是 rank(V) 给出 2 的原因。这种迭代方法的问题在于，向我们的集合中添加一个线性无关的向量实际上会导致 线性相关的集合 ！
想象一下这个算法选择第一个向量。尽管所有其他向量在数学上和数值上都是线性独立的，但子集 {V(:,1), V(:,5)} 没有数值等级 2，但可以由我们的函数选择。这就是为什么该算法只适用于相似范数的向量。
结论是：
- 数值数学很棘手。
- 想要稳定，就得付出计算成本的代价
rref是解决这个问题的方法。