如何以最快的方式计算一堆指数函数的乘积？

Question

这是 Fortran 代码。我想计算一个名为 pYq_i 的函数，它是 mi 指数函数的乘积，见下文，

左边 (LHS) 是 pYq_i，RHS 是它的表达式。在代码中，theta(1) 是 ki，theta(2) 是 Di。然而，ki 和 Di 是随机变量。 sigma, D, 只是常数，它们可以是任何值。 Yji 是数据的值，tj 是时间。例如，t1 - t5 是 0.1,0.2,0.3,0.4,0.5.

这是代码，

function pYq_i(theta,i)
integer(kind=i8) :: i,j
real(kind=r8) :: theta(dim_p),mean,sigma,pYq_i,fact_mean
fact_mean = D/theta(2)
pYq_i = normpdf_factor_pYq_i
do j = 1,mi
    mean = fact_mean*exp(-theta(1)*t(j))
    sigma = abs(sig*mean)   
    pYq_i = pYq_i/sigma*exp(-half*((Yji(j,i)-mean)/sigma)**2)   
enddo  
return
end function pYq_i

我想以最快的速度计算这个 pYq_i。

我首先想到的是我可以做一些事情（不是静态的，这里只是说明），

mean(:) = xxxxxxx
sigma(:) = xxxxxx
pYq(:) = 1.0/sigma(:)*exp(-half*((Yji(:,i)-mean(:))/sigma(:))**2)
pYq_i = product(pYq)

但是上面的代码比do循环代码慢。我想在 do 循环中，事情会被保存在缓存或其他东西中。

所以我做了个小把戏，我重写了do循环中的pYq_i表达式，这样我现在计算指数的对数，然后求和，所以我在循环中保存了一些exp计算.最后我只需要做一次exp。看下面的代码，pYq_i_detail和pYq_i,

完全一样

function pYq_i_detail(theta,i) ! This should be faster than pYq_i.
integer(kind=i8) :: i,j
real(kind=r8) :: theta(dim_p),sigma_inv,pYq_i_detail,product_sigma_inv,log_pYq_i,fact_mean,xk
xk = theta(1)
fact_mean = theta(2)/D
log_pYq_i = zero
product_sigma_inv = normpdf_factor_pYq_i 
do j = 1,mi
    sigma_inv = sig_inv*fact_mean*exp(xk*t(j))
    log_pYq_i = log_pYq_i - (Yji(j,i)*sigma_inv-sig_inv)**2 
    product_sigma_inv = product_sigma_inv*sigma_inv
enddo
pYq_i_detail = abs(product_sigma_inv)*exp(half*log_pYq_i)   
return
end function pYq_i_detail

顺便说一句，在上面的所有代码中，如果有任何变量没有在函数内定义，那么它就是在函数外定义的。它们只是变量，你可以给它们任何你想要的值。比如normpdf_factor_pYq_i，就是1/((2*pi)**mi)。 sig_inv 也只是一个常数。

单次计算成本低。但是，问题在于必须在代码中多次评估此函数。因此这个功能的评估似乎是一个瓶颈和关键。

现在我想问一下，有没有人对计算这个 pYq_i 的更快版本有建议？谢谢！

Answer 1

正如 kvantour 指出的那样，由于分子中的指数是平方的，因此您可以将术语 abs( D/Vi*exp(-ki*t(j)) ) 替换为 D/Vi*exp(-ki*t(j))。这允许将此项分解出来，给出更简单的指数。

然后，用总和的指数替换指数的乘积，并尽可能多地退出循环，我得到

function p(Yi,t,Vi,ki,D,sigma,mi)
  real(r8), intent(in) :: Yi(:)
  real(r8), intent(in) :: t(:)
  real(r8), intent(in) :: Vi
  real(r8), intent(in) :: ki
  real(r8), intent(in) :: D
  real(r8), intent(in) :: sigma
  integer,  intent(in) :: mi
  real(r8)             :: p
  
  real(r8) :: factor
  
  factor = D/Vi
  
  p = exp( -1/(2*sigma**2) * sum((Yi*exp(ki*t)/factor-1)**2) + ki*sum(t) ) &
  & / (sigma*abs(factor)*sqrt(2*pi))**mi
end function

我认为这与因式分解的速度一样快。为了更快地进行，您需要深入研究优化，并为此提供帮助，我们理想情况下需要 minimum example 典型用例。

如何以最快的方式计算一堆指数函数的乘积？

How to calculate the product of a bunch of exponential functions in the fastest way possible?

fortran