4

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

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

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

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

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

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

4

2 回答 2

6

在 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 ，并且尽管如此，您还是选择使用它。

于 2010-08-22T03:12:43.357 回答
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=AQ 的运算）。但是，如果您想知道是否存在，math.stackexchange.com 将是一个很好的询问地点。

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

于 2010-08-21T22:00:19.877 回答