python - Slow substitution of symbolic matrix with sympy -


i’m working sympy on symbolic jacobian matrix j of size qxq. each coefficient of matrix contains q symbols, f[0] f[q-1]. i’d substitute every symbol in every coefficient of j known values g[0] g[q-1] (which no more symbols). fastest way found follows:

for k in range(q):     j = j.subs(f[k], g[k]) 

however, find "basic" operation long! example, mcve:

import sympy import numpy 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() k in range(q):     j = j.subs(f[k], g[k]) print("--- %s seconds ---" % (time.time() - start_time)) 

the substitution takes 5s on computer, while computation of jacobian matrix takes 0.6s. on (longer) code, substitution takes 360s q=37 (while 20s jacobian computation)!

moreover, when @ running processes, can see python process stops working during matrix substitution.

  1. does know might come from?
  2. is there way make operation faster?

you might want try theano that. implements jacobian function parallel , faster sympy.

the following version achieves speedup of 3.88! substitution time inferior second.

import numpy np import sympy sp import theano 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)      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 

Comments

Popular posts from this blog

Is there a better way to structure post methods in Class Based Views -

performance - Why is XCHG reg, reg a 3 micro-op instruction on modern Intel architectures? -

jquery - Responsive Navbar with Sub Navbar -