How to avoid repetitive calculation when doing numerical integral?

Question

Suppose I have a function f[x] which is very complicated, together with a function g[f[x]]+h[x] to integrate. That is:

NIntegrate[{f[x], g[f[x]]+h[x]}, {x,0,1}]

I suppose mathematica will calculate f[x] and g[f[x]]+h[x] seperately, thus calculated f[x] twice. How can I speed up the calculation by telling mathematica calculate f[x] only once?

A concrete example suggested by the comments:

ClearAll["Global`*"]

f[x_, y_, z_] := Exp[Sin[x]] + Cos[y + z];

NIntegrate[{f[x, y, z], Sqrt[f[x, y, z]] + x}, {x, 0, 10}, {y, 0, 10}, {z, 0, 10}, 
Method -> {"LocalAdaptive", "SymbolicProcessing" -> 0}, 
PrecisionGoal -> 7] // Timing

{8.85938, {1429.54, 6081.95 + 59.0571 I}}

g[x_, y_, z_] := g[x, y, z] = Exp[Sin[x]] + Cos[y + z];

NIntegrate[{g[x, y, z], Sqrt[g[x, y, z]] + x}, {x, 0, 10}, {y, 0, 10}, {z, 0, 10}, 
Method -> {"LocalAdaptive", "SymbolicProcessing" -> 0}, 
PrecisionGoal -> 7] // Timing

{8.90625, {1429.54, 6081.95 + 59.0571 I}}

In this example the memoization seems not speed up the calcuation...

Answer 1

This answer shows how to define a new NIntegrate rule that evaluates f in the list of two integrands {f[x],g[f[x]]+h[x]} only once per sampling point. The answer can be also easily modified into an answer of "NIntegrate over a list of functions".

The definition of the NIntegrate rule LessEvaluationsRule given below is also aimed to be didactic and conceptually simple. The design (options and plug-in mechanism utilization) can be made to be more robust.

How to use

In order to compute the integrals in the standard NIntegrate call:

NIntegrate[{Sin[x], Sqrt[Sin[x]] + x^2}, {x, 0, 1}]
(* {0.459698, 0.976311} *)

we specify the function of one argument Sin and the function of two arguments (Sqrt[#2] + #1^2 &):

NIntegrate[1, {x, 0, 1}, 
 Method -> {"GlobalAdaptive", "SingularityHandler" -> None, 
   Method -> {LessEvaluationsRule, "BaseFunction" -> Sin, 
     "DependantFunction" -> (Sqrt[#2] + #1^2 &)}}, PrecisionGoal -> 6]
(* {0.459698, 0.976311} *)

Note that the rule definition below does not allow the (correct) application of NIntegrate's singularity handling so it is prevented with "SingularityHandler" -> None.

Comparison

Another consequence of specifying the integrated functions as option values is that we cannot use EvaluationMonitor. We can see though that the timing is twice smaller. (Or just examine IRuleEstimate defined below.)

In[214]:= RepeatedTiming[
 NIntegrate[{Sin[x], Sqrt[Sin[x]] + x^2}, {x, 0, 1}, PrecisionGoal -> 4]]

Out[214]= {0.0062, {0.459698, 0.976315}}


In[215]:= RepeatedTiming[
 NIntegrate[1, {x, 0, 1}, 
  Method -> {"GlobalAdaptive", "SingularityHandler" -> None, 
    Method -> {LessEvaluationsRule, "BaseFunction" -> Sin, 
      "DependantFunction" -> (Sqrt[#2] + #1^2 &)}}, PrecisionGoal -> 4]
 ]

Out[215]= {0.0023, {0.459698, 0.976315}}

Definition of LessEvaluationsRule

Clear[LessEvaluationsRule];
Options[LessEvaluationsRule] = {"Method" -> "GaussKronrodRule", 
   "BaseFunction" -> Sin, "DependantFunction" -> (Sqrt[#2] + #1 &)};

LessEvaluationsRuleProperties = Part[Options[LessEvaluationsRule], All, 1];

LessEvaluationsRule /: 
  NIntegrate`InitializeIntegrationRule[LessEvaluationsRule, nfs_, ranges_, 
   ruleOpts_, allOpts_] :=

  Module[{t, methodSpec, baseFunc, depFunc, pos, absc, weights, errweights},

   t = NIntegrate`GetMethodOptionValues[LessEvaluationsRule, 
     LessEvaluationsRuleProperties, ruleOpts];
   If[t === $Failed, Return[$Failed]];
   {methodSpec, baseFunc, depFunc} = t;

   t = NIntegrate`MOptionValue[methodSpec, nfs, ranges, allOpts];

   If[t === $Failed, Return[$Failed]];
   {absc, weights, errweights} = t[[1]];

   LessEvaluationsRule[{absc, weights, errweights}, baseFunc, depFunc]
   ];

Clear[IRuleEstimate]
IRuleEstimate[f_, g_, {a_, b_}, {absc_, weights_, errweights_}] :=
  Block[{fvals, integral1, error1, integral2, error2, xs},
   xs = Rescale[absc, {0, 1}, {a, b}];
   fvals = f /@ xs;
   {integral1, error1} = (b - a) {fvals.weights, fvals.errweights};
   fvals = g @@@ Transpose[{xs, fvals}];
   {integral2, 
     error2} = (b - a) {fvals.weights, fvals.errweights}; {integral1, 
    integral2, Abs[error1], Abs[error2]}];

LessEvaluationsRule[{absc_, weights_, errweights_}, baseFunc_, depFunc_][
   "ApproximateIntegral"[region_]] := 
  Block[{a, b, fvals, f2expr, integral1, integral2, error1, error2},

   {a, b} = region["Boundaries"][[1]];

   (* Integrals calculation*)
   {integral1, integral2, error1, error2} = 
    IRuleEstimate[baseFunc, depFunc, {a, b}, {absc, weights, errweights}];

   (* Assuming 1D integrals *)
   {{integral1, integral2}, 
    Max[error1, error2], 1}
   ];

Answer 2

First, NIntegrate[f1[x], {x, xmin, xmax}] usually proceeds by constructing an Experimental`NumericalFunction from the expression for f1[x]. This will circumvent an attempt to memoize f1 in the OP's manner, f1[x_] := f1[x] =.... One can prevent this by memoizing the function with ?NumericQ checks via f2[x_?NumericQ] := f2[x] = .... One thing to consider is that the NumericalFunction constructed in each case is different, and if one of them is more efficient, I would bet it's the first one; see below for evidence of this. Anton Antonov alluded to this in a comment.

Second, as Mr.Wizard has pointed out, NIntegrate on a list of integrands just call NIntegrate on each integrand. To expand further, the reason there is no speed-up is that the first integral samples only 116721 points, whereas the second samples 1063425. The time it takes to do the second integral dwarfs the first. Memoizing the first integral isn't going to help much, even if all the saved values are reused. (The reason for the difference is that the Sqrt creates some singularities, as well as complex values, which makes the second integral more difficult to compute.)

Code for testing the sampling:

ClearAll[g];
g[x_, y_, z_] := g[x, y, z] = Exp[Sin[x]] + Cos[y + z];
{pts1, pts2} = 
   Last@Last@Reap@
       NIntegrate[#, {x, 0, 10}, {y, 0, 10}, {z, 0, 10}, 
         Method -> {"LocalAdaptive", "SymbolicProcessing" -> 0}, 
         PrecisionGoal -> 7, 
         EvaluationMonitor :> Sow[{x, y, z}]] & /@
     {g[x, y, z], Sqrt[g[x, y, z]] + x}; // AbsoluteTiming
(*  {10.162, Null}  *)

ClearAll[g];
g[x_, y_, z_] := g[x, y, z] = Exp[Sin[x]] + Cos[y + z];
ptsall = Last@Last@Reap@
      NIntegrate[{Sqrt[g[x, y, z]] + x, g[x, y, z]},
       {x, 0, 10}, {y, 0, 10}, {z, 0, 10}, 
       Method -> {"LocalAdaptive", "SymbolicProcessing" -> 0}, 
       PrecisionGoal -> 7, 
       EvaluationMonitor :> Sow[{x, y, z}]]; // AbsoluteTiming
(*  {10.1304, Null}  *)

Length /@ {pts1, pts2}
Total@%
(*
  {116721, 1063425}
  1180146
*)

Length@ptsall    (* same as the preceeding Total *)
(*  1180146  *)

Code for estimating how much time could be saved by memoizing ~120K sample values, if it worked as expected (about half a second):

ClearAll[g];
g[x_, y_, z_] := g[x, y, z] = Exp[Sin[x]] + Cos[y + z];
Table[g[x, y, z], {x, 0., 10, 1./8}, {y, 0., 10, 1./4}, {z, 0., 10, 1./4}]
  // Flatten // Length // AbsoluteTiming
Table[Sqrt[g[x, y, z]] + x, {x, 0., 10, 1./8}, {y, 0., 10, 1./4}, {z, 0., 10, 1./4}]
  // Flatten // Length // AbsoluteTiming
(*
  {0.81421, 136161}
  {0.243886, 136161}
*)

Code to compare the Experimental`NumericalFunction of the symbolic integrand g and the ?NumericQ version:

ClearAll[g];
g[x_, y_, z_] := g[x, y, z] = Exp[Sin[x]] + Cos[y + z];
NIntegrate[g[x, y, z], {x, 0, 10}, {y, 0, 10}, {z, 0, 10}, 
   Method -> {"LocalAdaptive", "SymbolicProcessing" -> 0}, 
   PrecisionGoal -> 7, 
   IntegrationMonitor :> ((numfnSymbolic = First[#]["NumericalFunction"]) &)];

ClearAll[g];
g[x_?NumericQ, y_?NumericQ, z_?NumericQ] := g[x, y, z] = Exp[Sin[x]] + Cos[y + z];
NIntegrate[g[x, y, z], {x, 0, 10}, {y, 0, 10}, {z, 0, 10}, 
   Method -> {"LocalAdaptive", "SymbolicProcessing" -> 0}, 
   PrecisionGoal -> 7, 
   IntegrationMonitor :> ((numfnNumeric = First[#]["NumericalFunction"]) &)];

Table[numfnSymbolic[x, y, z], {x, 0., 10, 1./8}, {y, 0., 10, 1./8}, {z, 0., 10, 1./8}]
  // Flatten // Length // AbsoluteTiming
Table[numfnNumeric[x, y, z], {x, 0., 10, 1./8}, {y, 0., 10, 1./8}, {z, 0., 10, 1./8}]
  // Flatten // Length // AbsoluteTiming
(*
  {0.747187, 531441}
  {3.70936, 531441}
*)

Note this suggests that the difference between the OP's g and a truly memoized g with ?NumericQ should be around 6 sec. Well, it is:

ClearAll[g];
g[x_?NumericQ, y_?NumericQ, z_?NumericQ] := g[x, y, z] = Exp[Sin[x]] + Cos[y + z];
NIntegrate[{g[x, y, z], Sqrt[g[x, y, z]] + x}, {x, 0, 10}, {y, 0, 10}, {z, 0, 10}, 
   Method -> {"LocalAdaptive", "SymbolicProcessing" -> 0}, 
   PrecisionGoal -> 7]; // AbsoluteTiming
(*  {16.2901, Null}  *)

One caveat: In the Table[] comparison, I tested just one numerical function (from one integration subregion). I cannot say authoritatively that all numerical functions are the same. (Help, anyone?)

Answer 3

I got around to actually evaluating your code and I realized that g is not remembering its values; DownValues[g] only has a length of three. The "solution" is to restrict the function to numeric values, per The difference between "SymbolicProcessing" -> 0 and restricting the function definition to numeric values only, but doing that actually makes the integration much slower; the memoization simply has too large an overhead:

SetAttributes[numArgsQ, HoldFirst]
numArgsQ[_[___?NumericQ]] := True

mem : h[x_, y_, z_]?numArgsQ := mem = Exp[Sin[x]] + Cos[y + z];

NIntegrate[{h[x, y, z], Sqrt[h[x, y, z]] + x}, {x, 0, 10}, {y, 0, 10}, {z, 0, 10}, 
  Method -> {"LocalAdaptive", "SymbolicProcessing" -> 0}, 
  PrecisionGoal -> 7] // Timing

{9.12606, {1429.54, 6081.95 + 59.0571 I}}

(My timing for f or g is only ~ 4.83 seconds.)

Extensive memoization has taken place as indicated by:

DownValues[h] // Length

1108835

For my use of mem and numArgsQ please see:

• Quick way to use conditioned patterns when defining multi-argument function?
• What does the construct f[x_] := f[x] = ... mean?

Answer 4

Testing my comment does indicate that one can cut the number of evaluations of f[] in half easily.

f[x_, y_, z_] := Module[{},
    totCalls++;
    Exp[Sin[x]] + Cos[y + z]
]

totCalls = 0; 
NIntegrate[
    {f[x, y, z], Sqrt[f[x, y, z]] + x}, 
    {x, 0, 10}, {y, 0, 10}, {z, 0, 10}, 
    Method -> {"LocalAdaptive", "SymbolicProcessing" -> 0}, 
    PrecisionGoal -> 7] // RepeatedTiming
totCalls

(* {7.23, {1429.54, 6081.95 + 59.0571 I}} *)
(* 16 *)

totCalls = 0; 
NIntegrate[
    {#, Sqrt[#] + x} &[f[x, y, z]], 
    {x, 0, 10}, {y, 0, 10}, {z, 0, 10}, 
    Method -> {"LocalAdaptive", "SymbolicProcessing" -> 0}, 
    PrecisionGoal -> 7] // RepeatedTiming
totCalls

(* {7.281, {1429.54, 6081.95 + 59.0571 I}} *)
(* 8 *)

This tells me one very important thing: Your slow numerical evaluation is not caused by calling f[] often. Instrumenting f[] a little more turns up some very odd behaviour.

f[x_, y_, z_] := Module[{},
    totCalls++;
    Print[{x, y, z, Exp[Sin[x]] + Cos[y + z]} // RepeatedTiming];
    Exp[Sin[x]] + Cos[y + z]
]

totCalls = 0; 
NIntegrate[
    {f[x, y, z], Sqrt[f[x, y, z]] + x}, 
    {x, 0, 10}, {y, 0, 10}, {z, 0, 10}, 
    Method -> {"LocalAdaptive", "SymbolicProcessing" -> 0}, 
    PrecisionGoal -> 7] // RepeatedTiming
totCalls

(*  {1.9*10^-6 , {10, 10, 10, E^Sin[10]+Cos[20]}}  *)
(*  {1.35*10^-6, {x , y , z , E^Sin[x]+Cos[y+z]}}  *)
(*  ...  essentially repeat those two lines seven more times  ...  *)
(*  {7.34, {1429.54, 6081.95 + 59.0571 I}}  *)
(*  16  *)

So, adding a bunch of overhead to calls to f[] doesn't much change the total calculation time. Also, the function only ever seems to get evaluated at {x,y,z} = {10,10,10}, which is very odd.