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.
- does know might come from?
- 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
Post a Comment