8

Given that the documentation for FindMinimum uses ReplaceAll for its various examples of what to do with the found solution to an optimization problem, it would seem to be the right tool for the job. However, the following MWE highlights the speed problem I ran into using ReplaceAll for such a situation. Is using Block the appropriate tool in this case and why?

vars = Table[f[i], {i, 1, 10000}];
vals = RandomReal[1, Length[vars]];
soln = MapThread[#1 -> #2 &, {vars, vals}];
Total[vars] /. soln // Timing
(*{1.395787, 4980.51}*)

Block[{f},
 ReleaseHold[
  MapThread[Hold[#1 = #2] &, Transpose[List @@ # & /@ soln]]]; 
  Total[vars]] // Timing
(*{0.033995, 4980.51}*)

A slight variation on the Block approach, which I expect would be equivalent, is also slightly slower (here both have been repeated 1000 times to show the less-significant difference of ~20%). Is there another approach that would achieve the same but faster?

Do[Block[{f}, 
   ReleaseHold[
    MapThread[Hold[#1 = #2] &, Transpose[List @@ # & /@ soln]]]; 
   Total[vars]], {1000}]; // Timing
(*{26.072036, Null}*)

Do[Block[{f}, (Set @@ # &) /@ (List @@ # & /@ soln); 
    Total[vars]], {1000}]; // Timing
(*{31.651188, Null}*)

Further testing has found that the problem is exasperated by "bigger" expressions than Total, in particular by the use of ListInterpolation. This appears to be resolved by only Replacing at the right level, which works in this particular example because all the variables have the same levelspec ({-2} or {3}).

vars = Table[f[i], {i, 1, 10000}];
vals = RandomReal[{-1, 1}, Length[vars]];
soln = MapThread[#1 -> #2 &, {vars, vals}];
i = ListInterpolation[RandomReal[1, 40], {-1, 1}, 
   "ExtrapolationHandler" -> {(0 &), "WarningMessage" -> False}];
expr = Plus @@ 
   MapThread[Times, RotateLeft[i@vars, #] & /@ Range[0, 4]];

expr /. Dispatch[soln] // Timing
(*{2.504619, 309.407}*)

Replace[expr, Dispatch[soln], {-2}] // Timing
(*{0.124981, 309.407}*)

Block[{f}, 
  ReleaseHold[
   MapThread[Hold[#1 = #2] &, Transpose[List @@ # & /@ soln]]]; 
  expr] // Timing
(*{0.146978, 309.407}*)

In the actually problem I am facing, not all the variables have the same levelspec, so it seems to require the range {-2,-1}, which ends up taking slightly longer than Block.

Joel Bosveld
  • 257
  • 1
  • 10

1 Answers1

13

Original question

As noted in the comments the use of Dispatch is the easiest way to make this replacement operation much faster. However taking this as an opportunity to explore other optimizations here are some examples for you to consider:

(dsp = Dispatch[soln]) // RepeatedTiming // First
0.0030
Total[vars] /. dsp // RepeatedTiming
{0.0038, 5025.19}

It is faster if we focus the replacement on one level alone:

Replace[Total[vars], dsp, {-2}] // RepeatedTiming
{0.0032, 5025.19}

Interestingly levelspec {1} appears slightly faster than {-2}:

Replace[Total[vars], dsp, {1}] // RepeatedTiming
{0.0030, 5025.19}

Moving the Total outside prevents a wasted summation attempt:

Total @ Replace[vars, dsp, {1}] // RepeatedTiming
{0.0019, 5025.19}

Lookup on an Association is faster however:

asc = Association[soln];

Total @ Lookup[asc, vars] // RepeatedTiming
{0.0017, 5025.19}

Finally, Tr is faster than Total on packed arrays, and packing has low overhead here:

Tr @ Developer`ToPackedArray @ Lookup[asc, vars] // RepeatedTiming
{0.00088, 5025.19}

Updated question

In your update expr is a large expression containing many InterpolatingFunction subexpressions. Each of these is significant:

expr[[1, 1]] // LeafCount    (* 153 *)

Further they each contain multiple packed arrays:

expr[[1, 1, 0]] /. x_List?Developer`PackedArrayQ :> Print[Short@x];
{40}

{4}

{-1.,-0.948718,<<37>>,1.}

{0,1,2,3,4,5,6,7,8,<<23>>,32,33,34,35,36,37,38,39,40}

{0.947691,<<38>>,0.386701}

All of this is wasted space to search for a pattern match and every packed array must be first unpacked incurring additional overhead. Since you only want to replace expressions in the form f[_] and since these all occur in the same place in each subexpression it will be far better to target the replacement accordingly. That is exactly what you did with Replace[expr, Dispatch[soln], {-2}] which is faster than Block in this case. (Indicientally Block avoids this problem because the evaluator will have flagged the InterpolatingFunction expressions as already evaluated and not depending on f and therefore not requiring updating, as part of its optimizations.)

Since the long InterpolatingFunction expressions appear as heads you could also use Heads -> False with almost the same efficiency even with Infinity as the depth:

Replace[expr, Dispatch[soln], {-2}] // Timing
Replace[expr, Dispatch[soln], ∞, Heads -> False] // Timing
{0.121, 319.564}

{0.145, 319.564}

But at least for this example it brings up another question: why aren't you doing the replacement on vars first and then building expr?

With[{ivars = i @ Replace[vars, Dispatch[soln], {1}]},
  Times @@ NestList[RotateLeft, ivars, 4] // Total
] // Timing // First
0.0237

By the way if all your replacements are of the form shown there is a simpler way to use Block:

Block[{f},
  DownValues[f] = soln;
  expr
] // Timing
{0.143, 319.564}
Mr.Wizard
  • 271,378
  • 34
  • 587
  • 1,371
  • Dispatch does improve the speed significantly, it wasn't as good as Block but once I had an Evaluate in the right spot it was just as fast (the evaluation order seems to be my biggest hurdle).

    I had forgotten about Dispatch as last time I used it, I didn't notice any improvement; something else must have been the bottleneck in that situation (or, I remember reading somewhere that in some situations it is used automatically?).

    I can't test the Lookup method at the moment, but may be worth going through the effort of upgrading to v10 to get the further improvement.

    – Joel Bosveld May 20 '15 at 04:50
  • Actually, for the case I am particularly interested in, Dispatch is still 20x slower than Block. I'll have to give it some more testing it appears... – Joel Bosveld May 20 '15 at 05:16
  • @Joel I more representative example would be helpful in that case. – Mr.Wizard May 20 '15 at 05:53
  • I have found that a more complicated expression causes trouble (although usually not more than 2 or 3 times). What really seems to aggravate it is when I put in a ListInterpolation; then Block is 10 or 20 times faster than Dispatch (depending on how big the expression is). Is that a small enough change to allow an edit/addition, or do I need to ask a new question? – Joel Bosveld May 20 '15 at 08:05
  • @Joel Thanks for asking. Mine is the only answer so I am the only one affected by a "moving goalposts" question, and I encourage you to make the edit. – Mr.Wizard May 20 '15 at 08:07
  • Done. Using the right levelspec in Replaces does seem to do the trick, although in the actual problem I am facing, it still takes slightly more time. – Joel Bosveld May 20 '15 at 08:42
  • @Joel I updated my answer. – Mr.Wizard May 20 '15 at 09:28
  • Thanks very much. For the actually problem, the Heads->False and {-2,-1} methods are quite similar in time (although still a bit slower than Block) and much faster than ReplaceAll/Dispatch. As for the last question, the building of the expression is a lot more involved than shown in the MWE, though it may be worth trying. Also expr gets passed to FindMinimum (to optimize over vars), so the work of building it has been done and I just want to evaluate it (e.g. at the solution, initial condition). – Joel Bosveld May 20 '15 at 10:10