3
$\begingroup$

Note: this is a clarification of 2 earlier unanswered questions, Compile issues, scoping and order of evaluation and https://mathematica.stackexchange.com/questions/60549/optimising-a-simple-computation-with-vectors .

The problem

I have a list of expressions, e.g.

expressions = {a b + x^7 a^2 + x^30 Ma^2, a b c - 2 Mb^2 x + Mb Mc, a b^2 - Mb Md - x^10 d^2};

where x is a known fixed vector, M is a known fixed matrix, Ma stands for M.a, etc. and the rest are unknown vectors.

x = {1, 2};
M = {{0, 1}, {1, 0}};

The first is a function of a,b , the second of a,b,c and the last of a,b,d.

Note: this is a minimal example, in the actual problem each of these expressions is a list of expressions like it (of fixed length), they are longer and there's more of them, and they need to be evaluated thousands to millions of times.

A code to evaluate these quickly is something like

eval1 = Compile[{{a, _Real, 1}, {b, _Real, 1}},
   Module[{Ma = M.a}, a b + x^7 a^2 + x^30 Ma^2], CompilationTarget -> "C"];

eval2 = Compile[{{a, _Real, 1}, {b, _Real, 1}, {c, _Real, 1}},
Module[{Ma = M.a, Mb = M.b, Mc = M.c}, a b c - 2 Mb^2 x + Mb Mc], 
  CompilationTarget -> "C"];

What I want is to automate the writing of these functions, i.e. I want to write something like

Do[eval[i] = Compile[variables[[i]], Module[Mvariables[[i]], expressions[[i]]], 
   CompilationTarget -> "C"],{i, 1, Length[expressions]}]

where variables is the argument list, Mvariables is the assignment of M.a to Ma etc.

It is also important that all that can be evaluated in expression beforehand, i.e. without knowing the arguments, is evaluated once and for all and not every time. (for example the x^30 shouldn't be computed every time the function is evaluate but only when defining the function)

Partial solution

variables = Union@Cases[#, x_ /; Not[FreeQ[{a, b, c, d}, x]], Infinity] & /@ expressions
Do[exprs[i] = expressions[[i]], {i, 1, Length[expressions]}]

loopCode[i_] := Compile[Evaluate@({#, _Real, 1} & /@ variables[[i]]),
   Module[{Ma = M.a, expr = exprs[i]}, expr], CompilationTarget -> "C",
    CompilationOptions -> "InlineExternalDefinitions" -> True]

Do[finalCode[i] = loopCode[i], {i, 1, Length[expressions]}]

This gives,

finalCode[1][{1, 2}, {2, 2}]
{{7., 9.}, {36., 20.}}

This solution has 2 problems.

  1. I cheated with the initialisation in module.

I want it to figure out for itself which ones are needed. Doing something like

Mvariables = Union@Cases[#, x_ /; Not[FreeQ[{Ma, Mb, Mc, Md}, x]], Infinity] &/@expressions;
Mvariables2 = Map[ToExpression[StringTake[ToString[#], -1]] &, Mvariables, {2}];

Module[Evaluate@(Mvariables[[i]]=(M*#&)/@variables[[i]]),...]

doesn't work as it still isn't evaluated.

  1. The expression (exprs[i]) isn't evaluated as much as possible.

To see this we can check:

M = IdentityMatrix[30];
x= 1/Range[30]//N;
c1=x^7;
c2=x^30;
newCode=Compile[{{a,_Real,1},{b,_Real,1}},
         Module[{Ma=M.a}, a b + c1 a^2 + c2 Ma^2, CompilationTarget -> "C",
    CompilationOptions -> "InlineExternalDefinitions" -> True];
atest = RandomReal[{0, 1}, 30];
btest = RandomReal[{0, 1}, 30];

Now timing,

Do[finalCode[1][atest, btest], {10^5}] // Timing
{0.565053, Null}
Do[newCode[atest, btest], {10^5}] // Timing
{0.365492, Null}

I think both issues have to do with evaluation inside module, but I don't see how to fix it.

$\endgroup$
4
$\begingroup$

I suggest that you use local rules for your fixed values. Then, the following seems to do what you ask for:

ClearAll[dotQ, init, vars];
dotQ[s_String, MSymbol_Symbol] := 
   StringMatchQ[s, ToString[MSymbol] ~~ LetterCharacter];

vars[MSymbol_] := Select[Not@*(dotQ[#, MSymbol] &)@*ToString]@*Variables;

init[expr_, MSymbol_Symbol] :=
   Thread[#, Hold] & @ Union @ Cases[
      expr,
      s_Symbol :>
         With[{eval = ToString[s]},
           ToExpression[
             StringDrop[eval, 1],
             StandardForm,
             Function[sym, Hold[s = MSymbol.sym], HoldAll]
           ] /; dotQ[eval, MSymbol]
         ],
      Infinity
   ];

ClearAll[compile];
compile[expr_, MSymbol_, rules_] :=
  With[{
     i = init[expr, MSymbol] /. rules,
     computed =
       Hold[
         Evaluate[expr /. {Plus -> Hold[Plus], Times -> Hold[Times]} /. rules]
       ] /. {Hold[Plus] -> Plus, Hold[Times] -> Times}
     },
     Prepend[
       i /. Hold[x_] :> (computed /. Hold[y_] :> Hold[Hold[Module][x, y]]) /. Hold[Module] -> Module,
       {#, _Real, 1} & /@ vars[MSymbol][expr /. rules]
     ] /. Hold[args__] :> Compile[args, "CompilationTarget" -> "MVM"] 
 ];

You use this as:

ClearAll[M,a,b,c,d,x];
eval = Map[compile[#, M, {x -> {1, 2}, M -> {{0, 1}, {1, 0}}}] &]@expressions

Some examples to illustrate what the functions do:

vars[M]@expressions[[1]]

(* {a, b, x} *)

init[expressions[[2]], M]

(* Hold[{Mb = M.b, Mc = M.c}] *)

This could've been done in a less cryptic manner without using the tricks I use, but that would require some more code.

Note that, to my mind, your choice of using symbols like Ma etc needlessly complicate things. Note also that I am using operator forms in a few places, so this would have to be modified to work on V9 and earlier.

$\endgroup$
  • $\begingroup$ Thanks, I have to study this to see how it works, but for me it doesn't give the result I want. When i apply the first function to the arguments {1,2},{3,4} I get a 2x2 matrix, whereas it should be a vector. Maybe I should have mentioned this: the multiplication in the expressions should be interpreted as element wise multiplication of lists. $\endgroup$ – Jansen Sep 25 '14 at 16:26
  • $\begingroup$ And the reason that I'm using these Ma variables is that in the actual expressions they occur multiple times, so this limits the number of matrix multiplications to the minimum. $\endgroup$ – Jansen Sep 25 '14 at 16:37
  • $\begingroup$ @Jansen Ok, the problem there was that on one hand, you want it to be partially evaluated, on another hand, Plus and Times are Listable, so one has to fight that during the partial computation of your expression. I didn't have the time to make the code pretty, but that part has now been corrected. At this point, it becomes clear that your problem is complex enough that a custom expression parser / code generator would be much cleaner than what I wrote. No time for this at the moment though. But I do realize that the code above can be hard to digest. $\endgroup$ – Leonid Shifrin Sep 25 '14 at 16:53
  • $\begingroup$ Thanks, it works exactly as I wanted. Don't know what you mean by custom expression parser/code generator though. Btw, also thanks for writing the book on mathematica, learned a lot from it :) $\endgroup$ – Jansen Sep 26 '14 at 8:30
  • $\begingroup$ @Jansen Glad I could help, thanks for the accept. Re: book - good to know that it is useful! $\endgroup$ – Leonid Shifrin Oct 2 '14 at 9:07

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.