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

Question

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

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

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

Вот моя матрица образца, которая является такой же, как wikipedia example: Метод

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

(found here и here):

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

Это, кажется, работает хорошо для меня:

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

Answer 2

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

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

Answer 3

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

Answer 4

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

Answer 5

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