Мне нужно запустить этот код тысячи раз подряд:
def update_posterior(y, x, prior_mu, prior_V, prior_a, prior_b):
# Building blocks used to keep following calculation cleaner
prior_cov_inverse = np.linalg.inv(prior_V)
x_transpose = x.transpose()
n = len(y)
residuals = y - np.dot(x, prior_mu.transpose())
# Calculation of posterior parameters
V_posterior = np.linalg.inv((prior_cov_inverse + x_transpose * x))
mu_posterior = (V_posterior * (prior_cov_inverse * prior_mu.transpose() + x_transpose * y)).transpose()
a_posterior = prior_a + n/2
b_posterior = np.asscalar(prior_b + (residuals.transpose() * np.linalg.inv((np.identity(n) + x * prior_V * x_transpose)) * residuals)/2)
return mu_posterior, V_posterior, a_posterior, b_posterior
То, как он работает, заключается в том, что вывод функции возвращается в него, и mu_posterior становится prior_mu,V_posterior становится prior_V, a_posterior становится prior_a, а b_posterior становится prior_b для следующего вызова.Y и x различны в каждом вызове.
Эта функция ужасно медленная - мне требуется около 8 секунд для запуска.Это из-за масштаба.У меня ~ 5000 параметров, поэтому prior_mu равно (1, 5000), prior_V равно (5000,5000) и симметрично positive-definite, а prior_a и prior_b - скаляры.у - скаляр, а х - (1, 5000).
Ниже приводится разбивка по времени на строку:
3.75s: prior_cov_inverse = np.linalg.inv(prior_V)
3.86s: V_posterior = np.linalg.inv((prior_cov_inverse + x_transpose * x))
0.13s: b_posterior = np.asscalar(prior_b + (residuals.transpose() * np.linalg.inv((np.identity(n) + x * prior_V * x_transpose)) * residuals)/2)
Есть идеи, как мне это ускорить?Я пытался использовать разложение Холецкого, но оно еще медленнее ?!Я предполагаю, что есть более эффективный способ реализовать декомпозицию Холецкого в Python.
prior_cov_inverse2 = np.linalg.inv(np.linalg.cholesky(prior_V))
prior_cov_inverse2 = np.dot(prior_cov_inverse2.transpose(), prior_cov_inverse2)
РЕДАКТИРОВАТЬ: Вот несколько примеров данных, чтобы проиллюстрировать проблему ....
import numpy as np
prior_mu = np.asmatrix(np.full((1, 5040), 5))
prior_V = np.diagflat(np.asmatrix(np.full((1, 5040), 30))) #usually not diagonal, but always symmetric positive definitive
a = 2
b = 2
y = np.asmatrix([10])
x = np.asmatrix(np.concatenate(([1], np.zeros(5039))))
print(update_posterior(y, x, prior_mu, prior_V, a, b))
РЕДАКТИРОВАТЬII:
Мне удалось снизить это значение с ~ 8 с / прогон до ~ 1,4 с / прогон, убрав инверсии матриц в пользу решений, а также используя формулу Шермана Моррисона.Вот мой текущий код.Если у кого-то есть идеи, как его ускорить, поделитесь!:)
def update_posterior(y, x, prior_mu, prior_V, prior_a, prior_b, I):
# Building blocks used to keep following calculation cleaner
x_transpose = x.transpose()
n = len(y)
residuals = y - np.dot(x, prior_mu.transpose())
# Calculation of posterior parameters
# Below is equivalent to np.linalg.inv(prior_V_inverse + np.dot(x_transpose, x)) but significantly faster
V_posterior = prior_V - np.true_divide(np.linalg.multi_dot((prior_V, x_transpose, x, prior_V)), 1 + np.matmul(np.matmul(x, prior_V), x_transpose))
# Below is equivalent to mu_posterior = np.dot(V_posterior, (np.matmul(prior_V_inverse, prior_mu.transpose()) + np.matmul(x_transpose, y))).transpose() but significantly faster
mu_posterior = np.dot(V_posterior, np.linalg.solve(prior_V, prior_mu.transpose()) + np.matmul(x_transpose, y)).transpose()
a_posterior = prior_a + n/2
b_posterior = np.asscalar(prior_b + (np.matmul(np.matmul(residuals.transpose(), np.linalg.inv((np.identity(n) + np.matmul(np.matmul(x, prior_V), x_transpose)))), residuals))/2)
return mu_posterior, V_posterior, a_posterior, b_posterior