Я пытаюсь получить LU состав матрицы. Поскольку предполагается, что scipy.linalg.lu может быть не слишком точным, я решил убедиться, что максимальная абсолютная ошибка намного меньше минимального элемента матрицы, с которой я работаю.
Вот код фрагмент, где я пытаюсь LU разложить транспонирование A:
print min([min(r[np.nonzero(r)], key=abs) for r in A], key=abs), "min"
P, L, U = scipy.linalg.lu(np.transpose(A))
B=np.matmul(np.matmul(P,L),U)
diff=np.absolute(B-np.transpose(A))
self.error.append([C, max([max(r) for r in diff])])
print max([max(r) for r in diff]), "max"
Вывод выглядит следующим образом:
1,8 мин
3.552713678800501e-15 макс
, который идеально подходит для меня. Однако есть проблема.
$$ PA ^ T = LU $$
$$ A ^ T = P ^ {- 1} LU $$
Пока что в третьей строке фрагмента кода я не умножаюсь на $ P ^ {- 1} $, я умножаю его на $ P $. По какой-то причине это работает отлично, хотя и дает мне огромные ошибки, когда я заменяю P на $ P ^ {- 1} $.
Есть что-то, чего я не знаю о np.matmul или scipy.linalg. л? Scipy.linalg.lu возвращает P или P ^ {- 1}? Или в моем коде есть ошибка?
Извинения, если этот пост не соответствует стандартам stackoverflow.