¿Cómo puedo obtener las mismas soluciones ‘especiales’ para sistemas lineales no determinados que el operador `A \ b` (mldivide) de Matlab devuelve usando numpy / scipy?

Encontré un enlace donde se muestra con un ejemplo de que el operador Matlab mldivide ( \ ) ofrece soluciones “especiales” cuando el sistema de ecuaciones lineales tiene infinitas soluciones.

Por ejemplo:

 A = [1 2 0; 0 4 3]; b = [8; 18]; c_mldivide = A \ b c_pinv = pinv(A) * b 

da la salida:

 c_mldivide = 0 4 0.66666666666667 c_pinv = 0.918032786885245 3.54098360655738 1.27868852459016 

La solución es “especial” en el sentido de que el número de entradas distintas de cero en la solución c_mldivide es igual al rank(A) (en este caso 2). Intenté lo mismo en numpy usando numpy.linalg.lstsq , lo que da un resultado idéntico a c_pinv .

¿Hay alguna manera de lograr la solución c_mldivide en Python?

Había otra pregunta muy similar aquí , pero supongo que la explicación de la palabra “especial” no era lo suficientemente clara. Otra pregunta fue sobre el funcionamiento interno del operador mldivide , pero la respuesta aceptada no parece abordar este comportamiento.

Edición 1: código numpy

 In [149]: test_A = np.array([[1,2,0],[0,4,3]]) test_b = np.array([[8],[18]]) np.linalg.lstsq(test_A,test_b) Out[149]: (array([[ 0.918 ], [ 3.541 ], [ 1.2787]]), array([], dtype=float64), 2, array([ 5.2732, 1.4811])) 

Edición 2: utilizando scipy.optimize.nnls

 In[189]: from scipy.optimize import nnls nnls(test_A,test_b) Out[190]: ValueError Traceback (most recent call last)  in () 1 from scipy.optimize import nnls 2 ----> 3 nnls(test_A,test_b) C:\Users\abhishek\Anaconda\lib\site-packages\scipy\optimize\nnls.py in nnls(A, b) 43 raise ValueError("expected matrix") 44 if len(b.shape) != 1: ---> 45 raise ValueError("expected vector") 46 47 m, n = A.shape ValueError: expected vector 

Los mínimos cuadrados no negativos ( scipy.optimize.nnls ) no son una solución general para este problema. Un caso trivial en el que fallará es si todas las soluciones posibles contienen coeficientes negativos:

 import numpy as np from scipy.optimize import nnls A = np.array([[1, 2, 0], [0, 4, 3]]) b = np.array([-1, -2]) print(nnls(A, b)) # (array([ 0., 0., 0.]), 2.23606797749979) 

En el caso de que A · x = b esté indeterminado,

 x1, res, rnk, s = np.linalg.lstsq(A, b) 

Escogerá una solución x ‘ que minimice || x || L2 sujeto a || A · xb || L2 = 0 . Esto sucede no para ser la solución particular que estamos buscando, pero podemos transformarla linealmente para obtener lo que queremos. Para hacer eso, primero calcularemos el espacio nulo correcto de A , que caracteriza el espacio de todas las soluciones posibles para A · x = b . Podemos obtener esto usando una descomposición de QR que revela rango :

 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() Z = qr_null(A) 

Z es un vector (o, en el caso de que n – rnk ( A )> 1 , un conjunto de vectores base que abarcan un subespacio de A ) tal que A · Z = 0 :

 print(A.dot(Z)) # [[ 0.00000000e+00] # [ 8.88178420e-16]] 

En otras palabras, la (s) columna (s) de Z son vectores que son ortogonales a todas las filas en A. Esto significa que para cualquier solución x ‘ a A · x = b , entonces x’ = x + Z · c también debe ser una solución para cualquier factor de escala arbitrario c . Esto significa que al elegir un valor apropiado de c , podemos establecer cualquier n – rnk ( A ) de los coeficientes en la solución a cero.

Por ejemplo, digamos que queremos establecer el valor del último coeficiente en cero:

 c = -x1[-1] / Z[-1, 0] x2 = x1 + Z * c print(x2) # [ -8.32667268e-17 -5.00000000e-01 0.00000000e+00] print(A.dot(x2)) # [-1. -2.] 

El caso más general donde n – rnk ( A ) ≤ 1 es un poco más complicado:

 A = np.array([[1, 4, 9, 6, 9, 2, 7], [6, 3, 8, 5, 2, 7, 6], [7, 4, 5, 7, 6, 3, 2], [5, 2, 7, 4, 7, 5, 4], [9, 3, 8, 6, 7, 3, 1]]) x_exact = np.array([ 1, 2, -1, -2, 5, 0, 0]) b = A.dot(x_exact) print(b) # [33, 4, 26, 29, 30] 

Obtenemos x ‘ y Z como antes:

 x1, res, rnk, s = np.linalg.lstsq(A, b) Z = qr_null(A) 

Ahora, para maximizar el número de coeficientes de valor cero en el vector de solución, queremos encontrar un vector C tal que

x ‘ = x + Z · C = [x’ 0 , x ‘ 1 , …, x’ rnk (A) -1 , 0, …, 0] T

Si los últimos n – rnk ( A ) coeficientes en x ‘ han de ser ceros, esto impone que

Z {rnk (A), …, n} · C = – x {rnk (A), …, n}

Por lo tanto, podemos resolver para C (exactamente, ya que sabemos que Z[rnk:] debe ser de rango completo):

 C = np.linalg.solve(Z[rnk:], -x1[rnk:]) 

y calcular x ‘ :

 x2 = x1 + Z.dot(C) print(x2) # [ 1.00000000e+00 2.00000000e+00 -1.00000000e+00 -2.00000000e+00 # 5.00000000e+00 5.55111512e-17 0.00000000e+00] print(A.dot(x2)) # [ 33. 4. 26. 29. 30.] 

Para ponerlo todo junto en una sola función:

 import numpy as np from scipy.linalg import qr def solve_minnonzero(A, b): x1, res, rnk, s = np.linalg.lstsq(A, b) if rnk == A.shape[1]: return x1 # nothing more to do if A is full-rank Q, R, P = qr(AT, mode='full', pivoting=True) Z = Q[:, rnk:].conj() C = np.linalg.solve(Z[rnk:], -x1[rnk:]) return x1 + Z.dot(C) 
 np.array([[8],[18]]).shape 

es

 (2,1) 

pero tu quieres

 (2,) #!/usr/bin/env python3 import numpy as np from scipy.optimize import nnls test_A = np.array([[1,2,0],[0,4,3]]) try: test_b = np.array([[8],[18]]) # wrong print(nnls(test_A,test_b)) except Exception as e: print(str(e)) test_b = np.array([8,18]) # sic! print(nnls(test_A,test_b)) 

salida:

 expected vector (array([ 0. , 4. , 0.66666667]), 0.0)