numpy - organizar - ordenar matriz por columna c++

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

Question

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 en 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 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, 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 <= 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 partir del año pasado (2017), scipy ahora tiene un método null_space en el módulo scipy.linalg ( docs ).

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

Answer 2

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 valores propios casi nulos, 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

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 4

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

Answer 5

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

Answer 6

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