Slow symbolic matrix substitution with sympy

Question

Im working with sympy on a jxobian J character matrix of size QxQ . Each coefficient of this matrix contains Q symbols, from f[0] to f[Q-1] . What Id like to do is replace each character in each coefficient J known values g[0] to g[Q-1] (which are no longer characters). The fastest way I've found is this:

 for k in range(Q): J = J.subs(f[k], g[k]) 

However, I find this "main" operation for a very long time! For example, using this MCVE:

 import sympy import numpy as np import time Q = 17 f0, f1, f2, f3, f4, f5, f6, f7, f8, f9, f10, f11, f12, f13, f14, f15, f16 = \ sympy.symbols("f0 f1 f2 f3 f4 f5 f6 f7 f8 f9 f10 f11 f12 f13 f14 f15 f16") f = [f0, f1, f2, f3, f4, f5, f6, f7, f8, f9, f10, f11, f12, f13, f14, f15, f16] e = np.array([0., 1., 0., -1., 0., 1., -1., -1., 1., 2., -2., -2., 2., 3., 0., -3., 0.]) u = np.sum(f * e) / np.sum(f) function = e * np.sum(f) * (1. + u * e + (u * e)**2 - u * u) F = sympy.Matrix(function) g = e * (1. + 0.2 * e + (0.2 * e)**2) start_time = time.time() J = F.jacobian(f) print("--- %s seconds ---" % (time.time() - start_time)) start_time = time.time() for k in range(Q): J = J.subs(f[k], g[k]) print("--- %s seconds ---" % (time.time() - start_time)) 

substitution takes about 5 seconds on my computer, and calculating the jacobian matrix takes only 0.6 s. In another (longer) code, the replacement takes 360 with Q=37 (while 20 seconds for the calculation in the Jacobian language)!

Also, when I look at my running processes, I see that the Python process sometimes stops working during matrix replacement.

• Does anyone know where this might come from?
• Is there any way to make this operation faster?

Answer 1

You might want to try Theano . It implements a jacobian function that is parallel and faster than sympy .

The next version achieves an acceleration of 3.88! Now the substitution time is inferior to the second.

 import numpy as np import sympy as sp import theano as th import time def main_sympy(): start_time = time.time() Q = 17 f = sp.symbols(('f{} ' * Q).format(*range(Q))) e = np.array([0., 1., 0., -1., 0., 1., -1., -1., 1., 2., -2., -2., 2., 3., 0., -3., 0.]) u = np.sum(f * e) / np.sum(f) ue = u * e phi = e * np.sum(f) * (1. + ue + ue*ue - u*u) F = sp.Matrix(phi) J = F.jacobian(f) g = e * (1. + 0.2*e + (0.2*e)**2) for k in range(Q): J = J.subs(f[k], g[k]) print("--- %s seconds ---" % (time.time() - start_time)) return J def main_theano(): start_time = time.time() Q = 17 f = th.tensor.dvector('f') e = np.array([0., 1., 0., -1., 0., 1., -1., -1., 1., 2., -2., -2., 2., 3., 0., -3., 0.]) u = (f * e).sum() / f.sum() ue = u * e phi = e * f.sum() * (1. + ue + ue*ue - u*u) jacobi = th.gradient.jacobian(phi, f) J = th.function([f], jacobi) g = e * (1. + 0.2*e + (0.2*e)**2) Jg = J(g) # evaluate expression print("--- %s seconds ---" % (time.time() - start_time)) return Jg J0 = np.array(main_sympy().tolist(), dtype='float64') J1 = main_theano() print(np.allclose(J0, J1)) # compare results