Python (NumPy, SciPy), нахождение нулевого пространства матрицы

Question

Я пытаюсь найти нулевое пространство (пространство решения Ax = 0) данной матрицы.Я нашел два примера, но я не могу заставить их работать.Более того, я не могу понять, что они делают, чтобы попасть туда, поэтому я не могу отлаживать.Я надеюсь, что кто-то сможет провести меня через это.

Страницы документации (numpy.linalg.svd и numpy.compress) для меня непрозрачны.Я научился делать это, создавая матрицу C = [A|0], находя сокращенную форму ряда строк и решая переменные по строкам.Я не могу понять, как это делается в этих примерах.

Спасибо за любую помощь!

Вот мой пример матрицы, такой же, как Википедияпример :

A = matrix([
    [2,3,5],
    [-4,2,3]
    ])  

Метод ( найден здесь и здесь ):

import scipy
from scipy import linalg, matrix
def null(A, eps=1e-15):
    u, s, vh = scipy.linalg.svd(A)
    null_mask = (s <= eps)
    null_space = scipy.compress(null_mask, vh, axis=0)
    return scipy.transpose(null_space)

Когда я пытаюсь это сделать, явернуть пустую матрицу:

Python 2.6.6 (r266:84292, Sep 15 2010, 16:22:56) 
[GCC 4.4.5] on linux2
Type "help", "copyright", "credits" or "license" for more information.
>>> import scipy
>>> from scipy import linalg, matrix
>>> def null(A, eps=1e-15):
...    u, s, vh = scipy.linalg.svd(A)
...    null_mask = (s <= eps)
...    null_space = scipy.compress(null_mask, vh, axis=0)
...    return scipy.transpose(null_space)
... 
>>> A = matrix([
...     [2,3,5],
...     [-4,2,3]
...     ])  
>>> 
>>> null(A)
array([], shape=(3, 0), dtype=float64)
>>> 

Answer 1

Симпи делает это просто.

>>> from sympy import Matrix
>>> A = [[2, 3, 5], [-4, 2, 3], [0, 0, 0]]
>>> A = Matrix(A)
>>> A * A.nullspace()[0]
Matrix([
[0],
[0],
[0]])
>>> A.nullspace()
[Matrix([
[-1/16],
[-13/8],
[    1]])]

Answer 2

Вы получаете SVD-разложение матрицы A. s - вектор собственных значений. Вас интересуют почти нулевые собственные значения (см. $ A * x = \ lambda * x $, где $ \ abs (\ lambda) <\ epsilon $), который задается вектором логических значений <code>null_mask.

Затем вы извлекаете из списка vh собственные векторы, соответствующие почти нулевым собственным значениям, и это именно то, что вы ищете: способ заполнить нулевое пространство. По сути, вы извлекаете строки, а затем транспонируете результаты, чтобы получить матрицу с собственными векторами в виде столбцов.

Answer 3

Кажется, у меня все в порядке:

A = matrix([[2,3,5],[-4,2,3],[0,0,0]])
A * null(A)
>>> [[  4.02455846e-16]
>>>  [  1.94289029e-16]
>>>  [  0.00000000e+00]]

Answer 4

Более быстрый, но менее численно устойчивый метод заключается в использовании ранжирующего разложения QR , такого как scipy.linalg.qr с pivoting=True:

import numpy as np
from scipy.linalg import qr

def qr_null(A, tol=None):
    Q, R, P = qr(A.T, mode='full', pivoting=True)
    tol = np.finfo(R.dtype).eps if tol is None else tol
    rnk = min(A.shape) - np.abs(np.diag(R))[::-1].searchsorted(tol)
    return Q[:, rnk:].conj()

Например:

A = np.array([[ 2, 3, 5],
              [-4, 2, 3],
              [ 0, 0, 0]])
Z = qr_null(A)

print(A.dot(Z))
#[[  4.44089210e-16]
# [  6.66133815e-16]
# [  0.00000000e+00]]

Answer 5

По состоянию на прошлый год (2017), у scipy теперь есть встроенный метод null_space в модуле scipy.linalg ( документы ).

Реализация следует за канонической декомпозицией SVD и довольно мала, если у вас более старая версия scipy и вам нужно реализовать ее самостоятельно (см. Ниже).Однако, если вы в курсе, это для вас.

def null_space(A, rcond=None):
    u, s, vh = svd(A, full_matrices=True)
    M, N = u.shape[0], vh.shape[1]
    if rcond is None:
        rcond = numpy.finfo(s.dtype).eps * max(M, N)
    tol = numpy.amax(s) * rcond
    num = numpy.sum(s > tol, dtype=int)
    Q = vh[num:,:].T.conj()
    return Q

Answer 6

Ваш метод почти правильный.Проблема заключается в том, что форма s, возвращаемая функцией scipy.linalg.svd, имеет вид (K,), где K = min (M, N).Таким образом, в вашем примере s имеет только две записи (сингулярные значения первых двух сингулярных векторов).Следующее исправление вашей нулевой функции должно позволить ей работать для матрицы любого размера.

import scipy
import numpy as np
from scipy import linalg, matrix
def null(A, eps=1e-12):
...    u, s, vh = scipy.linalg.svd(A)
...    padding = max(0,np.shape(A)[1]-np.shape(s)[0])
...    null_mask = np.concatenate(((s <= eps), np.ones((padding,),dtype=bool)),axis=0)
...    null_space = scipy.compress(null_mask, vh, axis=0)
...    return scipy.transpose(null_space)
A = matrix([[2,3,5],[-4,2,3]])
print A*null(A)
>>>[[  4.44089210e-16]
>>> [  6.66133815e-16]]
A = matrix([[1,0,1,0],[0,1,0,0],[0,0,0,0],[0,0,0,0]])
print null(A)
>>>[[ 0.         -0.70710678]
>>> [ 0.          0.        ]
>>> [ 0.          0.70710678]
>>> [ 1.          0.        ]]
print A*null(A)
>>>[[ 0.  0.]
>>> [ 0.  0.]
>>> [ 0.  0.]
>>> [ 0.  0.]]