Substitution Optimization Question

[wflynny]

Hi!

For my current project, I start with a pretty complicated symbolic
expression that I want to numerically evaluate for different values of
different variables. Say I have an expression

expr(x,y,z) = A*x+B*y*z+DiracDelta(x-y+z)+...

and an array of numerical values for each variable, x,y,z:

x = np.array([1.0,2.0,3.0,4.0,5.0,6.0,7.0,8.0,9.0,10.0])
y = np.array([0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1.0])
z = np.array([1.0,1.5,2.0,2.5,3.0,3.5,4.0,4.5,5.0,5.5]).

Currently I want to evaluate expr(x,y,z) for each combination of
x,y,z, which yields an array of 1000 numerical values of expr, in this
example. This means a set of for loops scanning over all 1000
combinations. My question:

Is there is an easier (less-cpu-cycles-) way to do this kind of mass
substitution in sympy?

Can we map in an efficient way to do this? Our problem is that we want
to evaluate our expression at tens of thousands of points and the for
loops just eat up too much time.

Any suggestions would be greatly appreciated!

Bill

[Aaron S. Meurer]

In the end, you have to do 1000 computations one way or another. Generally, implicit maps are faster than for loops in Python, so maybe three nested map() functions would be faster.

Also, it would depend on what the expression is, but in some cases you might be able to get faster results by changing the order of substitution. For example, if you wanted to map 0, 1, and 2 onto x, y and z in x*y*z, you should first substitute x for each number, then y, then z, because as soon as you do the first substitute for 0, the expression will vanish and the substitutions for y and z will be trivial. If anything can make the expression evaluate simpler like this, you would want it in the outer loop so sympy only has to evaluate the simpler version once.

Maybe someone else here will have a more exact solution for you though.

Aaron Meurer

> --
>
> You received this message because you are subscribed to the Google Groups "sympy" group.
> To post to this group, send email to sy...@googlegroups.com.
> To unsubscribe from this group, send email to sympy+un...@googlegroups.com.
> For more options, visit this group at http://groups.google.com/group/sympy?hl=en.
>
>

[Ronan Lamy]

There's no optimised way of doing this inside sympy but it's an obvious
use case for numpy. So you should use sympy.lambdify and then build your
array in numpy via brodcasting. Here's an example session to explain
what I mean:

In [1]: from numpy import arange, newaxis

In [2]: expr = x**2 + 3*sin(y*z) + cosh(x-y+z)

In [3]: x_vals = arange(10.)

In [4]: y_vals = arange(0, 100, 10)

In [6]: z_vals = 1 - (x_vals/10)**3

# This returns a function that uses numpy for its computations, and
# returns expr with the values of x, y, z replaced by its arguments.
In [7]: expr_np = lambdify((x, y, z), expr, modules = "numpy")

# I think that a nested list comprehension is the best (i.e. clearest
# and fastest) way to build the array in "pure" sympy.
In [9]: res_sym = [[[expr.subs(((x,a), (y,b), (z,c))) for c in z_vals]
for b in y_vals] for a in x_vals]

# x_vals, y_vals, z_vals are each mapped as a row/column along a
# different axis. Then, magic (AKA broadcasting) happens!
# The result is a 3D array equivalent to res_sym.
In [11]: res_np = expr_np(x_vals[:, newaxis, newaxis],
y_vals[newaxis, :, newaxis], z_vals[newaxis, newaxis, :])

In [12]: res_sym[0][0]
Out[12]: [cosh(1.0), cosh(0.999), cosh(0.992), cosh(0.973), cosh(0.936),
cosh(0.875), cosh(0.784), cosh(0.657), cosh(0.488), cosh(0.271)]

In [13]: map(N, res_sym[0][0])
Out[13]:
[1.54308063481524, 1.54190620496611, 1.5337283038256, 1.51190883426005,
1.47097771056081, 1.4078686568228, 1.32339583962698, 1.223700435959
02, 1.1214538615077, 1.03694578339483]

In [14]: res_np[0,0]
Out[14]:
[ 1.54308063 1.5419062 1.5337283 1.51190883 1.47097771 1.40786866
1.32339584 1.22370044 1.12145386 1.03694578]

In [15]: %timeit [[[expr.subs(((x,a), (y,b), (z,c))) for c in z_vals]
for b in y_vals] for a in x_vals]
10 loops, best of 3: 35.4 ms per loop

In [16]: %timeit expr_np(x_vals[:, newaxis, newaxis], y_vals[newaxis, :,
newaxis], z_vals[newaxis, newaxis, :])
1000 loops, best of 3: 313 us per loop


Bottom line: it's 100 times faster with numpy.

However, not all sympy functions can be mapped to numpy, which might
cause some difficulties for you.

Ronan

[wflynny]

Hey!

Thank you so much for your comments! They were really helpful but also
leave me with some more questions. For example, say you're given an
expression

>>> expr = A*x**2 + B*x + C

and you are also given arrays of values

>>> A = np.array([...])
>>> B = np.array([...])
>>> C = np.array([...])

but no values for x. Thus x is unknown and you should just get an
array of polynomials in x. Is there any way to handle this with
lambdify? I run into problems if I only try to halfway evaluate an
expression.

Here's an example

In [9]: from numpy import newaxis as NAX

In [10]: A,B,C,x = sp.symbols('ABCx')

In [11]: AA,BB,CC = np.linspace(1,10,5),np.linspace(1,10,5),np.linspace
(1,10,5)

In [12]: expr = A*x**2 + B*x + C

In [13]: eval_expr = sp.lambdify((A,B,C),expr,modules="numpy")

In [14]: np_eval = eval_expr(AA[:,NAX,NAX],BB[NAX,:,NAX],CC
[NAX,NAX,:])
---------------------------------------------------------------------------
NameError Traceback (most recent call
last)

C:\Documents and Settings\wflynn\My Documents\Python\<ipython console>
in <module>()

C:\Documents and Settings\wflynn\My Documents\Python\<string> in
<lambda>(A, B,
C)

NameError: global name 'x' is not defined

[wflynny]

Nevermind! I just figured it out:

In [9]: from numpy import newaxis as NAX

In [10]: A,B,C,x = sp.symbols('ABCx')

In [11]: AA,BB,CC = np.linspace(1,10,5),np.linspace(1,10,5),np.linspace
(1,10,5)

In [12]: expr = A*x**2 + B*x + C

...

In [17]: eval_expr2 = sp.lambdify((A,B,C,x),expr,modules="numpy")

In [18]: np_eval = eval_expr2(AA[:,NAX,NAX],BB[NAX,:,NAX],CC
[NAX,NAX,:],sp.Symbol('x'))

In [22]: np_eval[0][0]
Out[22]:
array([1.0 + x + x**2, 3.25 + x + x**2, 5.5 + x + x**2, 7.75 + x +
x**2,
10.0 + x + x**2], dtype=object)

Thank you so much for your help guys!

Bill