solve is too slow for small systems of equations

[mtomassoli]

The following problem is solved almost instantaneously in both SageMath and Matlab, but it takes much longer in sympy. The problem has 8 solutions. If I use manual=True, solve is very fast, but it finds only 2 solutions.
If I try a bigger problem with the same structure, it takes forever in sympy but only a few seconds in SageMath and Matlab.

from sympy import *

init_printing()

t, tp, tq, ptq2, D, eps = symbols('t tp tq ptq2 D epsilon', nonnegative=True)
a, b, s = symbols('a b s', nonnegative=True)
mu_a, mu_b, mu_s = symbols('mu_a mu_b, mu_s', nonnegative=True)

L = (D*s*s + a*a + b*b - 2*s*t - 2*a*tp - 2*b*tq + 2*s*a + 2*s*b + 2*a*b*ptq2)
L += eps*(a-b)*(a-b)
L += -mu_a*a - mu_b*b - mu_s*s

primal_vars = [a, b, s]
dual_vars = [mu_a, mu_b, mu_s]

eqs = [diff(L, var) for var in primal_vars]
eqs += [dv * pv for pv, dv in zip(primal_vars, dual_vars)]

solve(eqs, primal_vars + dual_vars)

[oscarbenjamin]

Here are the equations and unknowns:

In [63]: for e in eqs: pprint(e)
2⋅a + 2⋅b⋅ptq₂ + ε⋅(2⋅a - 2⋅b) - μₐ + 2⋅s - 2⋅tp
2⋅a⋅ptq₂ + 2⋅b + ε⋅(-2⋅a + 2⋅b) - μ_b + 2⋅s - 2⋅tq
2⋅D⋅s + 2⋅a + 2⋅b - μₛ - 2⋅t
a⋅μₐ
b⋅μ_b
μₛ⋅s

In [64]: pprint(primal_vars + dual_vars)
[a, b, s, μₐ, μ_b, μₛ]

If I read this correctly the first 3 equations are linear. Then the remaining equations are fairly trivially solved by saying that one or other variable must be zero. Presumably the 8 solutions come from making that choice in each of the last 3 equations.

Although in principle this is a multivariate polynomial system it's easy to see heuristics that would be efficient for this kind of problem. Perhaps recognising equations of the form x*y = 0 where x and y are both unknowns would be a useful strategy...

[mtomassoli]

Yes, you're correct: that's equivalent to solving 8 linear problems.
I add non-negativity constraints to avoid quadratic cubic terms in the lagrangian.

Here's a problem with quadratic cubic terms where sympy only finds a trivial solution:

from sympy import *

t, tp, tq, tpq, ptq, D = symbols('t tp tq tpq ptq D')
a, b, s = symbols('a b s')

dtd = a*a + b*b + 2*a*b*ptq
L = D*s*s + dtd*dtd - 2*s*t - 2*a*a*tp - 2*b*b*tq - 4*a*b*tpq + 2*s*dtd

primal_vars = [a, b, s]

eqs = [diff(L, var) for var in primal_vars]
solve(eqs, primal_vars)

I'll have to solve it "manually" by exploiting its structure...

edit: Sorry for the many edits. This was the last one.

[smichr]

the 8 solutions come from making that choice in each of the last 3 equations

explicitly:

>>> list(cartes((a,ma),(b,mb),(s,ms)))  # list of sym-dependent factors of equations that are products
[(a, b, s), (a, b, ms), (a, mb, s), (a, mb, ms), (ma, b, s), (ma, b, ms), (ma, mb, s), (ma, mb, ms)]
>>> len(_)
8