Python (NumPy, SciPy), encontrar el espacio nulo de una matriz

Question

Estoy tratando de encontrar el espacio nulo (espacio solución de Ax = 0) de una matriz dada. Encontré dos ejemplos, pero parece que tampoco puedo trabajar. Además, no puedo entender lo que están haciendo para llegar allí, así que no puedo depurar. Espero que alguien pueda ayudarme a superar esto.

Las páginas de documentación (numpy.linalg.svd y numpy.compress) son opacas para mí. Aprendí a hacer esto creando la matriz C = [A|0], encontrando la forma de escala escalonada reducida y resolviendo las variables por fila. Parece que no puedo seguir cómo se está haciendo en estos ejemplos.

¡Gracias por toda la ayuda!

Aquí es mi matriz de la muestra, que es el mismo que el wikipedia example: Método

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

(found here y 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) 

Cuando lo intento, le regreso una lata vacía matriz:

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

Parece estar funcionando bien para mí:

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

Obtiene la descomposición SVD del matriz A. s es un vector de valores propios. Usted está interesado en casi sin valores propios (véase $ A * x = \ lambda * x $ donde $ \ ABS (\ lambda) < \ epsilon $), que viene dada por el vector de valores lógicos null_mask.

Luego, extrae de la lista vh los vectores propios correspondientes a los autovalores casi cero, que es exactamente lo que está buscando: una forma de abarcar el espacio nulo. Básicamente, extrae las filas y luego transpone los resultados para que obtenga una matriz con vectores propios como columnas.

Answer 3

Su método es casi correcta. El problema es que la forma de s devuelta por la función scipy.linalg.svd es (K,) donde K = min (M, N). Por lo tanto, en su ejemplo, s solo tiene dos entradas (los valores singulares de los primeros dos vectores singulares). La siguiente corrección a su función nula debería permitirle trabajar para cualquier matriz de tamaño.

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 hace esto sencillo.

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

Un método más rápido pero menos estable numéricamente es utilizar a rank-revealing QR decomposition, como scipy.linalg.qr con 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() 

Por ejemplo:

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