Ускорить/распараллелить multivariate_normal.pdf
У меня есть несколько точек Nx3, и я последовательно генерирую новое значение для каждой из соответствующего многомерного гауссиана, каждая со средним значением 1x3 и cov 3x3. Итак, у меня есть массивы: массив точек Nx3, массив средних Nx3 и массив covs Nx3x3.
Я вижу только, как это сделать с помощью классического цикла for:
import numpy as np
from scipy.stats import multivariate_normal
# Generate example data
N = 5 # Small number for minimal example, can be increased for real use case
points = np.random.rand(N, 3)
means = np.random.rand(N, 3)
covs = np.array([np.eye(3) for _ in range(N)]) # Identity matrices as example covariances
# Initialize an array to store the PDF values
pdf_values = np.zeros(N)
# Loop over each point, mean, and covariance matrix
for i in range(N):
pdf_values[i] = multivariate_normal.pdf(points[i], mean=means[i], cov=covs[i])
print("Points:\n", points)
print("Means:\n", means)
print("Covariances:\n", covs)
print("PDF Values:\n", pdf_values)
Есть ли способ ускорить это? Я пытался передать все непосредственно в multivariate_normal.pdf, а также из документации, которая, похоже, не поддерживается (в отличие от более простого случая генерации значений для точек Nx3, но с тем же средним значением и ковариацией.
Может быть какая-то реализация не от scipy?
Возможно, я слишком надеялся, но почему-то я надеюсь, что есть более простой способ ускорить это и избежать итерации с этим циклом for в большом массиве данных непосредственно с помощью цикла Pythonic.
@JérômeRichard конечно, решением для резервного копирования было бы переписать его с помощью некоторых инструментов/языков, которые вы упомянули, просто я едва ли единственный, кому это нужно, поэтому я понял, что кто-то, возможно, уже сделал это :)
Ответы 1
Вот решение, включающее разложение Холецкого.
метод 1:
import numpy as np
x = points - means
p = x.shape[1]
res = np.exp(-0.5*(x*np.linalg.solve(covs, x)).sum(1))
res/(2*np.pi)**(p/2)/np.linalg.det(covs)**0.5
array([0.03356053, 0.03167125, 0.08042358, 0.04351325, 0.1328082 ])
метод 2:
y = np.linalg.solve(LU := np.linalg.cholesky(covs), points - means)
np.exp(-0.5* (y**2+np.log(2*np.pi)).sum(1))/np.einsum('kii->ki',LU).prod(1)
array([0.03356053, 0.03167125, 0.08042358, 0.04351325, 0.1328082 ])
метод 3:
multivariate_normal.pdf(y, np.zeros(p))/np.einsum('kii->ki',LU).prod(1)
array([0.03356053, 0.03167125, 0.08042358, 0.04351325, 0.1328082 ])
метод 4: Использование svd
U,S,V = np.linalg.svd(covs)
res_svd = np.einsum('ki,kij,kj,kjl,kl->k', x, U, 1/S, V, x)
np.exp(-0.5 * (res_svd + np.log(2 * np.pi * S).sum(1)))
array([0.03356053, 0.03167125, 0.08042358, 0.04351325, 0.1328082 ])
Сгенерировать данные:
import numpy as np
from scipy.stats import multivariate_normal, wishart
# Generate example data
N = 5 # Small number for minimal example, can be increased for real use case
p = 3 # dimension
np.random.seed(10)
points = np.random.rand(N, p)
means = np.random.rand(N, p)
covs = wishart.rvs(p, np.eye(p),N, 2) # Identity matrices as example covariances
# Initialize an array to store the PDF values
pdf_values = np.zeros(N)
# Loop over each point, mean, and covariance matrix
for i in range(N):
pdf_values[i] = multivariate_normal.pdf(points[i], mean=means[i], cov=covs[i])
pdf_values
array([0.03356053, 0.03167125, 0.08042358, 0.04351325, 0.1328082 ])
Отличная работа! Обратите внимание, что использование optimize=True может немного ускорить работу np.einsum. Метод 2 кажется самым быстрым на моей машине для N=500 (с входами OP).
@JérômeRichard, сначала ты мог бы даже оптимизировать с помощью np.einsum_path. Также я ожидаю, что метод 2 будет самым быстрым. 3-й раствор дважды уничтожает холестерин
Это потрясающе и именно то, что я искал. Я взял версию номер 2. Обратите внимание, что по каким-то причинам с numpy 2.0 есть проблемы, как и со многими другими вещами — я решил придерживаться numpy 1.X, и это сработало.
@ Валерия, ты можешь использовать scipy.linalg. хотя не уверен, что он поддерживает 3D. Также, чтобы выжать больше скорости, используйте numba с scipy.linalg.lapack.
@Onyambu на самом деле, эта и некоторые другие оптимизации, которые я сделал для кода, помогли настолько, что теперь около 90-95% времени выполнения занимает просто чтение вводимых данных... :D Так что сейчас я сосредоточусь на этом, но спасибо за подсказки для дальнейших действий.
pdfизвестно, что функции Scipy довольно медленные. Переписывание его на родном языке (например, C/C++) или с помощью таких инструментов, как Cython или Numba, часто помогает значительно повысить производительность, особенно когда код оптимизирован под ваши конкретные нужды (вот пример). Это также помогает легко распараллелить цикл. Кроме того, обратите внимание, что Nx3x3, как правило, не является лучшим расположением памяти с точки зрения производительности. 3x3xN часто лучше.