эффективность обращения матрицы в numpy с разложением Холецкого

У меня есть симметричная положительно-определенная матрица (например, матрица ковариации), и я хочу вычислить ее обратную. Что касается математики, я знаю, что для инвертирования матрицы более эффективно использовать разложение Холецкого, особенно если ваша матрица большая. Но я не был уверен, как работает «numpy.lianlg.inv ()». Скажем, у меня есть следующий код:

import numpy as np

X = np.arange(10000).reshape(100,100)
X = X + X.T - np.diag(X.diagonal()) #  symmetry 
X = np.dot(X,X.T) # positive-definite

# simple inversion:
inverse1 = np.linalg.inv(X) 

# Cholesky decomposition inversion:
c = np.linalg.inv(np.linalg.cholesky(X))
inverse2 = np.dot(c.T,c)

Какой из них более эффективен (обратный1 или обратный2)? Если второй более эффективен, почему numpy.linalg.inv () вместо этого не использует его?


python numpy performance matrix matrix-inverse
person Babak    schedule 03.06.2017    source источник
comment
Что касается вашего последнего вопроса - numpy не знает, что ваша матрица симметрична, поэтому не может использовать последний метод. Проверка симметричности матрицы будет медленной.   -  person Eric    schedule 04.06.2017
comment
Обратите внимание, что inv также не использует тот факт, что cholesky треугольный, он не использует DTRTRI лапака.   -  person Eric    schedule 04.06.2017
comment
Ваш код не запускается для меня, утверждая, что X не является положительно определенным (предположительно из-за переполнения)   -  person Eric    schedule 04.06.2017
comment
В общем, плохая идея инвертировать матрица. inv является дорогостоящим и нестабильным в числовом отношении. Обычно вы хотите умножить обратное на вектор, то есть хотите решить систему уравнений. Во всех таких случаях лучше просто решить систему, используя что-то вроде linalg.solve (указание solve, что матрица симметрична и положительно определена, заставит solve использовать Холецкого). Если вы хотите использовать обратное несколько раз, вычислите и сохраните разложение для дальнейшего использования.   -  person Praveen    schedule 04.06.2017
comment
если вы хотите инвертировать коэффициент холецкого, используйте scipy.linalg.lapack.dtrtri   -  person Bananach    schedule 23.09.2020

Ответы (1)


arrow_upward
1
arrow_downward

При следующей настройке:

import numpy as np

N = 100
X = np.linspace(0, 1, N*N).reshape(N, N)
X = 0.5*(X + X.T) + np.eye(N) * N

Я получаю следующие тайминги с IPython's %timeit:

In [28]: %timeit np.linalg.inv(X)
255 µs ± 30.9 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)

In [29]: %timeit c = np.linalg.inv(np.linalg.cholesky(X)); np.dot(c.T,c)
414 µs ± 15.4 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)
person Eric    schedule 03.06.2017