2
$\begingroup$

I have a grid of points $x_i$ and corresponding function values $y_i=f(x_i)$. I'm interesting in something like the cumulant of $f$, but it has an awkward prefactor. The desired quantity we'll call

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

And this quantity is needed for each grid point, $z_i=z(x_i)$

Let us assume for numerical sake that $\infty$ is really a sufficiently large $q$ such that $f(q)$ goes smoothly to zero, and that $f(q)$ is a nicely behaved, non-negative function for any $q$ we might care about. In this particular context, its a probability distribution. Note that I don't always have a functional form for $f$, so an analytical solution isn't particularly helpful here.

Were it not for the $e^x$ prefactor, integrating this would be rather simple, numerically. Either use scipy.integrate.cumulative_trapezoid or do the trapezoidal rule yourself by taking

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

where dx is the grid spacing for the $x_i$. Done this way, you get $z_i$ in roughly $O(N)$ scaling, which is great.

The problem arises when $x$ gets sufficiently large that either $e^x$ or $e^{-x}$ are up against the limits of numerical precision, even though the product of the two terms is a reasonably sized number. Doing this in a loop is the current state of the "answer":

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

This ensures that there's no really huge numbers multiplied by really small numbers, and is instead a sum of reasonably sized numbers. But it also involves a lot of double counting and $O(N^2)$-ness, which is a bit of a pain when this is something that's being calculated as integrating an ODE (via solve_ivp).

Is there a trick to doing this sort of computation, or am I just doomed to $O(N^2)$?

$\endgroup$
2
  • 1
    $\begingroup$ Would arbitrary precision arithmetic help? The individual operations would likely be slower, but it would allow you to use your proposed O(N) algorithm. I suppose it would depend how large a system you are considering whether this would be worth it. $\endgroup$
    – Tyberius
    Jan 11 at 2:32
  • $\begingroup$ @Tyberius I suppose that's the obvious answer. I totally didn't think of that, and lo and behold it actually helps quite a bit. The overhead for the arbitrary precision hurts a little at small system sizes but at the number of points I need for $x_i$ its a 30% increase in speed! $\endgroup$ Jan 12 at 21:36

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.