python - 两个不同大小的数组的乘积求和

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

可以通过将 newaxis 广播到数组来轻松矢量化。因此，瓶颈 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)

正如您所看到的，这似乎是极少数情况之一，对代码进行矢量化不会提高性能。 scipy 的 norm 模块中似乎还存在另一个瓶颈，这会影响在多维数组上操作时的性能。

因此，提高代码性能的唯一可能是实现循环的并行执行(将 += 运算符替换为分配给固定数组并随后求和) .