从满秩非方矩阵获取可逆方矩阵的 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} 呢？这可以通过奇异值分解来实现吗？如果能提供一些代码示例，那就太好了。

编辑：为了让这个问题更具体，考虑以下的 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

我首先想到的办法是尝试从 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 问问。

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

Answer 2

在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的前四列表示要使用的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）来判断一个矩阵是否是奇异矩阵，这绝对是最糟糕的方法。

应该使用秩（rank）来判断。

基本上，除非你明白为什么使用行列式是个坏主意，否则在MATLAB中几乎不要使用它，除非你是故意选择这样做。