Суммирование произведений двух массивов разных размеров

Question

Мне интересно представить функцию логарифмического правдоподобия двух переменных.

Вот небольшой код для этого:

import numpy as np
from scipy.stats import norm
import matplotlib.pyplot as plt

#Generate fake data
x = np.random.normal(0,1,size = 100)
y = 2*x+1 + np.random.normal(size = x.size)

#Create a grid to visualize the log-likelihood
grid = np.linspace(-5,5,101)
B0,B1 = np.meshgrid(grid,grid)

#Compute the log likelihood
LogLik = 0
for xs,ys in zip(x,y):
    LogLik+= norm.logpdf(ys, loc = B0+B1*xs)

plt.contourf(B0,B1,LogLik)

Узким местом в этом небольшом фрагменте кода являетсявычисление логарифмической вероятности, а именно:

for xs,ys in zip(x,y):
    LogLik+= norm.logpdf(ys, loc = B0+B1*xs)

Если длина x или y очень велика, то это занимает больше времени, чем я забочусь.Есть ли способ векторизовать создание среднего распределения (то есть B0+B1*xs), а также оценку logpdf?

Answer 1

Это можно легко векторизовать, передавая новые оси массивам.В результате узкое место norm.logpdf будет выполнено только один раз:

log_lh = norm.logpdf(y, loc=B0[..., None] + B1[..., None] * x[None, None, :]).sum(axis=2)

# comparison with LogLik:
np.allclose(LogLik, log_lh)
# Out: True

Преобразование этого в функцию позволит рассчитать время выполнения:

def loglik(x, y, B0, B1):
    return norm.logpdf(y, loc=B0[..., None] + B1[..., None] * x[None, None, :]).sum(axis=2)

def loglik_loop(x, y, B0, B1):
    LogLik = 0
    for xs, ys in zip(x, y):
        LogLik+= norm.logpdf(ys, loc=B0+B1*xs)

%timeit loglik(x, y, B0, B1)
# Out: 94.1 ms ± 1.51 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)
%timeit loglik_loop(x, y, B0, B1)
# Out: 54 ms ± 4.25 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)

Как видите,похоже, это один из редких случаев, когда векторизация кода не улучшает производительность .Кажется, в модуле norm модуля scipy есть еще одно узкое место, которое снижает производительность при работе с многомерными массивами.

В результате ваша единственная возможность улучшить производительность вашего кода будет заключаться в следующем:реализовать параллельное выполнение цикла (замена оператора += назначением фиксированного массива и последующим суммированием).