How to ensure efficient order of operations for rules, slots, substitutions?

Question

I need to calculate several functionals Q[[t]] of a pdf f[x].

maxt = 5;
integrands = Table[Expand[f[x]*D[Log[f[x]], {x, 2*t}]/(2^t)], {t, 1, maxt, 1}];
fSub1[expr_] := {f[x] -> expr // N, Derivative[i_][f][x] -> D[expr, {x, i}] // N};

Later, the program will receive candidate pdfs c[x] for which to calculate values of the Q[[t]]'s:

gauMix[means_, vars_] := Total[Apply[(1/(Sqrt[2*Pi*#2]*Length[vars]))*E^-(((x - #1)^2)/(2*#2)) &, Transpose[{means, vars}], 1]];
c[x] = gauMix[{1, 2, 3, 4, 5}, {1, 2, 4, 2, 1}];
Q = Table[0, {t, 1, maxt, 1}];
Table[Q[[t]] = NIntegrate[integrands[[t]] /. fSub1[c[x]], {x, -Infinity, Infinity}, Method -> {Automatic, "SymbolicProcessing" -> False}], {t, 1, maxt, 1}] // AbsoluteTiming

The "//N" in fSub1 and the "SymbolicProcessing->False" in the NIntegrate generate the fastest results I've achieved for this problem.

But I'm wondering if the substitutions are being done in an optimally efficient manner. The f[x] and each of its derivs appears multiple times in every integrand, for example integrands[[2]] =

$$\frac{1}{4} f^{(4)}(x)-\frac{3 f''(x)^2}{4 f(x)}-\frac{3 f'(x)^4}{2f(x)^3}-\frac{f^{(3)}(x) f'(x)}{f(x)}+\frac{3 f'(x)^2 f''(x)}{f(x)^2} $$

So the order of operations I HOPE is happening is:

1) The equation forms of derivatives of c[x] are cacluated once each.

2) The value of each derivative of c[x] is calculatece once each for every x.

3) These values are substituted by fSub1 in for the derivs of f[x] in the integrands.

4) All the integrands and integrated simultaneously so f[x] and derivs won't have to be recalculated over the domain of x for every t.

But I have no idea how to test if this is happening or not. If it's happening, how can I know? It it's not, how can I force it? Is this the most efficient way to get this done?

I attempted to force the order of ops to work this way but using a different substitution pattern, but this substitution did not work, and I have no idea why.

dfdx = Table[0, {t, 1, maxt, 1}];
Do[dfdx[[i]] = D[c[x], {x, i}], {i, 1, maxt, 1}];
fSub2 := {f[x] -> c[x], Derivative[j_][f][x] -> dfdx[[j]]};
Table[Q[[t]] = NIntegrate[integrands[[t]] /. fSub2, {x, -Infinity, Infinity}, Method -> {Automatic, "SymbolicProcessing" -> False}], {t, 1, maxt, 1}] // AbsoluteTiming

Thanks!

BTW the methods I used to speed up the integration as much as I did came from THIS very helpful discussion: How to use slots and patterns to reduce repetitive calculation?

Answer 1

Answer the last question first. In a word, it's a simple mistake:

Cases[integrands[[-1]], Derivative[_], Infinity, Heads -> True] // Sort // Last

Derivative[10]

As you see, the highest order of derivative in integrands isn't tmax which equals 3.

Another small mistake that leads to some warnings but doesn't influence the result is that Derivative[j_][f][x] -> dfdx[[j]] should be Derivative[j_][f][x] :> dfdx[[j]], which is just a matter of calculation order.

Fixed code:

maxt = 5;
integrands = Table[f[x] D[Log@f[x], {x, 2 t}]/(2^t), {t, maxt}];

dfdx = With[{maxd = 10}, Table[D[c[x], {x, i}], {i, maxd}]];
fSub2 = {f[x] -> c[x], Derivative[j_][f][x] :> dfdx[[j]]};
Q = NIntegrate[integrands /. fSub2, {x, -Infinity, Infinity}, 
    Method -> {Automatic, "SymbolicProcessing" -> False}] // AbsoluteTiming

{0.343200, {-0.166669, 0.00805105, -0.0068687, 0.00741254, -0.0116914}}

Then let me deal with your 4 HOPEs. The first 3 of them are all about substitution, let alone whether substitution will influence the performance or not, actually you can omit all the substitution in your code if you like:

f[x_] = gauMix[{1, 2, 3, 4, 5}, {1, 2, 4, 2, 1}];
Q = NIntegrate[integrands, {x, -Infinity, Infinity}, 
   Method -> {Automatic, "SymbolicProcessing" -> False}] // AbsoluteTiming

{0.327600, {-0.166669, 0.00805105, -0.0068687, 0.00741254, -0.0116914}}

In this case the performance is almost the same, with some negligible oscillation. Maybe the difference of performance will be significant when maxt is larger. (Anyway, ReplaceAll is actually slow for large expressions.) But I think a large maxt will lead to some precision issue that one will have to turn to high precision calculation (WorkingPrecision -> (*some integers*) etc.) that will significantly slow down the integration, to guarantee the correctness of the result so the time consuming of substitution will no longer be a issue.

Last, just talk a little about the 4th HOPE. I believe, as one of the most complicated function in Mathematica, NIntegrate has done almost all that one (at least one belonging to average users) can imagine to optimize the integration. (As you see in this and the former post, the form of your expression rarely influence the timing of NIntegrate.) Maybe later someone more knowledgeable will come with some magic more powerful, but this is my conclusion.