Optimizing Expressions

Using the Optimize.m package

Terry Robb
<tdr@vaxc.cc.monash.edu.au>

There is sometimes a need to improve, or otherwise optimize, the way in which some expressions are evaluated.

Consider for example a function definition like h[x_] := Erf[x] + Erf[x]^2 and the way the expression on the right is evaluated for a given numerical x. The Erf function will be evaluated twice (Share[] not-with-standing), which is quite unnecessary in this case. If however we had defined the function as h[x_] := Block[{e=Erf[x]}, e+e^2], then Erf would be evaluated just once and that is a more optimized evaluation.

The package Optimize.m implements a command called Optimize to automate the above type of expression optimization.

Needs["Optimize`"];

What Optimize does is replace an expression with a HoldBlock object, which is a device for creating a block of code that will be held unevaluated. The idea is that HoldBlock will ultimately be used inside a function or a Compile definition etc, where the hold can be released and the optimized code inserted.

Here are two simple examples which illustrate a typical use of Optimize:

f[x_Real] := Optimize[Sin[x] + Cos[Sin[x]]];
?f

and

g = Compile[x, Optimize[Sin[x] + Cos[Sin[x]]]]

In each case Optimize evaluates itself (an Evaluate is not required here, because  Optimize insists on that anyway in this context) and optimizes the expression Sin[x] + Cos[Sin[x]] so as to avoid the unnecessary double evaluation of Sin[x]. Thus when you call f[2.2] or g[2.2], the Sin[2.2] will be evaluated just once and not twice.

Note that only basic structural optimization is performed – in other words only if two or more sub-expressions have literally the same structural form will any optimization take place. No mathematical knowledge is therefore used in this procedure. See example 6 below for a simple way to effectively add some extra knowledge though.

Commands

There are three routines defined in Optimize.m, and here's what they do:

?Optimize

?HoldBlock

?Cost

Example 1

Here's a fairly long expression:

expr = Nest[(#+Sin[#])&, x, 10];

Here's how we would define an optimized function to evaluate it:

f1[x_Real]  :=  Optimize[expr];

And here's how we would compile the optimized expression:

f2 = Compile[x, Optimize[expr]];

Now let's use the un-optimized expression in a command to see how fast it is:

Plot[expr//Evaluate, {x, 0, 0.01},
     DisplayFunction->Identity]; // Timing

Now let's use the optimized expression:

Plot[f1[x], {x, 0, 0.01},
     DisplayFunction->Identity]; // Timing

And finally let's use the compiled optimized expression:

Plot[f2[x], {x, 0, 0.01},
     DisplayFunction->Identity]; // Timing

In each case the timings get better.

To see what optimization has taken place use:

?f1

Note the cost of the original expression is:

Cost[expr]

while the cost of the optimized expression is just:

Cost[Optimize[expr]]

The optimization has reduced the 1023 Sin evaluations down to just 10, by finding 17 rules.

Example 2

Here's a big expression which we optimize:

expr = Expand[(Sin[x*y] + Log[x] + Cosh[x^-y]^5)^10];
oexpr = Optimize[expr];

The time to evaluate the expression at a typical point is:

(expr        /. {x->1.1, y->0.5}) // Timing

But the optimized expression evaluates a little faster:

Normal[oexpr /. {x->1.1, y->0.5}] // Timing

The expression costs indicate why the latter is faster:

Cost[expr]

Cost[oexpr]

Example 3

You can use Optimize to show structural information that is not otherwise obvious in the output from some commands such as Solve, D, and Integrate etc. For example:

Optimize[z /.
       Solve[z^4 - 12*z^3 - 16*z^2 + 4*z + c == 0, z]
        ]//InputForm

Example 4

The symbols e, f and g which are defined below, each optimize the expression Sin[x] + Cos[Sin[x]] in slightly different contexts:

e = Optimize[Sin[x] + Cos[Sin[x]]]

f[x_Real] := Optimize[Sin[x] + Cos[Sin[x]]]

g = Compile[x, Optimize[Sin[x] + Cos[Sin[x]]]]

The three methods can be evaluated at the point say x->2.2 as follows:

{Normal[e /. x->2.2], f[2.2], g[2.2]}

The evaluation cost of the un-optimized code is:

Cost[Sin[x] + Cos[Sin[x]]]

While the evaluation cost of the optimized code is:

Cost[Optimize[Sin[x] + Cos[Sin[x]]]]

Example 5

Big expressions like the following are a common application of Optimize

expr = Sum[(Sin[x] + Cos[Sin[x]])^n, {n, 1, 100}];

The initial cost is large

Cost[expr]

But the optimized cost is much better

Cost[Optimize[expr]]

Example 6

Here's a simple expression:

expr = (1+x*y)/(1-x*y)

Let's try optimize it:

Optimize[expr]

No optimization however takes place here because Optimize only sees structural stuff, and so doesn't understand anything about the mathematics of multiplication. By looking at the FullForm of the expression we see why Optimize doesn't optimize the expression:

FullForm[expr]

No structures are explicitly duplicated in the above, and so no optimization is performed. We can however tell Optimize one fact about multiplication that will let it see more structure. By turning Times into a right associative binary operator (rather than a flat orderless operator), extra structural information appears.

To do this we repeatedly apply the rule Times[a_,b__] :> TIMES[a,Times[b]] to turn Times into a right associative binary operator named TIMES. Then we call Optimize on that new structure, and at the end replace TIMES with Times to get back the normal Times operator. Thus

Optimize[expr //. Times[a_, b__] :> TIMES[a, Times[b]]
        ] /. TIMES -> Times

The above has now extracted out the common x*y term.

To make the above easier a second argument named Times->True is provided with Optimize, so that you can get the same effect by simply typing

Optimize[expr, Times->True]

In order to save a small amount of initilization time, this feature has not been implemented as a normal option however, and so Options[Optimize] is not modifiable.

Caveats

There is one slight caveat with using Optimize. It doesn't understand about expressions that contain anything with a HoldAll (or HoldFirst or HoldRest) attribute. Tweak the source code yourself if you desparately need to optimize such objects.

Also this implementation of Optimize works well if a smallish number (say less than 500) of the optimization variables (the O$ variables) are ultimately needed. The worst case is if the number of variables is near half the number of basic components in a very large expression, but such cases are rare. If efficiency were a major issue then Optimize would be best implemented in C. The Mathematica version does however work well for the common cases.