# How to speed up calculations with large number of replacement rules?

In my program I have an expression which looks like:

AnalyticalExpression = 2^(-1-n-np)R^(-2+n+np)(8np1f[-2+n+np,-1+p1+p1p,1+p2+p2p,1+q1+q1p,3+q2+q2p]-8np1f[-2+n+np,-1+p1+p1p,3+p2p2p,1+q1+q1p,1+q2+q2p]-4np1f[-2+n+np,p1+p1p,p2+p2p,q1+q1p,4+q2+q2p]+...)


I gave here a short example; full expression is few pages long. From my point of view, this expression serves as a function of the integer variables n,np,p1,p1p,q1,q1p,... and two real variables R,x. My program evaluates the above expression by using a replacement rule, for instance:

tmp = AnalyticalExpression /. {n->2,np->3,p1->0,p2->4,...}


The problem is that the program has to evaluate this expression hundreds of thousands times with different values of these integer variables. This makes the whole program painfully slow. Is there a tricky way to perform this substitution efficiently? Tabulation of AnalyticalExpression for every combination of the integer variables seems quite weird to me...

Michal

-
It is still not exactly clear what you do and whether the slow down comes from the large expression or the substitution. Are you changing the values in the list of replacement rules? Have you seen the function Dispatch? – halirutan Sep 4 '13 at 0:16
Why don't you make it a function with named arguments? – Sjoerd C. de Vries Sep 4 '13 at 6:05
I can. Would it change the cost of the evaluation? – michal Sep 5 '13 at 0:00

As halirutan comments Dispatch will speed the application of long lists of rules:

SetAttributes[timeAvg, HoldFirst]
timeAvg[func_] := Do[If[# > 0.3, Return[#/5^i]] & @@ Timing@Do[func, {5^i}], {i, 0, 15}]

n = 1500;
big = Sum[Expand[(RandomInteger[99] + a[i])^RandomInteger[9]], {i, n}];

vals = RandomInteger[9999, n];
rules = Thread[Array[a, n] -> vals];

big /. rules           // timeAvg
big /. Dispatch[rules] // timeAvg


0.936

0.008608

Somewhat faster, you could convert your entire expression into an anonymous Function one time, and Apply it to your list of values:

big2 = Function @@ {big} /. Dispatch@Thread[Array[a, n] -> Array[Slot, n]];

big2 @@ vals // timeAvg


0.004496

For variety, if all of the expressions to replace are Symbols you could use a variation of listWith from this post. It is between the two others in speed in this test.

syms = Table[Symbol["a" <> ToString@i], {i, n}];

big3 = big /. Dispatch@Thread[Array[a, n] -> syms];

SetAttributes[listBlock, HoldAll];
listBlock[(set : Set | SetDelayed)[L_, R_], body_] :=
Inner[set, L, R, Hold] /. _[x__] :> Block[{x}, body]

listBlock[syms = vals, big3] // timeAvg


0.005368

This has the advantage of keeping your expression in the normal form rather than turning it into a Function.

## Optimizing the expression

In light of the nature of the actual expression which contains a number of repeated sub-expressions we can improve the performance a bit by exploiting this commonality. An internal function exists that attempts this automatically: ExperimentalOptimizeExpression. You can search the site for "OptimizeExpression" for more details and other examples of use.

With your full expression assigned to fullAE:

vars = Variables@Level[fullAE, {-1}];
rules = MapIndexed[Function[{a, b}, a -> Slot[b[[1]]]], vars];

func = Function @@ {fullAE} /. rules;

optimized = ExperimentalOptimizeExpression[fullAE];
funcOptimized = Function @@ optimized /. rules;

values = RandomInteger[99, 12];


Timings for the incrementally faster methods:

fullAE /. Thread[vars -> values]          // timeAvg
fullAE /. Dispatch@Thread[vars -> values] // timeAvg
func @@ values                            // timeAvg
funcOptimized @@ values                   // timeAvg


0.00324

0.0022208

0.001048

0.0006752

-
Thank you for your answers, halirutan and Mr. Wizard. Unfortunately, Dispatch[] does not improve the timings significantly because I change the integer values in the replacement rule at every evaluation. At the very beginning of the evaluation the replacement rule looks like {n->0,np->0,p1->0,p2->0,...} and the last may be {n->10,np->10,p1->10,p2->10,...}. Therefore, I have to generate a dispatch table as many times as the replacement rule is applied. This is clearly not the same situation as in the halirutan's example. – michal Sep 4 '13 at 0:49
@michal Actually, I specifically included the creation of the Dispatch table within the timing loop to address that situation. (Of course I may have made a mistake). How long are your replacement rule lists? I chose 1500; if yours are short Dispatch may not help much. What is the LeafCount of your expression? – Mr.Wizard Sep 4 '13 at 0:55
My replacement rule is about 20 entries long. Now I realised that yours is 1500 entries long. LeafCount[] of my expression is 11562. – michal Sep 4 '13 at 0:57
@michal Could you put the entire expression on pastebin.com ? If it can't be compiled there may not be a lot more than can be done beyond the options I have given. – Mr.Wizard Sep 4 '13 at 1:05
@michal I updated my answer with an automatic method. There may be further optimizations to be had, but I need to move on to other things. I'll try to return to this problem but I'll admit up front I have a poor record when it comes to returning to problems. – Mr.Wizard Sep 4 '13 at 1:40