在numpy或matlab中从一个满秩的非方阵矩阵中获取可逆方阵。

Question

在numpy或matlab中从一个满秩的非方阵矩阵中获取可逆方阵。

4

假设您有一个满秩的NxM矩阵A，其中M>N。如果我们用Ci（维度为Nx1）表示列，则可以将矩阵写成：

A = [C_1, C_2, ..., C_M]

如何获取原始矩阵 A 的第一个线性独立列，以便您可以构造一个新的可逆矩阵 B，其行列式不为零。

B = [C_i1, C_i2, ..., C_iN]

你如何在matlab或python numpy中找到指数{i1，i2，...，iN}？这可以使用奇异值分解完成吗？非常欢迎提供代码片段。

编辑：为了更具体地说明，考虑以下python代码

from numpy import *
from numpy.linalg.linalg import det

M = [[3, 0, 0, 0, 0],
     [0, 0, 1, 0, 0],
     [0, 0, 0, 0, 1], 
     [0, 2, 0, 0, 0]]
M = array(M)

I = [0,1,2,4]
assert(abs(det(M[:,I])) > 1e-8)

因此，对于给定的矩阵M，需要找到一组线性无关列向量的索引，其中列向量的数量为N。

- D R

2个回答

1

我的第一个想法是尝试从M列中选择N列的每种可能组合。可以使用以下Python代码实现：

import itertools
import numpy.linalg

# 'singular' returns whether a matrix is singular.
# You could use something more efficient than the determinant
# (I'm not sure what options there are in NumPy)
singular = lambda m: numpy.linalg.det(m) == 0

def independent_square(A):
    N,M = A.shape
    for colset in itertools.combinations(xrange(M), N):
        B = A[:,colset]
        if not singular(B):
            return B

如果你想要列索引而不是结果矩阵，只需将return B替换为return colset。或者你可以使用return colset,B同时得到两者。

我不知道SVD在这里有任何帮助的方式。事实上，我无法想出任何纯数学运算可以将A转换为B（甚至找出MxN列选择矩阵Q，使得B=A.Q），除了通过试错。但如果你想找出是否存在这样的方法，math.stackexchange.com是一个很好的问答平台。

如果你只需要计算方法，上述代码应该足够了。

- David Z

谁给踩了，请说明为什么你认为这个过程是错误的。 - David Z

我没有给它点踩，但是它真的很贵。行列式计算很昂贵：如果numpy提供了行约简，那么就会好得多（我不知道）。一种更好的算法是归纳的：假设n-1列是线性独立的，然后继续尝试下一个未使用的列，直到找到另一个线性独立的列为止。重复这个过程。请记住，即使您可能不必搜索其中的所有组合（除非您这样做），但您愿意搜索的组合有M!/(N!(M - N)!) 种。 - aaronasterling

我不知道Numpy是否也提供行降阶... 行列式只是我能想到的一个确定列的线性独立性的方法，因为我找不到一个rank函数（好吧，确实有一个，但它并不像你想象的那样工作）。我想我应该编辑一下，说任何用于确定列是否线性独立的过程都可以在那里发挥作用。 - David Z

1

你想知道为什么这是一个糟糕的解决方案吗？这个解决方案的问题在于，它可能非常低效，具体取决于M和N的值。例如，当N = 100且M = 50时，有超过1e29种可能的列组合。通过列的暴力搜索就太疯狂了。考虑到QR分解中有一个非常好的解决方案，使用像暴力搜索这样的低效算法是愚蠢的。更糟糕的是，我认为你最初的回答建议使用行列式来确定奇异性。对任何矩阵使用det都是一件可怕的事情。 - user85109

我想知道为什么它是“错误”的，而不是“差劲的”。此外，在Python中QR方法不起作用，因为NumPy/SciPy的qr方法没有从您的解决方案中提供E矩阵。 - David Z

1

没有人说这在数学上是错误的。我建议不要提供这种解决问题的可怕建议，并给出了一个很好的理由。如果Python中没有旋转QR，那么仍然有其他选择。你可以建议一种Gram-Schmidt的变体，它与QR并没有太大区别。因此，编写一个基于Gram-Schmidt的旋转形式，它仅仅是减去指向到目前为止找到的前面列方向的分量，然后在剩余幅度最大的列上进行旋转。（这将非常容易编写。） - user85109

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- user85109 · Accepted Answer

在MATLAB中，使用QR分解来解决问题非常容易。特别是，使用带有枢轴的QR分解会更加方便。

M = [3 0 0 0 0;
     0 0 1 0 0;
     0 0 0 0 1; 
     0 2 0 0 0]

[Q,R,E] = qr(M)
Q =
     1     0     0     0
     0     0     1     0
     0     0     0     1
     0     1     0     0

R =
     3     0     0     0     0
     0     2     0     0     0
     0     0     1     0     0
     0     0     0     1     0

E =
     1     0     0     0     0
     0     1     0     0     0
     0     0     1     0     0
     0     0     0     0     1
     0     0     0     1     0

E的前4列指定要使用M的哪些列，即[1,2,3,5]列。如果你想要M的列，只需将M与E相乘即可。

M*E
ans =
     3     0     0     0     0
     0     0     1     0     0
     0     0     0     1     0
     0     2     0     0     0

顺便说一下，使用行列式来确定矩阵是否奇异绝对是你能做到的最糟糕的方式。

请改用秩进行判断。

基本上，在MATLAB中，除非你了解为什么使用行列式是这样一个坏事，并且你选择尽管如此使用它，否则你几乎永远不应该使用行列式。