Optimize evaluation of symbolic expressions

[Harrie Kraai]

Hello all,

With a specific example I countered today I want to introduce a question
that I have had for a longer time.

From some composition of functions (and a subsequent FullSimplify) I
get the following expression:

{ArcTan[-1 + 2/(1 + x^2 + y^2), 2 Sqrt[(x^2 + y^2)/(1 + x^2 + y^2)^2]],
ArcTan[x/(1 + x^2 + y^2), y/(1 + x^2 + y^2)]}

which I need to evaluate many times inside a Manipulate.

Now, how does Mathematica evaluate this?

If I would program this myself, I would, of course, evaluate the
subexpression (1 + x^2 + y^2) first and use this result 4 times in the
subsequent calculations.

It looks to me like Mathematica is not able to or at least does not make
this optimization. (Find out using 'Trace'). Perhaps this is difficult
to do in a general sense because any subexpression may have
side-effects. However, in this case it should be possible to instruct
Mathematica to collect common subexpressions first. It would surely make
calculations much faster.
The symbolic calculations have led to an expression that is
"Simplify"-ed in terms of reading (perhaps) but not in terms of
evaluation. Is there, or should there be a function that translates
expressions to a form that is optimized for evaluation?

This appears to me as a frequently encountered problem. How do we go
about it?

Thanks for your advice and/or thoughts.

Harrie

[Roland Franzius]

Harrie Kraai schrieb:

> Hello all,
>
> With a specific example I countered today I want to introduce a question
> that I have had for a longer time.
>
> From some composition of functions (and a subsequent FullSimplify) I
> get the following expression:
>
> {ArcTan[-1 + 2/(1 + x^2 + y^2), 2 Sqrt[(x^2 + y^2)/(1 + x^2 + y^2)^2]],
> ArcTan[x/(1 + x^2 + y^2), y/(1 + x^2 + y^2)]}
>
> which I need to evaluate many times inside a Manipulate.

Define something like

f[x_,y_]:=
({ArcTan[-1+2/#, 2 Sqrt[(#-1)/#^2], ArcTan[x/#,y/#]}&)[1+x^2+y^2]

but for an interpreting language the gain will be moderate if more than
this single function only has to be evaluated during an inner loop.

--

Roland Franzius

[Bob Hanlon]

You need to help it

expr = {ArcTan[-1 + 2/(1 + x^2 + y^2),

2 Sqrt[(x^2 + y^2)/(1 + x^2 + y^2)^2]],

ArcTan[x/(1 + x^2 + y^2), y/(1 + x^2 + y^2)]};

Note that if x and y are Real, this can be further simplified

Simplify[expr, Element[{x, y}, Reals]]

{ArcTan[-x^2 - y^2 + 1, 2*Sqrt[x^2 + y^2]], ArcTan[x, y]}

In the general case,

Simplify[expr, t == x^2 + y^2 + 1]

{ArcTan[2/t - 1, 2*Sqrt[(t - 1)/t^2]], ArcTan[x/t, y/t]}

Manipulate[Module[{t = x^2 + y^2 + 1},
{ArcTan[2/t - 1, 2*Sqrt[(t - 1)/t^2]],
ArcTan[x/t, y/t]}],
{x, -1, 1, .01, Appearance -> "Labeled"},
{y, -1, 1, .01, Appearance -> "Labeled"}]

Bob Hanlon

---- Harrie Kraai <hak...@xs4all.nl> wrote:

=============
Hello all,

With a specific example I countered today I want to introduce a question
that I have had for a longer time.

From some composition of functions (and a subsequent FullSimplify) I
get the following expression:

{ArcTan[-1 + 2/(1 + x^2 + y^2), 2 Sqrt[(x^2 + y^2)/(1 + x^2 + y^2)^2]],
ArcTan[x/(1 + x^2 + y^2), y/(1 + x^2 + y^2)]}

which I need to evaluate many times inside a Manipulate.

Now, how does Mathematica evaluate this?

[DrMajorBob]

expr1 = {ArcTan[-1 + 2/(1 + x^2 + y^2),

2 Sqrt[(x^2 + y^2)/(1 + x^2 + y^2)^2]],
ArcTan[x/(1 + x^2 + y^2), y/(1 + x^2 + y^2)]};

expr2 = expr1 /. (1 + x^2 + y^2 -> s)

{ArcTan[-1 + 2/s, 2 Sqrt[(x^2 + y^2)/s^2]], ArcTan[x/s, y/s]}

f[x_, y_] = Block[{s = 1 + x^2 + y^2}, expr2];
f[a, b]

{ArcTan[-1 + 2/(1 + a^2 + b^2),
2 Sqrt[(a^2 + b^2)/(1 + a^2 + b^2)^2]],
ArcTan[a/(1 + a^2 + b^2), b/(1 + a^2 + b^2)]}

Bobby

--
DrMaj...@bigfoot.com

[Andrzej Kozlowski]

This has been frequently asked and answered, usually (or maybe always)
by Daniel Lichtblau. His latest post on this topic was, I think,

http://groups.google.com/group/comp.soft-sys.math.mathematica/msg/10fed4b8bd60b762

Here is a brief supplementary comment. The function
OptimizeExpression, which is in the Experimantal context (so it can be
modified in the future or even disappear altogether) will attempt to
do what you wish. However, it is really intended for efficient
compiling of numerical functions and may not be convenient to use for
other purposes.

Here is how it works in the case of your example:

Define a compiled function by:

g = Compile[{x, y},
Evaluate[
Experimental`OptimizeExpression[{ArcTan[-1 + 2/(1 + x^2 + y^2),

2 Sqrt[(x^2 + y^2)/(1 + x^2 + y^2)^2]],

ArcTan[x/(1 + x^2 + y^2), y/(1 + x^2 + y^2)]}, OptimizationLevel
-> 2]]]

evaluate the above and look at the output, you will see how the
expression became optimized. Now you can use it for efficient
computation:

g[1, 2]
{2.30052, 1.10715}

Andrzej Kozlowski

[Stonewall Ballard]

On Jul 25, 4:18 am, Harrie Kraai <hakr...@xs4all.nl> wrote:

> ...

> If I would program this myself, I would, of course, evaluate the
> subexpression (1 + x^2 + y^2) first and use this result 4 times in the
> subsequent calculations.
>
> It looks to me like Mathematica is not able to or at least does not make
> this optimization. (Find out using 'Trace'). Perhaps this is difficult
> to do in a general sense because any subexpression may have
> side-effects. However, in this case it should be possible to instruct
> Mathematica to collect common subexpressions first. It would surely make
> calculations much faster.
> The symbolic calculations have led to an expression that is
> "Simplify"-ed in terms of reading (perhaps) but not in terms of
> evaluation. Is there, or should there be a function that translates
> expressions to a form that is optimized for evaluation?

>...

I have read that internally, Mathematica recognizes common
subexpressions and evaluates them only once. That doesn't help me when
I use Mathematica to create C code, so I wrote an explicit common
subexpression hoister that rewrites an expression into one with
variables holding multiply-used parts. It's not perfect, but it works
for me.
<http://stoney.sb.org/wordpress/2009/06/converting-symbolic-
mathematica-expressions-to-c-code/>

Maybe this would be useful to you too.

- Stoney

[biost...@gmail.com]

Hi, I'm having a similar problem... I have generated a set of expressions in Mathematica, in a format like:
(a X1 + b X2 + c X3 + d) / (e X1 + f X2 + g X3 + h) + (i X1 + j X2 + k X3 + l)(m X1 + n X2 + o X3 + p)/(q X1 + r X2 + s X3 + t)^2 ,
here is an example:
>(-C1 (-k13 Cos[r1] Cos[r2] Cos[r3] - k11 Cos[r2] Sin[r1]) +
> X1 (-k13 Cos[r1] Cos[r2] Cos[r3] - k11 Cos[r2] Sin[r1]) -
> C3 (k11 Cos[r1] - k13 Cos[r3] Sin[r1]) +
> X3 (k11 Cos[r1] - k13 Cos[r3] Sin[r1]) -
> C2 (-k13 Cos[r1] Cos[r3] Sin[r2] - k11 Sin[r1] Sin[r2]) +
> X2 (-k13 Cos[r1] Cos[r3] Sin[r2] - k11 Sin[r1] Sin[r2]))/(-C3 Cos[
> r1] Cos[r3] + X3 Cos[r1] Cos[r3] -
> C2 (-Cos[r3] Sin[r1] Sin[r2] - Cos[r2] Sin[r3]) +
> X2 (-Cos[r3] Sin[r1] Sin[r2] - Cos[r2] Sin[r3]) -
> C1 (-Cos[r2] Cos[r3] Sin[r1] + Sin[r2] Sin[r3]) +
> X1 (-Cos[r2] Cos[r3] Sin[r1] +
> Sin[r2] Sin[r3])) - ((C1 Cos[r1] Cos[r2] Cos[r3] -
> X1 Cos[r1] Cos[r2] Cos[r3] + C3 Cos[r3] Sin[r1] -
> X3 Cos[r3] Sin[r1] + C2 Cos[r1] Cos[r3] Sin[r2] -
> X2 Cos[r1] Cos[r3] Sin[
> r2]) (-C3 (k13 Cos[r1] Cos[r3] + k11 Sin[r1]) +
> X3 (k13 Cos[r1] Cos[r3] + k11 Sin[r1]) -
> C2 (k11 Cos[r1] Sin[r2] +
> k13 (-Cos[r3] Sin[r1] Sin[r2] - Cos[r2] Sin[r3])) +
> X2 (k11 Cos[r1] Sin[r2] +
> k13 (-Cos[r3] Sin[r1] Sin[r2] - Cos[r2] Sin[r3])) -
> C1 (k11 Cos[r1] Cos[r2] +
> k13 (-Cos[r2] Cos[r3] Sin[r1] + Sin[r2] Sin[r3])) +
> X1 (k11 Cos[r1] Cos[r2] +
> k13 (-Cos[r2] Cos[r3] Sin[r1] + Sin[r2] Sin[r3]))))/(-C3 Cos[
> r1] Cos[r3] + X3 Cos[r1] Cos[r3] -
> C2 (-Cos[r3] Sin[r1] Sin[r2] - Cos[r2] Sin[r3]) +
> X2 (-Cos[r3] Sin[r1] Sin[r2] - Cos[r2] Sin[r3]) -
> C1 (-Cos[r2] Cos[r3] Sin[r1] + Sin[r2] Sin[r3]) +
> X1 (-Cos[r2] Cos[r3] Sin[r1] + Sin[r2] Sin[r3]))^2
and I would like to make a C source code of it too... the syntax is not a problem, the optimisation is. C's, k's and r's are just scalars, so I can precompute all of them beforehand. On the other hand, X's are large arrays. I looked all over the internet, played with OptimizeExpression, Hold,..., but nothing helped. I tried your solution then. It looks promising, but there are two problems:
1) k12 and k13 just disappeared... I would believe they just cancelled out as they're kinda additive constants, but 1a) they didn't disappear earlier (and neither from my solution in another system (which is too slow anyway)) 1b) the result is just wrong, and
2) how could I set higher weight/cost (e.g. 1000x) for operations with X's? I guess I need to change the opCost function, but could you give me a hint please?

Thanks a lot in advance! Karel

[biost...@gmail.com]

Update: we got it working, in the end it was just that our problem was so large that the number of aux "vXX" variables was more than 100 so we just needed to extend it to "vXXX"... Now it works and it's quite faster... The other question remains, can I somehow make it worth more operations with X's?

Anyway it is a great piece of code (and an awesome documentation too!), thanks a lot ;)