# Fast evaluation

I have a particularly nasty expression (many of them),

Full simplify is also not an option.

Ideally, I would like to evaluate this function (preferably fast) for around 1500 different values its argument. The problem is that, if I first construct the numerical list of the arguments and then use

func./list


the evaluation takes too long. Also, alternatively, I've tried to evaluate the function in a table, but also much too long.

EDIT:

Frenorm[Px_, Py_, Pz_] :=
a1s (Px^2 + Py^2) + a3s Pz^2 + a11s (Px^4 + Py^4) + a33s Pz^4 +
a13s (Px^2 Pz^2 + Py^2 Pz^2) +
a12s (Px^2 Py^2) + alpha111 (Px^6 + Py^6 +
Pz^6) + alpha112 (Px^4 (Py^2 + Pz^2) + Pz^4 (Px^2 + Py^2) +
Py^4 (Px^2 + Pz^2)) + alpha123 Px^2 Py^2 Pz^2 + mu

acPhaseSolSet = {Sqrt[p1], Sqrt[p3]} /. Solve[{Expand[D[Frenorm[Sqrt[p1], 0, Sqrt[p3]], p1]] == 0, Expand[D[Frenorm[Sqrt[p1], 0, Sqrt[p3]], p3]] == 0}, {p1, p3}]

acPhaseEnergy =  Table[Frenorm[Px, 0, Pz] /. {Px -> acPhaseSolSet[[j]][[1]], Pz -> acPhaseSolSet[[j]][[2]]}, {j, 1, Length[acPhaseSolSet]}]


Here, acPhaseEnergy, yields a result (there are five solutions) that is so long (above) that won't let you efficiently (and rapidly) input a list value of {a1s, a3s, ...}. I would like to put in 1500 lists of different numbers of a1s, etc and get a grid of energy values.

Ie, for example, maybe these values would be

list = {a1s -> -4.01483*10^8, a3s -> -1.92393*10^8, a11s -> 3.3416*10^8,
a33s -> 4.99091*10^7, a13s -> 8.10182*10^8, alpha111->-4.0*10^7, alpha112->35.0*10^6, alpha123->-3.0*10^7,
a12s -> 1.09218*10^9, mu -> 1.81818*10^7}


so one could try

acPhaseEnergy/.list

• Try Compile. I would have expanded on this if you only had provided copyable code. Commented Apr 17, 2018 at 13:43
• Maybe try Dispatch on list? Commented Apr 17, 2018 at 13:44
• Henrik, the copyable code is much too long. Sorry Commented Apr 17, 2018 at 14:08
• You could try using Block, as in Block[{a1s=1, a11s=2, ..}, func]. Commented Apr 17, 2018 at 15:04
• Not much to go on from this description. Commented Apr 17, 2018 at 16:05

ExperimentalOptimizeExpression gives a significant speed-up. After injecting the values, First will let the optimized expression evaluate:

opt = ExperimentalOptimizeExpression[acPhaseEnergy];
opt /. list // First // AbsoluteTiming
(*
{0.002094, {1.95313*10^9 + 3.446*10^9 I,
1.95313*10^9 - 3.446*10^9 I, -5.38374*10^7, 2.49635*10^12}}
*)

acPhaseEnergy /. list // AbsoluteTiming
(*
{7.05487, {1.95313*10^9 + 3.446*10^9 I,
1.95313*10^9 - 3.446*10^9 I, -5.38374*10^7, 2.49635*10^12}}
*)


To evaluate many times:

It helps if the values are packed arrays. This can be accomplished with DevelopToPackedArray[values]. Many functions create or preserve packed arrays. For instance RandomReal[{0.5, 1.5}, 1500] generates 1500 values in a packed array, so that using ToPackedArray in DeveloperToPackedArray[RandomReal[{0.5, 1.5}, 1500]] is unnecessary.

Block[{n = 1500},
list = {a1s -> -4.01483*10^8 RandomReal[{0.5, 1.5}, n],
a3s -> -1.92393*10^8 RandomReal[{0.5, 1.5}, n],
a11s -> 3.3416*10^8 RandomReal[{0.5, 1.5}, n],
a33s -> 4.99091*10^7 RandomReal[{0.5, 1.5}, n],
a13s -> 8.10182*10^8 RandomReal[{0.5, 1.5}, n],
alpha111 -> -4.0*10^7 RandomReal[{0.5, 1.5}, n],
alpha112 -> 35.0*10^6 RandomReal[{0.5, 1.5}, n],
alpha123 -> -3.0*10^7 RandomReal[{0.5, 1.5}, n],
a12s -> 1.09218*10^9 RandomReal[{0.5, 1.5}, n],
mu -> 1.81818*10^7 RandomReal[{0.5, 1.5}, n]}
];

opt /. list // First // Dimensions // AbsoluteTiming
(*  {0.229654, {4, 1500}}  *)


To get the values in their proper order use Transpose[opt /. list // First].

• ExperimentalOptimizeExpression was definitely the way to go. Commented Apr 18, 2018 at 12:31

As suggested in my comment, using Block instead of ReplaceAll provides a 4-fold speed up:

acPhaseEnergy /. list //AbsoluteTiming

Hold@@list /. Rule->Set /. Hold[s__] :> Block[{s}, acPhaseEnergy] //AbsoluteTiming


{6.36568, {1.95313*10^9 + 3.446*10^9 I, 1.95313*10^9 - 3.446*10^9 I, -5.38374*10^7, 2.49635*10^12}}

{1.56335, {1.95313*10^9 + 3.446*10^9 I, 1.95313*10^9 - 3.446*10^9 I, -5.38374*10^7, 2.49635*10^12}}

Update

Another factor of 6 speedup can be obtained by doing the replacements in stages. For example:

t = Hold@@list /. Rule->Set /. Hold[s__] :> Block[{s}, acPhaseSolSet]; //AbsoluteTiming

Table[Frenorm[p[[1]], 0, p[[2]]] /. list, {p, t}] //AbsoluteTiming


{0.256775, Null}

{0.000228, {1.95313*10^9 + 3.446*10^9 I, 1.95313*10^9 - 3.446*10^9 I, -5.38374*10^7, 2.49635*10^12}}

in agreement with the previous results.

Update 2

We can combine @MichaelE2's nice approach using ExperimentalOptimizeExpression and my Block and 2-step ideas to improve speed further:

opt = ExperimentalOptimizeExpression[acPhaseSolSet]; //RepeatedTiming

t = Hold@@list /. Rule->Set /. Hold[s__] :> Block[{s}, First @ opt]; //RepeatedTiming

Table[Frenorm[p[[1]], 0, p[[2]]] /. list, {p, t}] //RepeatedTiming


{0.182, Null}

{0.00032, Null}

{0.0002, {1.95313*10^9 + 3.446*10^9 I, 1.95313*10^9 - 3.446*10^9 I, -5.38374*10^7, 2.49635*10^12}}

Compare this to not using my Block and 2-step ideas:

opt = ExperimentalOptimizeExpression[acPhaseEnergy]; //RepeatedTiming

opt /. list //First //RepeatedTiming


{1.09, Null}

{0.0018, {1.95313*10^9 + 3.446*10^9 I, 1.95313*10^9 - 3.446*10^9 I, -5.38374*10^7, 2.49635*10^12}}

It takes more time to optimize the acPhaseEnergy expression, and using /. is about 3-4 times slower.

Assuming that

1. the lengthy code is contained in the variable prettylengthyandexpensivecode;

2. you have a list vars containing all symbols in prettylengthyandexpensivecode;

3. all parameters are supposed to be complex numbers; and

4. no special functions are involved in prettylengthyandexpensivecode,

you can try to compile. With the provided data, this could look as follows:

vars = {a1s, a3s, a11s, a33s, a13s, alpha111, alpha112, alpha123,
a12s, mu};
prettylengthyandexpensivecode = acPhaseEnergy
cf = With[{
ccode = N[prettylengthyandexpensivecode],
cvars = Map[var \[Function] {var, _Complex}, vars]
},
Compile[
Evaluate[cvars],
ccode,
CompilationTarget -> "C",
RuntimeAttributes -> {Listable},
Parallelization -> True
]
]


Then, if you have a matrix data of real numbers where the $i$th row contains the bazillion values of vars[[i]] you want to try, just execute

cf@@data


Using the example data from Michael E2, we can process about 1.5 million data sets per second:

Block[{n = 1500000},
list = {
a1s -> -4.01483*10^8 RandomReal[{0.5, 1.5}, n],
a3s -> -1.92393*10^8 RandomReal[{0.5, 1.5}, n],
a11s -> 3.3416*10^8 RandomReal[{0.5, 1.5}, n],
a33s -> 4.99091*10^7 RandomReal[{0.5, 1.5}, n],
a13s -> 8.10182*10^8 RandomReal[{0.5, 1.5}, n],
alpha111 -> -4.0*10^7 RandomReal[{0.5, 1.5}, n],
alpha112 -> 35.0*10^6 RandomReal[{0.5, 1.5}, n],
alpha123 -> -3.0*10^7 RandomReal[{0.5, 1.5}, n],
a12s -> 1.09218*10^9 RandomReal[{0.5, 1.5}, n],
mu -> 1.81818*10^7 RandomReal[{0.5, 1.5}, n]}
];

cf @@ list[[All, 2]]; // AbsoluteTiming // First


1.10661

However, another drawback of this method is that only computations in machine precision are possible within a LibraryFunction. I get numerical errors frequently with this random data. In these cases, the CompiledFunction cf is clever enough to switch to a standard Mathematica function as fallback (so the results should be accurate but not as fast as with a LibraryFunction).

By the way: As far as I know, ExperimentalOptimizeExpression is also used internally by Compile`.

• I've edited my question to include minimal working code Commented Apr 17, 2018 at 18:25