# How to speedup integration of interpolated function with logarithmized data?

Consider some dataset corresponding to the grid x1, x2, function[x1,x2]:

function[x1_, x2_] = Exp[-x2/25]*Cos[20*x1]^2*Sin[x1];
dataset =
Flatten[Table[{x1, x2, function[x1, x2]}, {x1,
Table[10^
x, {x, -5., Log10[1/20. Pi/2], (Log10[1/20. Pi/2] + 5.)/
20}]}, {x2,
Table[10^
x, {x, Log10[5.3], Log10[200.], (Log10[200.] - Log10[5.3])/
20}]}], {1, 2}];
{x1min, x1max} = MinMax[dataset[[All, 1]]]
{x2min, x2max} = MinMax[dataset[[All, 2]]]


I interpolate it in two ways:

int1[x1_, x2_] = 10^(Interpolation[{Log10[#[]], Log10[#[]], Log10[#[] + 10^-90]}&/@dataset, InterpolationOrder -> 1][Log10[x1], Log10[x2]]);
int2[x1_, x2_] = Interpolation[{#[], #[], #[]}&/@dataset, InterpolationOrder -> 1][x1, x2];


The integration of the int1 is much slower than int2:

testfunction[x1_, x2_] = Sqrt[(1 - Cos[x1])]*Exp[5/x2];
NIntegrate[int1[x1, x2], {x1, x1min, x1max}, {x2, x2min, x2max},
Method -> "InterpolationPointsSubdivision"] // AbsoluteTiming
NIntegrate[int2[x1, x2], {x1, x1min, x1max}, {x2, x2min, x2max},
Method -> "InterpolationPointsSubdivision"] // AbsoluteTiming
NIntegrate[
int1[x1, x2]*testfunction[x1, x2], {x1, x1min, x1max}, {x2, x2min,
x2max}, Method ->
"InterpolationPointsSubdivision"] // AbsoluteTiming
NIntegrate[
int2[x1, x2]*testfunction[x1, x2], {x1, x1min, x1max}, {x2, x2min,
x2max}, Method ->
"InterpolationPointsSubdivision"] // AbsoluteTiming


{0.960714, 0.0146339}

{0.0011156, 0.0189933}

{1.01974, 0.000414336}

{0.0718938, 0.000659965}

On the other hand, int1 has better behavior in case if function changes fast between the neighboring values of x1, x2 in the grid.

Switching to "AdaptiveMonteCarlo" integration method somehow improves the situation, but non-negligible difference in timing (up to a few times) remains (it's funny that it slows down much the integration of int2 without testfunction).

Could you please tell me how to speed up the integration with int1?

• What does int1 have to do with int2 ? Mar 15 at 17:53
• Of course int1 is slower: In addition to evaluating an interpolation function, it has also to evaluate 10^... and Log10. Mar 15 at 18:04
• What do you actually try to achieve here? Mar 15 at 18:05
• This question seems ill-defined, because you're relying on Mathematica to give you an interpolation of a function, and it could be computing the two interpolations totally differently. One might have $1000$ times as many points as the other one. That would be my first guess as to why one is slow vs. fast. I think you need to provide more context for the question to be meaningful to others. What real data are you given to work with, and what are your options for processing it? I assume that you're finding one option is a lot slower than the other, but otherwise preferable. Mar 15 at 19:43
• @HenrikSchumacher : I would like to get the accurate interpolation in the domain where the distribution function decreases fastly (large x1 and x2), keeping also good performance. Mar 15 at 21:13

I suggest that you try to use the substitution formula:

$$\int_a^b \int_c^d f( \varphi(x),\psi(y)) \, \varphi'(x) \, \psi'(y) \, \mathrm{d}x \, \mathrm{d}y = \int_{\varphi(a)}^{\varphi(a)} \int_{\psi(c)}^{\psi(d)} f(u,v) \, \mathrm{d}u \, \mathrm{d}v.$$

So you could apply this to $$\varphi(t) = \psi(t) = \log_{10}(t)$$ or $$\varphi(t) = \psi(t) = 10^{t}$$. (I don't know at the moment which direction would be suitable for you.)

If $$\varphi'(x) \neq 0$$ and $$\psi'(y) \neq 0$$ for all $$x \in [a,b]$$ and $$y \in [c,d]$$, then you may define

$$g(u,v) := f(u,v) \, \varphi'(\varphi^{-1}(u)) \, \psi'(\psi^{-1}(v)),$$ which leads to the following identity: $$\int_a^b \int_c^d g( \varphi(x),\psi(y)) \, \mathrm{d}x \, \mathrm{d}y = \int_{\varphi(a)}^{\varphi(a)} \int_{\psi(c)}^{\psi(d)} g(u,v) \, \frac{1}{\varphi'(\varphi^{-1}(u))} \, \frac{1}{\psi'(\psi^{-1}(v)))} \, \mathrm{d}u \, \mathrm{d}v.$$

The point is to apply this to the following function $$g$$:

ip = Interpolation[Transpose[{Log10[x], Log10[y], Log10[z + 10^-90]}],
InterpolationOrder -> 1];
g[x1_, x2_] := 10^ip[x1, x2];


Method -> "InterpolationPointsSubdivision" should allow for fast integration of the integral on the right. And if I understand correctly, it is the integral on the left that you want to compute in the end.

• Glückwunsch zu den 100k! Mar 16 at 5:23
• Hey, danke! War 'ne Menge Arbeit... =) Mar 16 at 6:25