从numpy或matlab中的满秩非方阵获取可逆方阵

Question

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

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

如何获得原始矩阵A的第一个线性独立列，以便可以构造一个新的NxN矩阵B，它是一个具有非零行列式的可逆矩阵。

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

如何在 matlab 或 python numpy 中找到索引 {i1, i2, ..., iN}？这可以使用奇异值分解来完成吗？代码 sn-ps 将非常受欢迎。

编辑： 为了使这一点更具体，请考虑以下 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 线性独立的列向量的索引。

Answer 1

在 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

顺便说一句，使用 det 来确定矩阵是否为奇异矩阵是绝对、肯定、绝对最糟糕的方法。

改用排名。

基本上，您几乎不应该在 MATLAB 中使用 det，除非您明白为什么这样做是一件糟糕的事情，并且尽管如此，您还是选择使用它。

Answer 2

我的第一个想法是尝试 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 将是一个很好的询问地点。

如果您只需要一种计算方式，上面的代码就足够了。