1
$\begingroup$

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

enter image description here

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
$\endgroup$
  • $\begingroup$ Try Compile. I would have expanded on this if you only had provided copyable code. $\endgroup$ – Henrik Schumacher Apr 17 '18 at 13:43
  • $\begingroup$ Maybe try Dispatch on list? $\endgroup$ – Chris K Apr 17 '18 at 13:44
  • $\begingroup$ Henrik, the copyable code is much too long. Sorry $\endgroup$ – John M Apr 17 '18 at 14:08
  • $\begingroup$ You could try using Block, as in Block[{a1s=1, a11s=2, ..}, func]. $\endgroup$ – Carl Woll Apr 17 '18 at 15:04
  • 1
    $\begingroup$ Not much to go on from this description. $\endgroup$ – Daniel Lichtblau Apr 17 '18 at 16:05
4
$\begingroup$

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].

$\endgroup$
  • $\begingroup$ Experimental`OptimizeExpression was definitely the way to go. $\endgroup$ – John M Apr 18 '18 at 12:31
3
$\begingroup$

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.

$\endgroup$
3
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ I've edited my question to include minimal working code $\endgroup$ – John M Apr 17 '18 at 18:25

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.