Fast evaluation

Question

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:

More information:

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

Answer 1

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

opt = Experimental`OptimizeExpression[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 Develop`ToPackedArray[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 Developer`ToPackedArray[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].

Answer 2

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 Experimental`OptimizeExpression and my Block and 2-step ideas to improve speed further:

opt = Experimental`OptimizeExpression[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 = Experimental`OptimizeExpression[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.

Answer 3

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, Experimental`OptimizeExpression is also used internally by Compile.