Получение обратимой квадратной матрицы из неквадратной матрицы полного ранга в numpy или matlab - PullRequest
Новая
       17

Получение обратимой квадратной матрицы из неквадратной матрицы полного ранга в numpy или matlab

4 голосов
/ 22 августа 2010

Предположим, у вас есть 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? Можно ли это сделать с помощью разложения по сингулярным числам? Фрагменты кода будут очень кстати.

EDIT: Чтобы сделать это более конкретным, рассмотрим следующий код 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 линейно независимых векторов столбцов.

Ответы [ 2 ]

6 голосов
/ 22 августа 2010

Легко, peasy в 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

Первые 4 столбца 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 для определения того, является ли матрица единственной, является абсолютно, положительно, определенно худшим способом сделать это.

Вместо этого используйте ранг.

По сути, вы практически НИКОГДА не должны использовать det в MATLAB, если вы не понимаете, почему это так плохо, и вы решаете использовать его, несмотря на этот факт.

1 голос
/ 22 августа 2010

Моей первой мыслью было бы попробовать каждую возможную комбинацию N из M столбцов.Это можно сделать следующим образом (в 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 будет хорошим местом для вопроса.

Если все, что вам нужно, это способ сделать это вычислительно, приведенного выше кода должно быть достаточно.

...