On Sun, 24 Aug 2025 at 12:01, Gábor Horváth <[email protected]> wrote:
>
>
> > Could you post an example of a large system?
>
> Attached one with 239 equations to this email. Not sure if attached files
> go through, though. If not, I can paste it directly into the email body,
> it's size is 17K.
Yes, it has come through.
Okay looking at this you definitely want a special algorithm for these
kinds of systems rather than just using a general purpose solver.
Ideally SymPy would have a reduce function and if it did we could
consider whether it should be able to handle larger systems of this
kind efficiently.
First find the Max expressions and any symbols that are in them:
maxes = Tuple(*eqs).atoms(Max)
special = Tuple(*maxes).free_symbols
Then replace all Max expressions with symbols themselves:
syms_max = symbols(f'm:{len(maxes)}')
rep_m = dict(zip(maxes, syms_max))
eqs2 = [eq.xreplace(rep_m) for eq in eqs]
Now you have a linear system of equations for the original variables
v0, v1 etc and also new variables m0, m1 representing the maxes. We
want to get the smallest possible system involving only the maxes, the
symbols that are in the maxes (special) and as few of the other
symbols as possible.
Choose the ordering of the symbols carefully and ask linsolve to
reduce the underdetermined system:
dontcare = [s for s in syms if s not in special]
unknowns = list(syms_max) + dontcare + list(special)
[sol] = linsolve(eqs2, unknowns)
In so far as possible this tries to solve for syms_max in terms of
everything else and to solve for everything in terms of the symbols in
special. It succeeds in expressing all of the maxes in terms of
everything else:
>>> sol.atoms(Symbol) - set(syms)
set()
Mostly this expresses everything in terms of the special symbols
although a few others still appear as free parameters:
>>> print(Tuple(*sol).free_symbols)
{v232, v238, v121, v223, v183, v197, v62, v3, v218, v189, v106, v214,
v216, v231, v11, v213, v6, v88, v71, v120, v111, v180, v21, v129,
v208, v44, v153, v192, v201, v184}
>>> print(Tuple(*sol).free_symbols - special)
{v238, v197, v129, v223, v192, v214}
Everything is now expressed in terms of these 30 symbols of which 24
are symbols appearing in the maxes and 6 are not. These 6 symbols
appearing in maxes have been solved purely from the linear equations
but are given in terms of the 30 symbols left over:
>>> print(special - Tuple(*sol).free_symbols)
{v25, v9, v154, v160, v229, v84}
I think at this point then you can translate the solution into 30
equations for the remaining 30 unknowns:
rep = dict(zip(unknowns[30:], sol[30:]))
eqs3 = [Eq(m.xreplace(rep), s) for s, m in zip(sol[:30], maxes)]
Now you have 30 equations for 30 unknowns involving 30 maxes and the
other equations and variables are all eliminated and solved
The other equations after these 30 can be inspected to give bounds on
the 30 variables left over. You have a parametric solution for the 209
other variables in terms of the 30 symbols. You can substitute that
solution into any auxiliary inequality constraints to further
constraint the values of the 30 symbols (or prove unsatisfiability).
These are the remaining Maxes:
In [222]: for m in Tuple(*eqs3).atoms(Max): print(m.n(3))
Max(0.802, v111 - v153 + v180)
Max(0.757, v62)
Max(0.597, 49.8*v153 - 30.8*v180 - 54.3*v183 + 210.0*v189 - 210.0*v201
+ 15.4*v21 - 130.0*v218 - 4.23*v3 - 21.7*v44 + 210.0*v88 - 33.9)
Max(0.879, v189 - v201 + v88)
Max(0.732, v21)
Max(0.735, v216)
Max(0.767, v44)
Max(0.645, v231)
Max(0.721, v184)
Max(0.789, v121)
Max(0.897, v189)
Max(0.741, v3)
Max(0.835, v153)
Max(0.827, v88)
Max(0.677, -100.0*v153 + 164.0*v180 + 46.8*v183 - 183.0*v189 +
183.0*v201 - 46.9*v21 + 128.0*v218 + 2.33*v3 - 11.2*v44 - 183.0*v88 +
1.07)
Max(0.846, v201)
Max(0.786, v183)
Max(0.783, v6)
Max(0.727, v213)
Max(0.719, v71)
Max(0.78, v111 - v153 + v218)
Max(0.658, v208)
Max(0.804, 2.92*v153 - 6.78*v180 - 0.734*v183 + 6.25*v189 - 1.72*v201
+ 1.16*v21 - 1.2*v218 - 0.0219*v3 - 0.123*v44 + 1.72*v88 - 0.473)
Max(0.703, v11)
Max(0.786, v120)
Max(0.683, v232)
Max(0.839, v218)
Max(0.837, v106)
Max(0.862, v180)
Max(0.777, v111)
These are given linearly in terms of the same symbols e.g.:
>>> print(eqs3[0].n(3))
Eq(Max(0.757, v62), 474.0*v106 + 0.503*v111 - 81.4*v120 - 312.0*v121
- 953.0*v153 + 1.8e+3*v180 + 506.0*v183 - 2.1e+3*v189 + 2.03e+3*v201 -
1.91*v208 - 543.0*v21 - 0.416*v216 + 1.53e+3*v218 + 7.85*v232 -
1.58*v238 - 71.3*v3 - 213.0*v44 - 13.0*v6 + 28.1*v62 - 2.03e+3*v88 -
61.2)
Maybe we are actually putting this the wrong way round. Keep the m
symbols and invert to solve for the remaining symbols in terms of the
Maxes:
eqs4 = [Eq(m, s) for s, m in zip(sol[:30], syms_max)]
[sol2] = linsolve(eqs4, syms2)
Yes, that's better. Now you have solved for 209 symbols in terms of
these 30 and you have solutions for these 30 in terms of the maxes.
Now each time you pick a branch for a Max it directly tells you the
value of some symbol which you can then eliminate from everything
else. There is no further need to call linsolve because the linear
equations are all solved and it is just a case of considering the
maxes in turn.
Or something like that...
I guess this is already where you are at and you just need a more
efficient way of pruning the 2^30 branches.
--
Oscar
--
You received this message because you are subscribed to the Google Groups
"sympy" group.
To unsubscribe from this group and stop receiving emails from it, send an email
to [email protected].
To view this discussion visit
https://groups.google.com/d/msgid/sympy/CAHVvXxSQFzjLo9WR6fMJna%3D4Zbq%3DVLKYtEkKKT2eevffzcN-_g%40mail.gmail.com.
