Speed up SymPy equation solver

Question

I am trying to solve a set of equations using the following python code (using SymPy of course):

def Solve(kp1, kp2):
    a, b, d, e, f = S('a b d e f'.split())
    equations = [
      Eq(a+b, 2.6),
      Eq(2*a + b + d + 2*f, 7),
      Eq(d + e, 2),
      Eq(a*e,kp2*b*d),
      Eq( ((b * (f**0.5))/a)*((a+b+d+e+f+13.16)**-0.5), kp1)
    ]
    return solve(equations)

The code solves the equations successfully but after approximately 35 seconds. The Solve() function is used inside an iterative block(about 2000 iterations) in another file so speed really matters to me.

Is there any way to speed up the solver? If not can you recommend another way to solve system of equations using python?

Answer 1

Your system of equations can be transformed to a polynomial system in terms of kp1 and kp2 as symbolic parameters. In this case it is advantageous to convert it to a polynomial system so that it can be partially solved in preparation for numerical solving for particular values of the parameters.

The technique her involves computing Groebner bases for which you will want to not have any floats in your equations so let's first make sure that we use exact rational numbers (e.g. use sqrt rather than **0.5 and Rational('2.6') rather than 2.6):

In [20]: kp1, kp2 = symbols('kp1, kp2')

In [21]: a, b, d, e, f = symbols('a, b, d, e, f')

In [22]:     equations = [
    ...:       Eq(a+b, Rational('2.6')),
    ...:       Eq(2*a + b + d + 2*f, 7),
    ...:       Eq(d + e, 2),
    ...:       Eq(a*e,kp2*b*d),
    ...:       Eq( ((b * sqrt(f))/a)/sqrt((a+b+d+e+f+Rational('13.16'))), kp1)
    ...:     ]

In [23]: for eq in equations:
    ...:     pprint(eq)
    ...: 
a + b = 13/5
2⋅a + b + d + 2⋅f = 7
d + e = 2
a⋅e = b⋅d⋅kp₂
              b⋅√f                   
─────────────────────────────── = kp₁
      _________________________      
     ╱                     329       
a⋅  ╱  a + b + d + e + f + ───       
  ╲╱                        25 

The final equation can be trivially transformed to a polynomial: multiply out the denominator and square both sides e.g.:

In [24]: den = denom(equations[-1].lhs)

In [25]: neweq = Eq((equations[-1].lhs*den)**2, (equations[-1].rhs*den)**2)

In [26]: neweq
Out[26]: 
 2      2    2 ⎛                    329⎞
b ⋅f = a ⋅kp₁ ⋅⎜a + b + d + e + f + ───⎟
               ⎝                     25⎠

In [27]: equations[-1] = neweq

In [28]: for eq in equations:
    ...:     pprint(eq)
    ...: 
a + b = 13/5
2⋅a + b + d + 2⋅f = 7
d + e = 2
a⋅e = b⋅d⋅kp₂
 2      2    2 ⎛                    329⎞
b ⋅f = a ⋅kp₁ ⋅⎜a + b + d + e + f + ───⎟
               ⎝                     25⎠

Now that we have a polynomial system we are ready to compute its Groebner basis. We will select a, b, d, e, f as the generators for the Groebner basis and leave kp1 and kp2 as part of the coefficient ring:

In [29]: basis = list(groebner(equations, [a, b, d, e, f]))

This now is a reduced system of equations for abdef with complicated looking coefficients that depend on kp1 and kp2. I'll show one of the equations as an example:

In [30]: basis[-1]
Out[30]: 
                                           ⎛       4    2          2    2          2    ⎞      ⎛       
 4 ⎛   4    2      2    2      2    ⎞    3 ⎜778⋅kp₁ ⋅kp₂    399⋅kp₁ ⋅kp₂    384⋅kp₁    1⎟    2 ⎜102481⋅
f ⋅⎝kp₁ ⋅kp₂  - kp₁ ⋅kp₂  - kp₁  + 1⎠ + f ⋅⎜───────────── - ───────────── - ──────── + ─⎟ + f ⋅⎜───────
                                           ⎝      25              25           25      5⎠      ⎝      6

   4    2            2    2         2               2      ⎞     ⎛             4    2           2    2 
kp₁ ⋅kp₂    15579⋅kp₁ ⋅kp₂    13⋅kp₁ ⋅kp₂   5148⋅kp₁     1 ⎟     ⎜  3799752⋅kp₁ ⋅kp₂    8991⋅kp₁ ⋅kp₂  
───────── + ─────────────── - ─────────── + ───────── + ───⎟ + f⋅⎜- ───────────────── - ────────────── 
25                500              5           125      100⎠     ⎝         3125              625       

          2                2⎞               4    2
  5772⋅kp₁ ⋅kp₂   15984⋅kp₁ ⎟   23853456⋅kp₁ ⋅kp₂ 
- ───────────── - ──────────⎟ + ──────────────────
       125           625    ⎠         15625  

Although these equations might look complicated they have a particular structure which means that once you have substituted values for kp1 and kp2 the final equation depends only on f and all other equations give symbols linearly in terms of f. This means that you can substitute your values and then use either sympy's real_roots or numpys np.roots to get either exact or approximate roots for f, substitute those into your other equations and solve a simple system of linear equations for the remaining unknowns. These are the only steps that would need to be done in the final loop and we can use lambdify to make them happen quickly. The final result then is

from sympy import *

# First solve the system as much as possible in terms of symbolic kp1, kp2
kp1, kp2 = symbols('kp1, kp2')

a, b, d, e, f = symbols('a, b, d, e, f')

# Don't use any floats (yet)
equations = [
    Eq(a+b, Rational('2.6')),
    Eq(2*a + b + d + 2*f, 7),
    Eq(d + e, 2),
    Eq(a*e,kp2*b*d),
    Eq( ((b * sqrt(f))/a)/sqrt((a+b+d+e+f+Rational('13.16'))), kp1)
]

# Convert to full polynomial system:
den = denom(equations[-1].lhs)
neweq = Eq((equations[-1].lhs*den)**2, (equations[-1].rhs*den)**2)
equations2 = equations[:-1] + [neweq]

# Compute the Groebner basis and split the linear equations for a, b, d, e
# from the polynomial equation for f
basis = list(groebner(equations2, [a, b, d, e, f]))
lineqs, eq_f = basis[:-1], basis[-1]

# Solve for a, b, d, e in terms of f, kp1 and kp2
[sol_abde] = linsolve(lineqs, [a, b, d, e])

# Use lambdify to be able to efficiently evaluate the solutions for a, b, d, e
sol_abde_lam = lambdify((f, kp1, kp2), sol_abde)

# Use lambdify to efficiently substitute kp1 and kp2 into the coefficients for eq_f
eq_f_lam = lambdify((kp1, kp2), Poly(eq_f, f).all_coeffs())


def solve_k(kp1val, kp2val):
    eq_f_coeffs = eq_f_lam(kp1val, kp2val)

    # Note that sympy's real_roots function is more accurate and does not need
    # a threshold. It is however slower than np.roots:
    #
    #    p_f = Poly(eq_f_coeffs, f)
    #    f_vals = real_roots(p_f)                  # exact (RootOf)
    #    f_vals = [r.n() for r in real_roots(p_f)] # approximate
    #
    f_vals = np.roots(eq_f_coeffs)
    f_vals = f_vals.real[abs(f_vals.imag) < 1e-10] # arbitrary threshold

    sols = []
    for f_i in f_vals:
        abde_vals = sol_abde_lam(f_i, kp1val, kp2val)
        sol = {f: f_i}
        sol.update(zip([a,b,d,e], abde_vals))
        sols.append(sol)

    return sols

Now you can generate the solutions for particular values of kp1 and kp2 in around 1 millisecond:

In [40]: %time solve_k(2.1, 3.8)
CPU times: user 1.59 ms, sys: 0 ns, total: 1.59 ms
Wall time: 1.51 ms
Out[40]: 
[{a: 49.77432869943, b: -47.174328699430006, d: -0.7687860254920669, e: 2.768786025492067, f: -22.30277
1336968966}, {a: 3.4616794447024692, b: -0.8616794447024684, d: 36.964491584342106, e: -34.964491584342
11, f: -18.01308551452229}, {a: -0.5231222050642267, b: 3.1231222050642273, d: -0.0922228459726158, e: 
2.092222845972616, f: 2.507672525518422}, {a: 0.34146680496125464, b: 2.258533195038746, d: 0.076528664
56901455, e: 1.9234713354309858, f: 1.9910022652348665}]

Caveats: the method above will work in most cases for a polynomial system like this although there are cases where the Groebner basis won't fully separate without the additional use of a splitting variable. Basically just introduce a new variable say t and a random equation like t = a + 2*b + 3*d + 4*e + 5*f. Then make t the last symbol for the Groebner basis and compute t instead of f with np.roots. Also in some cases computing the Groebner basis can be extremely slow but remember that it only ever needs to be computed once (you could save the result in a file if needed).

Another caveat is that squaring to transform into a polynomial system will introduce spurious solutions that don't satisfy the original radical equations. You can check by substituting into the original system. In this case only the last equation matters:

In [41]: sols =  solve_k(2.1, 3.8)

In [42]: equations[-1].subs(sols[0])
Out[42]: -2.10000000000001 = kp₁

In [43]: equations[-1].subs(sols[1])
Out[43]: -2.10000000000006 = kp₁

In [44]: equations[-1].subs(sols[2])
Out[44]: -2.09999999999995 = kp₁

In [45]: equations[-1].subs(sols[3])
Out[45]: 2.09999999999993 = kp₁

So we see that of the 4 solutions returned only the last is correct because it (approximately) satisfies kp1 = 2.1 which was the value passed in. This means a bit of post-processing is needed to get the solutions that you actually want.

Answer 2

Solve the 1st 4 linear equations for a,b,d,e. This generates two solutions. Substitute each of these into the 5th equation (in this form: Eq(b**2*f, d**2*kp1**2*(a + b + 2*d + f + 13.16))) to give two equations (e51 and e52). These two equations are nonlinear in f and, if you use unrad on them you will end up with 2 sixth-order polynomials in f -- which is likely why you get no solution since only quartics are exactly solvable in general. Call those two equations e51u, and e52u.

If you are interested in the real roots, you can use real_roots to give the roots to these polynomials after substituting in the desired values for kp1 and kp2 to leave only f as the unknown. For example,

>>> ans = solve(equations[:-1],a,b,d,e)
>>> l=equations[-1]  # modified as suggested
>>> l=l.lhs-l.rhs
>>> l
b**2*f - d**2*kp1**2*(a + b + 2*d + f + 13.16)
>>> e51=l.subs(ans[0]); e51u=unrad(e51)[0]
>>> e52=l.subs(ans[1]); e52u=unrad(e52)[0]
>>> import random
>>> for r1,r2 in [[random.random() for i in range(2)] for j in range(3)]:
...     print [i.n(2) for i in real_roots(e51u.subs(dict(kp1=r1,kp2=r2)))]
...     print [i.n(2) for i in real_roots(e52u.subs(dict(kp1=r1,kp2=r2)))]
...     print '^_r1,r2=',r1,r2
...     
[1.7, 2.9, 3.0, 8.2]
[1.7, 2.9, 3.0, 8.2]
^_r1,r2= 0.937748743197 0.134640776315
[1.3, 2.3, 4.7, 7.4]
[1.3, 2.3, 4.7, 7.4]
^_r1,r2= 0.490002815309 0.324553144174
[1.1, 2.1]
[1.1, 2.1]
^_r1,r2= 0.308803300429 0.595356213169

It appears that e51u and e52u both give the same solutions, so perhaps you only need to use one of them. And you should check the answers in the original equations to see which one(s) are true solutions:

>>> r1,r2
(0.30880330042869408, 0.59535621316941589)
>>> [e51.subs(dict(kp1=r1,kp2=r2,f=i)).n(2) for i in real_roots(e51u.subs(dict(kp1=r1,kp2=r2)))]
[1.0e-12, 13.]

So here you see that only the first solution (seen to be 1.1 from above) is actually a solution; 2.1 was a spurious solution.

Answer 3

You only need to solve the equation one time. After that you will have equations of the form:

• a = f_1(kp1, kp2)
• b = f_2(kp1, kp2)
• ...

so you can simply compute a, ..., e in dependency of kp1 and kp2. For example solving the first, third and fourth equations gives:

• b: -a + 2.6
• e: 2.0*kp2*(5.0*a - 13.0)/(5.0*a*kp2 - 5.0*a - 13.0*kp2),
• d: 10.0*a/(-5.0*a*kp2 + 5.0*a + 13.0*kp2)

Solving all five equations on my pc is too slow, but if it gives you an expression you just have to plug in (substitute) the values for kp1 and kp2, you don't have to solve the equations again. For substitution take a look at the sympy documentation.

So your loop should look like:

solutions = sympy.solve(eqs, exclude=[kp1, kp2])
for data_kp1, data_kp2 in data:
    for key, eq in solutions:
        solution = eq.subs([(kp1, data_kp1), (kp2, data_kp2)])
        final_solutions.append("{key}={solution}".format(key=key, solution=solution))