Left Matrix Division y Numpy Solve

Question

Estoy tratando de convertir el código que contiene el operador \ de Matlab (Octave) a Python. Código de muestra

B = [2;4] 
b = [4;4] 
B \ b 

Esto funciona y produce 1.2 como respuesta. El uso de esta página web

http://mathesaurus.sourceforge.net/matlab-numpy.html

traduje que a medida:

import numpy as np 
import numpy.linalg as lin 
B = np.array([[2],[4]]) 
b = np.array([[4],[4]]) 
print lin.solve(B,b) 

Esto me dio un error:

numpy.linalg.linalg.LinAlgError: Array must be square 

¿Cómo es que Matlab \ trabaja con matriz no cuadrada de B?

¿Alguna solución para esto?

Answer 1

De MathWorks documentation para la división matriz de la izquierda:

If A is an m-by-n matrix with m ~= n and B is a column vector with m components, or a matrix with several such columns, then X = A\B is the solution in the least squares sense to the under- or overdetermined system of equations AX = B. In other words, X minimizes norm(A*X - B), the length of the vector AX - B.

El equivalente en numpy es np.linalg.lstsq:

In [15]: B = np.array([[2],[4]]) 

In [16]: b = np.array([[4],[4]]) 

In [18]: x,resid,rank,s = np.linalg.lstsq(B,b) 

In [19]: x 
Out[19]: array([[ 1.2]]) 

Answer 2

Matlab realmente hará una serie de diferentes operaciones cuando se utiliza el \ operador, en función de la forma de las matrices involucradas (ver here para más detalles). En su ejemplo, Matlab está devolviendo una solución de mínimos cuadrados, en lugar de resolver directamente la ecuación lineal, como ocurriría con una matriz cuadrada. Para obtener el mismo comportamiento en numpy, haga esto:

import numpy as np 
import numpy.linalg as lin 
B = np.array([[2],[4]]) 
b = np.array([[4],[4]]) 
print np.linalg.lstsq(B,b)[0] 

que debería darle la misma solución que Matlab.

Answer 3

Puede formar la inversa por la izquierda:

import numpy as np 
import numpy.linalg as lin 
B = np.array([[2],[4]]) 
b = np.array([[4],[4]]) 

B_linv = lin.solve(B.T.dot(B), B.T) 
c = B_linv.dot(b) 
print('c\n', c) 

Resultado:

c 
[[ 1.2]] 

En realidad, puede simplemente ejecutar el solucionador de una vez, sin que se forme una inversa, así:

c = lin.solve(B.T.dot(B), B.T.dot(b)) 
print('c\n', c) 

Resultado:

c 
[[ 1.2]] 

.... como antes

¿Por qué? Porque:

Tenemos:

Multiplicar a través de B.T, nos da:

Ahora, B.T.dot(B) es cuadrada, de rango completo, tiene una inversa. Y, por lo tanto, podemos multiplicar por el inverso de B.T.dot(B), o usar un solucionador, como se indica arriba, para obtener c.