Управление левой матрицей и решение Numpy

Я пытаюсь преобразовать код, содержащий оператор\из Matlab (Octave) в Python. Пример кода

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

Это работает и дает 1.2 в качестве ответа. Использование этой веб-страницы

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

Я перевел это как:

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

Это дало мне ошибку:

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

Почему Matlab\работает с не квадратной матрицей для B?

Любые решения для этого?

Ответы

Ответ 1

От Документация MathWorks для разделения левой матрицы:

Если A - матрица m-by-n с m ~ = n и B - вектор-столбец с m компонентами или матрицей с несколькими такими столбцами, то X = A\B является решение в смысле наименьших квадратов к недоопределенному или переопределенному система уравнений AX = B. Другими словами, X минимизирует норму (A * X - B), длина вектора AX - B.

Эквивалент в numpy равен 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]])

Ответ 2

Matlab фактически выполнит ряд различных операций, когда используется оператор\в зависимости от формы используемых матриц (подробнее см. здесь). В вашем примере Matlab возвращает решение наименьших квадратов, а не прямое решение линейного уравнения, как это происходит с квадратной матрицей. Чтобы получить такое же поведение в numpy, сделайте следующее:

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]

который должен дать вам то же решение, что и Matlab.

Ответ 3

Вы можете сформировать левый обратный:

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)

Результат:

c
 [[ 1.2]]

Собственно, мы можем просто запустить решатель один раз, не образуя обратного, например:

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

Результат:

c
 [[ 1.2]]

.... как и раньше

Почему? Потому что:

Имеем:

Умножьте на B.T, дайте нам:

Теперь B.T.dot(B) является квадратным, полный ранг имеет инверсию. И поэтому мы можем умножить на инверсию B.T.dot(B) или использовать решатель, как указано выше, для получения c.