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

Estoy tratando de encontrar el espacio nulo (espacio de 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] , buscando 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 cualquier y toda la ayuda!

Aquí está mi matriz de muestra, que es la misma que el ejemplo de wikipedia :

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

Método ( encontrado aquí , y aquí ):

 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 bash, recupero una matriz vacía:

 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 >> A = matrix([ ... [2,3,5], ... [-4,2,3] ... ]) >>> >>> null(A) array([], shape=(3, 0), dtype=float64) >>> 

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

Sympy lo hace 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]])] 

Obtienes la descomposición SVD de la matriz A s es un vector de valores propios. Usted está interesado en valores propios casi cero (vea $ A * x = \ lambda * x $ donde $ \ abs (\ lambda) < \ epsilon $), que viene dado 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.

Tu método es casi correcto. 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.]] 

Un método más rápido pero menos numéricamente estable es usar una descomposición QR que revela el rango , 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(AT, 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]] 

A partir del año pasado (2017), scipy ahora tiene un método null_space en el módulo scipy.linalg ( docs ).

La implementación sigue la descomposición canónica de SVD y es bastante pequeña si tiene una versión anterior de scipy y necesita implementarla usted mismo (consulte a continuación). Sin embargo, si está actualizado, está ahí para usted.

 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