我有一个点网格$x_i$和对应的函数值$y_i=f(x_i)$. 我对诸如 cumulant 之类的东西很感兴趣$f$，但它有一个尴尬的前因。我们将调用的所需数量

$$z(x) = e^{x}\int_{x}^{\infty}e^{-q}f(q)dq$$

每个网格点都需要这个数量，$z_i=z(x_i)$

让我们为了数字的缘故假设$\infty$真的是一个足够大的$q$这样$f(q)$顺利归零，而且$f(q)$是一个表现良好的非负函数$q$我们可能会关心。在这个特定的上下文中，它是一个概率分布。请注意，我并不总是有一个函数形式$f$，因此分析解决方案在这里并不是特别有用。

如果不是为了$e^x$prefactor，在数字上整合它会相当简单。要么自己使用梯形规则，要么scipy.integrate.cumulative_trapezoid自己做梯形规则

integrand = np.exp(-x) * y
almostZ = dx * (np.cumsum(integrand[::-1])[::-1] - integrand/2 - integrand[-1]/2)
Z = almostZ * np.exp(x)

dx的网格间距在哪里$x_i$. 这样做，你得到$z_i$大致在$O(N)$缩放，这很棒。

问题出现时$x$变得足够大$e^x$要么$e^{-x}$即使这两项的乘积是一个合理大小的数字，也违反了数值精度的限制。在循环中执行此操作是“答案”的当前状态：

for i in range(len(x)):
    Z[i] = np.trapz(np.exp(x[i] - x[i:]) * y[i:], x[i:])

这确保了没有非常大的数字乘以非常小的数字，而是合理大小的数字的总和。但它也涉及大量重复计算和$O(N^2)$-ness，当这是被计算为积分 ODE（通过）的东西时，这有点痛苦solve_ivp。

做这种计算有诀窍吗，还是我注定要这样做$O(N^2)$?