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

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

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

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

Вот моя матрица образцов, которая совпадает с пример wikipedia:

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)
>>>

Ответ 1

Кажется, что все работает хорошо для меня:

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

Ответ 2

Sympy делает это просто.

>>> 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]])]

Ответ 3

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

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

Ответ 4

Ваш метод почти правильный. Дело в том, что форма 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.]]

Ответ 5

Более быстрый, но менее численный метод заключается в том, чтобы использовать рандомизирующую 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]]