Potential sympy.solve error?

Gábor Horváth

unread,

Aug 9, 2025, 8:26:28 AMAug 9

to sympy

Hi Team!

I'm using sympy to programmatically solve systems of equations, that are essentially linear, but some equations contain Max or Min in some easy manner. However, sympy.solve seems to throw a NotImplementedError in such cases. This makes sense, if currently Min and Max are not supported by sympy.solve.

That said, if I rewrite the equations to Piecwise, then it can solve the system in many cases. But sometimes it throws an Invalid NaN comparison somewhere in relational.py. See below for a minimal example that reproduces the issue (system is x=Max(1, y), y=Max(2, x)).

I'm not quite sure if this is the intended behaviour or not, but I don't know enough of the inner workings of sympy to identify the direct cause of the NaN comparison exception. The solution I'm looking for here is x=y>=2. Now if I only give one equation (say, x=Max(1,x)), then the Piecewise version throws me a NotImplementedError, but also tells me that the solution is likely containing the region x>1, which is helpful.

Finally, solveset corretly identifies the interval solution for x=Max(1,x) for reals, but I don't know is solveset can be directly called somehow on a system or not.

So my questions are:

1) Is the expected behaviour for sympy.solve for the Piecewise case to throw a NaN comparison exception, or is this a sign of some underlying issue in sympy.solve for Piecewise functions?

2) Is sympy.solve expected to rewrite Min and Max to Piecewise at some point in the future? I found some piece of code in solvers/solvers.py:944-957 which does some sort of rewriting for Abs, but didn't find any for Min or Max.

3) From some reading on the sympy docs page, I figured solveset is expected to be the future solver. Is there a direct way I can apply solveset on a system, rather than on one equation?

4) If the answer to 3) is yes, is there also a way which provides a stable output that I can use in a programmatic way?

Thanks for the help!

Gábor

#######################################

## Minimal working example reproducing the error: ##

#######################################

import sympy
x, y = sympy.symbols('x y', real=True)
sys = [sympy.Eq(x, sympy.Max(y, 1)), sympy.Eq(y, sympy.Max(2, x))]
variables = (x, y)
sympy.solve(sys, variables)
new_sys = [eq.rewrite(sympy.Piecewise) for eq in sys]
sympy.solve(new_sys, variables)

Oscar Benjamin

unread,

Aug 9, 2025, 8:43:57 AMAug 9

to sy...@googlegroups.com

On Sat, 9 Aug 2025 at 13:26, Gábor Horváth <hungabo...@gmail.com> wrote:
>
> So my questions are:
> 1) Is the expected behaviour for sympy.solve for the Piecewise case to throw a NaN comparison exception, or is this a sign of some underlying issue in sympy.solve for Piecewise functions?

No, that is just a bug.

> 2) Is sympy.solve expected to rewrite Min and Max to Piecewise at some point in the future? I found some piece of code in solvers/solvers.py:944-957 which does some sort of rewriting for Abs, but didn't find any for Min or Max.

It probably does make sense to rewrite Min and Max as Piecewise if
they contain the variables of interest.

> 3) From some reading on the sympy docs page, I figured solveset is expected to be the future solver. Is there a direct way I can apply solveset on a system, rather than on one equation?

This documentation needs to be changed. The solveset function is not
going to replace solve. It was originally intended as a replacement
for solve but was then misdesigned as something different that cannot
replace solve because it does not provide comparable functionality or
interface.

> 4) If the answer to 3) is yes, is there also a way which provides a stable output that I can use in a programmatic way?

Only solve with dict=True gives stable output that is suitable for
programmatic use but it is not really designed to handle systems of
the type that you have.

The solutions to the system x = Max(y, 1), y = Max(2, x) are y = x for
x >= 2. Given that the solution set is constrained by an inequality
I'm not sure what form of output you would want for programmatic use.

It is very difficult to make a solve function that can take arbitrary
inputs and give something useful as output, especially handling
underdetermined cases, inequality constraints and so on. However if
you have a specific class of problems like these with min and max it
is a lot easier to make something that handles a more limited case. In
this example you can split the equations into 4 linear systems and
solve those.

In general I think that what you should really want in a situation
like this is something like Mathematica's Reduce function rather than
Solve. Unfortunately SymPy has several versions of solve but nothing
like reduce.

--
Oscar

Gábor Horváth

unread,

Aug 9, 2025, 12:31:30 PMAug 9

to sy...@googlegroups.com

Hi Oscar,

Thanks for replying!

> On Sat, 9 Aug 2025 at 13:26, Gábor Horváth <hungabo...@gmail.com> wrote:
>>
>> So my questions are:
>> 1) Is the expected behaviour for sympy.solve for the Piecewise case to throw a NaN comparison exception, or is this a sign of some underlying issue in sympy.solve for Piecewise functions?
>
> No, that is just a bug.

Ok, thanks. Should I report it in some official way somewhere (if yes,
where and how), or has this email already taken care of the bugreport?
Maybe I should open an issue on https://github.com/sympy/sympy/issues ?
(As I don't fully understand the code behind the Abs rewrite, I don't feel
confident enough to work on this on my own and submit a PR.)

>> 2) Is sympy.solve expected to rewrite Min and Max to Piecewise at some point in the future? I found some piece of code in solvers/solvers.py:944-957 which does some sort of rewriting for Abs, but didn't find any for Min or Max.
>
> It probably does make sense to rewrite Min and Max as Piecewise if
> they contain the variables of interest.

Ok, thanks. Same question as above: should I open some issue on
https://github.com/sympy/sympy/issues or elsewhere, or this is already
taken care of by this email?

>> 3) From some reading on the sympy docs page, I figured solveset is expected to be the future solver. Is there a direct way I can apply solveset on a system, rather than on one equation?
>
> This documentation needs to be changed. The solveset function is not
> going to replace solve. It was originally intended as a replacement
> for solve but was then misdesigned as something different that cannot
> replace solve because it does not provide comparable functionality or
> interface.

Ah, I see. This was the page I found first:
https://lidavidm.github.io/sympy/modules/solvers/solveset.html
and only later did I find the official documentation:
https://docs.sympy.org/latest/modules/solvers/solveset.html
Now the two are pretty much (looks to be more than 90%) identical, but the
first has this section that confused me about solveset potentially
replacing solve:

"What will you do with the old solve?

There are still a few things solveset can’t do, which the old solve
can, such as solving non linear multivariate & LambertW type equations.
Hence, it’s not yet a perfect replacement for old solve. The ultimate goal
is to:

Replace solve with solveset once solveset is at least as powerful
as solve, i.e., solveset does everything that solve can do currently, and
eventually rename solveset to solve."

This section is missing from the official documentation page on
docs.sympy.org, but considering google search shows the second page, as
well, it might make sense to remove or rephrase this section there.

BTW, are there some old threads somewhere maybe where such discussions
can be read? I would be curious what you mean by misdesign, and what were
the discussion and decision points around solve an solveset.

There is also mention in the docs about FiniteSet handling tuples
corresponds to a finite subset of points in the multi-dimensional space.
That suggests that solveset was intended to be able to solve multi-variate
equations or even systems. Is that indeed the case, or only some special
cases (like linsolve) can really handle multi-variate equations/systems?

>> 4) If the answer to 3) is yes, is there also a way which provides a stable output that I can use in a programmatic way?
>
> Only solve with dict=True gives stable output that is suitable for
> programmatic use but it is not really designed to handle systems of
> the type that you have.
>
> The solutions to the system x = Max(y, 1), y = Max(2, x) are y = x for
> x >= 2. Given that the solution set is constrained by an inequality
> I'm not sure what form of output you would want for programmatic use.
>
> It is very difficult to make a solve function that can take arbitrary
> inputs and give something useful as output, especially handling
> underdetermined cases, inequality constraints and so on. However if
> you have a specific class of problems like these with min and max it
> is a lot easier to make something that handles a more limited case. In
> this example you can split the equations into 4 linear systems and
> solve those.
>
> In general I think that what you should really want in a situation
> like this is something like Mathematica's Reduce function rather than
> Solve. Unfortunately SymPy has several versions of solve but nothing
> like reduce.

Hm, all are very good points. I have some very crude code to solve these
systems exactly as you suggest, by splitting at the various Min and Max
occurrences, solving the underlying linear system (with linsolve), and
then check if the obtained solutions indeed solve the original system with
Min and Max.

That said, somehow my code is rather slow, especially when there are many
Max occurrences. I feel that sympy.solve is faster when it doesn't fail,
so somehow sympy handles Piecewise functions pretty fast. I guess I'll
have to do some profiling on where I can still gain on performance.

Your comment also made me realize that in my systems there are some
underlying assumptions, as well, as inequalities. For the above system I
quickly checked and the corresponding assumption would be 1<=Min(x,y) and
Max(x,y)<=2. With that assumption the solution becomes unique. I guess
I'll need to handle these assumptions in my code, as well. (sympy.solve
dies with the same NaN comparison error even with the extended
assumptions)

Thanks!
Gábor

Chris Smith

unread,

Aug 9, 2025, 3:43:18 PMAug 9

to sympy

Expressions with Max/Min can already be converted to Piecewise:

>>> (1+Max(x,y)).rewrite(Piecewise)

Piecewise((x, x >= y), (y, True)) + 1

/c

Gábor Horváth

unread,

Aug 9, 2025, 5:39:32 PMAug 9

to sympy

> Expressions with Max/Min can already be converted to Piecewise:
> >>> (1+Max(x,y)).rewrite(Piecewise)
> Piecewise((x, x >= y), (y, True)) + 1

I'm sorry, I don't understand what you're implying. I understand that
Min/Max can be rewritten to Piecwise manually. My point was that after the
rewrite sympy.solve throws a NaN comparison error, and was wondering if
that is intended there, or if I can do anything about it.

Thanks,
Gábor

Chris Smith

unread,

Aug 11, 2025, 6:37:58 AMAug 11

to sympy

My guess is that solve is not intended to solve a system of inequalities so it doesn't know what to do with the relationals that are seen in the system of Piecewise. But if you know the structure you can build the appropriate product of args and reconstruct a solution:

solns= []
for p in product(*[P.args[1].args for P in new_sys]):
e, c = zip(*p); s = solve(e,(x,y),dict=True)
if s:solns.append((s,And(*c)))

>>> solns

[([{x: 0, y: 0}], True)]

/c

Oscar Benjamin

unread,

Aug 11, 2025, 7:49:56 AMAug 11

to sy...@googlegroups.com

On Sat, 9 Aug 2025 at 17:31, Gábor Horváth <hungabo...@gmail.com> wrote:
> > On Sat, 9 Aug 2025 at 13:26, Gábor Horváth <hungabo...@gmail.com> wrote:
> >>
> >> So my questions are:
> >> 1) Is the expected behaviour for sympy.solve for the Piecewise case to throw a NaN comparison exception, or is this a sign of some underlying issue in sympy.solve for Piecewise functions?
> >
> > No, that is just a bug.
>
> Ok, thanks. Should I report it in some official way somewhere (if yes,
> where and how), or has this email already taken care of the bugreport?
> Maybe I should open an issue on https://github.com/sympy/sympy/issues ?

Yes, please open an issue there.

I think that the expected behaviour here is to raise an error but it
should be an explicit error like "equation type not supported".

> >> 2) Is sympy.solve expected to rewrite Min and Max to Piecewise at some point in the future? I found some piece of code in solvers/solvers.py:944-957 which does some sort of rewriting for Abs, but didn't find any for Min or Max.
> >
> > It probably does make sense to rewrite Min and Max as Piecewise if
> > they contain the variables of interest.
>
> Ok, thanks. Same question as above: should I open some issue on
> https://github.com/sympy/sympy/issues or elsewhere, or this is already
> taken care of by this email?

Yes, please open an issue. It can be the same issue.

> >> 3) From some reading on the sympy docs page, I figured solveset is expected to be the future solver. Is there a direct way I can apply solveset on a system, rather than on one equation?
> >
> > This documentation needs to be changed. The solveset function is not
> > going to replace solve. It was originally intended as a replacement
> > for solve but was then misdesigned as something different that cannot
> > replace solve because it does not provide comparable functionality or
> > interface.
>
> Ah, I see. This was the page I found first:
> https://lidavidm.github.io/sympy/modules/solvers/solveset.html

That is an old version hosted from someone's outdated fork. Can you
open a separate sympy github issue about that? We can tag the owner of
the fork and ask them to remove this. It isn't possible to open an
issue in their repo as it has issues disabled I think:
https://github.com/lidavidm/sympy

> BTW, are there some old threads somewhere maybe where such discussions
> can be read? I would be curious what you mean by misdesign, and what were
> the discussion and decision points around solve an solveset.

Some will be here on the mailing list which you can search. Others
will be on GitHub.

What I mean by misdesign here is not that solveset is misdesigned in
general but rather that if the intention was to replace solve then it
was misdesigned for that purpose. It is not possible to replace solve
with solveset because they just don't do the same thing.

The purpose of a function like solve is that it turns an implicit
representation of the solutions (equations to be solved) into an
explicit representation which is a list of assignments like x = ..., y
= .... Its return value is therefore a list of dicts each of which
could be passed to subs.

There are some limitations in this explicit format which mean that it
cannot always be fully mathematically correct/complete. For example
your system of equations has the solution set:

{(a, a) for a in [2, oo)}

The best solve can do to represent that solution set is to return

[{x: y}]

but that loses the information that y >= 2 is a constraint on the
solutions. It could be possible to extend solve to also return the
condition like

[({x: y}), y >= 2)]

but that would be an incompatible change in its output format.

With solveset the intention is that the output is always
mathematically correct/complete but then by design it does not
necessarily return an explicit representation of the solutions which
contradicts the basic purpose of solve. There is a tradeoff here where
you have to choose between having explicit representations but with
some attached explicit or implicit conditions vs not having an
explicit representation at all.

The purpose of the solve function is to have the explicit
representations subject to the fact that they may be incomplete or not
always valid whereas the purpose of solveset is to be always complete
and valid and therefore not necessarily give an explicit
representation of the solutions. These are contradictory goals and so
in a basic way solveset cannot replace solve.

There are other reasons why solveset cannot replace solve. Another
reason is that solveset returns output in an unstructured form that is
not easily usable even when the solutions can be given explicitly.
There are some places in the sympy codebase that try to use solveset
and they always have awkward scaffolding code that basically tries to
turn the output of solveset into something that looks like the output
of solve. This is because the sets that solveset returns are only
really nice just to "look at". If you actually want to use the
solutions of the equation for something then you want an explicit
representation i.e. you basically want solve rather than solveset.

An even more obvious reason why solveset cannot replace solve is that
it only handles single univariate equations rather than systems of
equations. It has not even been designed in a way that it could be
extended to handle systems of equations in future.

> There is also mention in the docs about FiniteSet handling tuples
> corresponds to a finite subset of points in the multi-dimensional space.
> That suggests that solveset was intended to be able to solve multi-variate
> equations or even systems. Is that indeed the case, or only some special
> cases (like linsolve) can really handle multi-variate equations/systems?

The linsolve and nonlinsolve functions do handle systems of equations
and are useful. They also return a somewhat awkward output format and
do not really handle infinite solution sets in a coherent way which
undermines the benefit of not simply returning a list of dicts like
solve does. I am not sure that there is any particular benefit in
using either of these functions rather than just using solve. They are
basically just new implementations of the way that solve handles
systems and perhaps have cleaner code but otherwise do not
fundamentally increase capability or usefulness.

> > In general I think that what you should really want in a situation
> > like this is something like Mathematica's Reduce function rather than
> > Solve. Unfortunately SymPy has several versions of solve but nothing
> > like reduce.
>
> Hm, all are very good points. I have some very crude code to solve these
> systems exactly as you suggest, by splitting at the various Min and Max
> occurrences, solving the underlying linear system (with linsolve), and
> then check if the obtained solutions indeed solve the original system with
> Min and Max.
>
> That said, somehow my code is rather slow, especially when there are many
> Max occurrences. I feel that sympy.solve is faster when it doesn't fail,
> so somehow sympy handles Piecewise functions pretty fast. I guess I'll
> have to do some profiling on where I can still gain on performance.
>
> Your comment also made me realize that in my systems there are some
> underlying assumptions, as well, as inequalities. For the above system I
> quickly checked and the corresponding assumption would be 1<=Min(x,y) and
> Max(x,y)<=2. With that assumption the solution becomes unique. I guess
> I'll need to handle these assumptions in my code, as well. (sympy.solve
> dies with the same NaN comparison error even with the extended
> assumptions)

Okay, so your actual system here is:

x = Max(y, 1)
y = Max(2, x)
1<=Min(x,y)
Max(x,y)<=2

If I wanted to make a solver for systems like this then I would make
something that turns it into separate systems of inequalities that
don't involve Min and Max. For example x = Max(y, 1) can be separated
into two disjoint cases represented by pairs of inequalities:

(1)
x = y
y > 1

(2)
x = 1
y <= 1

You can do the same for each Max and Min. That will give you 8 (2^3)
systems of linear equations and inequalities. To solve those first
solve the linear equations and then substitute into the inequalities.
Then Fourier-Motzkin can solve the remaining system of inequalities.

If you have a large number of Max/Min then the approach I just
described might be inefficient but it depends how many of the Max/Min
are actually distinct. In the above case Max(x,y) and Min(x,y) can
both switch in the same condition x>y. That is why we get 2^3 cases
rather than 2^4. If your 2^n here is large then a better approach is
something that works by progressively solving and substituting e.g.
for case (1) above you substitute x = y into all remaining
inequalities before splitting any of the other Min/Max at which point
they simplify significantly.

If SymPy were to have a solver for this kind of system it would be one
of the inequality solvers but I don't think any of those handles a
multivariate system as is needed here. Ideally SymPy would have a
reduce function that would be able to handle exactly this kind of
system but it does not. Much of the effort that went into solveset
would have been better spent on a reduce function which would have
then been a good foundation for building a new version of solve.

--
Oscar

Gábor Horváth

unread,

Aug 20, 2025, 10:45:00 AMAug 20

to sy...@googlegroups.com

Hi Oscar,

Sorry for the late reply, I opened the two issues:

https://github.com/sympy/sympy/issues/28331
https://github.com/sympy/sympy/issues/28332

Not quite sure if the format and the text is fine, let me know if I need
to modify anything.

Thanks very much for the detailed explanation on solve vs solveset. I
found it very interesting, and I agree with the points. For my usecase I
definitely need to use solve, because I need the solutions afterwards.

I'm writing an implementation for my very special and restricted usecase
of systems, where I might need to convert between dicts and sets, mainly
to keep track of the potential infinite solutions:
{x[0]: x[2], x[1]: x[2], x[2]: x[2]} == {x[0]: x[1], x[1]: x[1], x[2]:
x[1]} is False for dicts, but they both represent the same FiniteSet if
the domain for x[1] and x[2] are the same. Is there some code in sympy
which converts between these in certain cases, or I can just go with my
(rather naive) approach with something like:

def subs_to_fset(subs:list[dict],
variables:Iterator[sympy.Symbol]) -> sympy.FiniteSet:
if any(len(d) != len(variables) for d in subs):
raise SubstitutionSizeNotMatcingNumberOfVariablesError
return sympy.FiniteSet(*[tuple(d[v] for v in variables) for d in
subs])

def fset_to_subs(fset:sympy.FiniteSet,
variables:Iterator[sympy.Symbol]) -> tuple[tuple]:
if any(len(p) != len(variables) for p in fset):
raise FiniteSetDimensionNotMatcingNumberOfVariablesError
return tuple(tuple(zip(variables, p)) for p in fset)

You also mention that there is no benefit using linsolve and nonlinsolve
vs solve. Do you mean that for the performance, as well, or only for
capability? To me performance may matter, and the awkward conversion
between sets and list[dicts] might be worth it in case linsolve is
significantly faster than solve. But if they essentially have the same
underlying code, then solve is far superior for its output.

Thanks for the pointers on how to handle systems w Min/Max. For now I
ended up using an extra domain parameter that restricts the potential
solution onto an interval (in my special usecase that's the situation),
and since all non-MinMax equations are linear, I can use some heuristics
to compute upper and lower bounds for each variable based on what
subsolutions / substitutions I found so far, and intersect that with the
domain of interval of possible solutions. In certain cases (apparently a
lot of the cases for my systems) this provides a solution to a lot of
variables, just like in the case for the minimal example I provided. This
way I can significantly reduce the number of equations having MinMax in
them.
For some big systems there still are many of those MinMax equations,
though (many meaning ~20, but that is already >1mil different cases to
solve for). So for those I'll try your suggested approach of recursively
solve one equation and substitute back, but that seems to be less
efficient for now that simply enumerating all potential 2**20-many
possibilities. I must have some silly mistake in my code, which I'll have
to check.

In any case, thank you very much for all the help!

Thanks,
Gábor

> --
> 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 sympy+un...@googlegroups.com.
> To view this discussion visit https://groups.google.com/d/msgid/sympy/CAHVvXxQ%3D0t4J9hqjhW43wYezd5%2BBwxKN8aeG6oBWmjkJU9hapQ%40mail.gmail.com.
>

Chris Smith

unread,

Aug 21, 2025, 8:48:10 PMAug 21

to sympy

Could you post an example of a large system?

/c

Gábor Horváth

unread,

Aug 24, 2025, 7:01:20 AMAug 24

to sympy

> 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.

Variables are v0 to v238, all are real, and they should obtain values from
Interval(0, 1). Their definition was
v = sympy.symbols('v:{n}'.format(n=n), real=True)

Not sure how to load the system from this form, though. It must be in the
docs somewhere, but don't recall reading it. :-(

System has 239 equations, all of them have lhs as one of the variables.
115 has even the rhs as a variable.
74 has Max on rhs, which always have two arguments, a number and another
symbol.
Rest are linear in the variables. In fact, they are a convex combination
of the variables and a rational number.

This is a pretty good representation of the systems I'm working with.

What I did so far is to solve the linear part, substitute back into the
equations with Max, then use the domain conditions and some easy
consequences of the equations to further restrict the potential domains
for the variables.
For example, v12 equals to Max(r, v44), thus domain of v12 is
always at least r, and the domain of it can be updated from Interval(0, 1)
to Interval(r, 1).
Another example is v71 equals to a*v141 + b*v178 + c*v45 + d, and here
using the (potentially restriced) domain conditions on the variables, I
can further restrict the domain of v71.

With this approach I managed to decrease the number of equations with Max
to 20, which then results in traversing through 2**20 cases, taking too
long. I probably need to do this recursively, applying the domain
shrinking heuristic at every recursive step, but I haven't yet managed
to do that.

Thanks,
Gábor

> To view this discussion visithttps://groups.google.com/d/msgid/sympy/1d1c80c9-8b37-48b6-ac1e-1286d2eb83e
> en%40googlegroups.com.
>
>

minmax_sys_example.txt

Oscar Benjamin

unread,

Aug 24, 2025, 10:30:38 AMAug 24

to sy...@googlegroups.com

On Sun, 24 Aug 2025 at 12:01, Gábor Horváth <hungabo...@gmail.com> 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

Chris Smith

unread,

Aug 25, 2025, 1:14:08 PMAug 25

to sympy

That system is helpful and you've got some good advice already. I would add that this sort of system can be reduced to a linear programming form and sympy has a function that handles your system:

```

from sympy import *
from sympy.solvers.simplex import lpmin

def lpsys(eqs):
eqs = list(eqs)
maxes = Tuple(*eqs).atoms(Max)
syms_max = symbols(f'm0:{len(maxes)}', real=True)
rep_m = dict(zip(maxes, syms_max))
eqs2 = [e.xreplace(rep_m) for e in eqs] # linear now after replacing Max with symbols

max_eq=[Eq(*z) for z in zip(syms_max, maxes)]
return eqs2, max_eq

eqs2, max_eq =lpsys(eqs) # your list of 239 equations is eqs
obj = sum(ms.lhs for ms in max_eq) # minimize the sum of m variables used to represent the Max functions

cons = list(eqs2)

# m_i >= each argument for Max's
for meq in max_eq:
for a in meq.rhs.args:
cons.append(meq.lhs >= a)

opt_val, sol = lpmin(obj, cons)

```

This gives 22.9893360923961 for the opt_val and solutions for the 269 equations (30 maxes and 239 others).

/c

Gábor Horváth

unread,

Aug 26, 2025, 4:25:01 PMAug 26

to sympy

Hi Oscar, Chris,

Thank you for both for your suggestions. I very much like the ideas from
both of you, and while I was doing something similar, this is more
clear-cut and understandable. I can follow it easier that my original
thoughts on it. :-)

Chris,
about the linear programming: This is wizard! I haven't thought about
rewriting to a linear programming problem, and even though it likely is
not equivalent, as long as there is only one solution satisfying all of
these conditions, it will find it.

Your solution needed a bit of tweaking: even though it is super fast, but
it's output is not correct, because it doesn't satisfy the Max equations,
it only satisfies that all of the args of Max is at most the Max variable.
I think this is due to the fact that it's not an equivalent rewrite.
However, when I add also that variables are supposed to be from the
Interval(0,1), then it spits out an opt_val of 23.138180070450467, and
that provides a solution that does satisfy all of the equations, even the
Max ones.

I did then the same computation with lpmax (with the conditions of being
in Interval(0,1), and that provided the same solution dictionary, and the
same opt_val value. And it was reasonably fast (within minutes!). Doing
these together can probably be used, and if opt from lpmin is the same as
opt from lpmax, together with the solution, then I can be reasonably
confident that that is the only solution. (The only way to be fully sure
is to compute for 30 linearly independent objective functions, but that
might be a bit of an overkill. :-) )

Oscar,
I extended your approach with some heuristics, namely knowing if v_i is
from Interval(a_i, b_i) then sum(c_iv_i) + sum(-d_jv_j) is from
Interval(sum a_iv_i-d_jb_j, b_iv_i -d_ja_j), and also that Max(v_i, v_j)
is then from the Interval(Max(a_i, a_j), Max(b_i, b_j)). Applying this
over and over solves 10 out of the 30 Max variables. Then I'll need to
start branching, and now I'm looking if I can apply the heuristic at every
recursive branching step, then maybe I can quickly discard those branches
where no solution lies, and arrive to the final solution.

My gut feeling is that this should be faster than the linear programming
approach, as this is a tailored solution to this specific problem, but
seeing my previous failed attempts at the brancinh and how fast the linear
programming is, I'm not so sure. I'll check both approaches for timing,
and report back.

Thanks again for all the help and pointers!

Gábor

> --
> 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 sympy+un...@googlegroups.com.

> To view this discussion visithttps://groups.google.com/d/msgid/sympy/403209cd-1bf6-4284-aef8-bc3a313a23e
> bn%40googlegroups.com.
>
>

Chris Smith

unread,

Aug 26, 2025, 4:53:16 PMAug 26

to sympy

Here is an AI translation based on Oscar's approach that only takes a few seconds. Does it do what you need?

```python
from __future__ import annotations
from typing import Dict, List, Tuple as _T, Iterable
from dataclasses import dataclass
from sympy import *
eqs=S('''[Eq(v0, v39), Eq(v1, v21), Eq(v2, v160), Eq(v3, 125*v136/1296 + 143*v157/648 + 13*v158/162 + 35*v196/324 + 125*v230/1296 + 7*v67/162 + 5*v96/81 + 155/648), Eq(v4, v46), Eq(v5, v21), Eq(v6, 7*v103/162 + 13*v13/162 + 143*v138/648 + 5*v165/81 + 125*v188/1296 + 35*v56/324 + 145/432), Eq(v7, v21), Eq(v8, 125*v136/1296 + 143*v157/648 + 13*v158/162 + 125*v238/1296 + 35*v64/324 + 7*v67/162 + 13*v70/162 + 5*v96/81 + 103/648), Eq(v9, 125*v107/1296 + 5*v114/81 + 35*v169/324 + 125*v176/1296 + 143*v234/648 + 13*v77/162 + 155/648), Eq(v10, v152), Eq(v11, 5*v150/81 + 35*v182/324 + 125*v23/1296 + 143*v48/648 + 13*v69/162 + 145/432), Eq(v12, Max(11973429987871294032791/15614510052024077690880, v44)), Eq(v13, v9), Eq(v14, 143*v131/648 + 35*v164/324 + 125*v174/1296 + 13*v179/162 + 125*v224/1296 + 155/648), Eq(v15, v97), Eq(v16, Max(22153944431869014019426428453312493667994869174799325454923686985538985161903624297370947248977440887249367238742791250601673384171676607719/32454635176053367993998877873873950990459011341470372284598283919422255066489869933707729610512961907754883836467874521812063002516520960000, v232)), Eq(v17, v87), Eq(v18, v236), Eq(v19, Max(8153526396622056593187661712693/9765476659969457979924424224000, v153)), Eq(v20, v97), Eq(v21, 35*v104/324 + 7*v115/162 + 125*v163/1296 + 125*v190/1296 + 5*v65/81 + 13*v66/162 + 143*v89/648 + 155/648), Eq(v22, Max(7275710867969158635412531309174807/8437371834213611694654702529536000, v180)), Eq(v23, v9), Eq(v24, Max(2614555075632615128064129937148739196765976986107121040972735374362471663685097832269563604123454373834497709933777759023191743438472863355594852611/3252479130413442679612535606086476366411081822411880393493248150349958821111651172916031511388192323765980310338178294488167832084589834141696000000, v229)), Eq(v25, 5*v168/81 + 7*v202/162 + 13*v36/162 + 35*v43/324 + 35*v68/1296 + 143*v7/648 + 125*v78/1296 + 145/432), Eq(v26, Max(3854097996225308316203175283663385/5973659258623237079815529390911488, v231)), Eq(v27, v35), Eq(v28, Max(7082014151066409483273977986654807/8437371834213611694654702529536000, v218)), Eq(v29, v3), Eq(v30, v156), Eq(v31, v32), Eq(v32, 5*v100/81 + 143*v162/648 + 35*v186/324 + 125*v40/1296 + 125*v49/1296 + 13*v52/162 + 13*v59/162 + 103/648), Eq(v33, v9), Eq(v34, v145), Eq(v35, 125*v126/1296 + 13*v173/162 + 125*v181/1296 + 13*v206/162 + 35*v211/324 + 143*v98/648 + 103/648), Eq(v36, v97), Eq(v37, v39), Eq(v38, Max(435458653906733290600307/519997654931407615526400, v106)), Eq(v39, 125*v159/1296 + 13*v19/162 + 125*v74/1296 + 35*v82/324 + 143*v91/648 + 13*v94/162 + 103/648), Eq(v40, Max(226542105817153888727102225905133917442673378023193639810034013/282361652451162721818510415100150396932752441247124690042880000, v25)), Eq(v41, v9), Eq(v42, Max(2394733571844680723947486046892229995824719990932886027/3293635682949831743806860531669548694292913046945792000, v213)), Eq(v43, Max(7275710867969158635412531309174807/8437371834213611694654702529536000, v180)), Eq(v44, 5*v150/81 + 125*v23/1296 + 143*v48/648 + 13*v69/162 + 575/1296), Eq(v45, v147), Eq(v46, 143*v134/648 + 125*v204/1296 + 35*v63/324 + 125*v75/1296 + 13*v86/162 + 103/648), Eq(v47, Max(17774053714151175033989655690633134825729039/24642541563136281050949853769599826662195200, v184)), Eq(v48, v97), Eq(v49, v39), Eq(v50, v46), Eq(v51, v147), Eq(v52, v35), Eq(v53, v151), Eq(v54, Max(279093664883899921882792279729481055727296597875364409018874143/330075542795919385459138332466610996414444404698606408499200000, v201)), Eq(v55, Max(279093664883899921882792279729481055727296597875364409018874143/330075542795919385459138332466610996414444404698606408499200000, v201)), Eq(v56, Max(17162041845107476095157951754621/19530953319938915959848848448000, v84)), Eq(v57, Max(134845448633481464957323497536579332465563743749277/172297325954688833637102978220838496248844582912000, v6)), Eq(v58, v147), Eq(v59, Max(7275710867969158635412531309174807/8437371834213611694654702529536000, v180)), Eq(v60, v183), Eq(v61, v229), Eq(v62, 35*v104/324 + 143*v118/648 + 7*v123/162 + 13*v125/162 + 35*v144/1296 + 125*v163/1296 + 125*v18/1296 + 5*v93/81 + 155/648), Eq(v63, Max(175636969241485433817388790704239855296988799222569480849525214868889482808102263013043142576463500794632154622270708701/294314475806011059781427635267647592615090539835558997573966535841876742686892500252978044232249079792621628328850227200, v160)), Eq(v64, Max(5061423529750895019224713210229530919488626594511426252665941900089074860550587856332801957246398173667897955595061549032234894743025570361086403/6440552733491965702203040804131636369130855093884911670283659703663284794280497372110953487897410542110852099679560979184490756603148186419200000, v183)), Eq(v65, v156), Eq(v66, v152), Eq(v67, v46), Eq(v68, v147), Eq(v69, v87), Eq(v70, Max(17162041845107476095157951754621/19530953319938915959848848448000, v84)), Eq(v71, 125*v141/1296 + 143*v178/648 + 13*v45/162 + 575/1296), Eq(v72, Max(253802678761938692393280149165758641400287017359935136027687/326807468114771668771424091551099996449944955147135057920000, v111)), Eq(v73, v236), Eq(v74, v156), Eq(v75, Max(3854097996225308316203175283663385/5973659258623237079815529390911488, v231)), Eq(v76, Max(7275710867969158635412531309174807/8437371834213611694654702529536000, v180)), Eq(v77, v35), Eq(v78, v183), Eq(v79, Max(8850883678872758268947173651427717/9863131426569152559723668466240000, v189)), Eq(v80, Max(8153526396622056593187661712693/9765476659969457979924424224000, v153)), Eq(v81, v97), Eq(v82, Max(2358584054497148792710148400261053499897881214516493192068924801916534408118507179932080824564044765382257088703238080044388691345123909081442699/3220276366745982851101520402065818184565427546942455835141829851831642397140248686055476743948705271055426049839780489592245378301574093209600000, v21)), Eq(v83, Max(3636449405556959850095431375189982793388139699327/5173750886263588479209023698635022807381206630400, v11)), Eq(v84, 7*v112/162 + 13*v139/162 + 5*v155/81 + 35*v198/1296 + 125*v5/1296 + 143*v60/648 + 575/1296), Eq(v85, Max(11973429987871294032791/15614510052024077690880, v44)), Eq(v86, Max(16191931/22535817, v71)), Eq(v87, 125*v12/1296 + 35*v214/324 + 143*v50/648 + 155/648), Eq(v88, 7*v112/162 + 13*v139/162 + 5*v155/81 + 35*v198/1296 + 35*v205/324 + 125*v5/1296 + 143*v60/648 + 145/432), Eq(v89, v236), Eq(v90, v8), Eq(v91, v152), Eq(v92, Max(8850883678872758268947173651427717/9863131426569152559723668466240000, v189)), Eq(v93, v152), Eq(v94, v46), Eq(v95, Max(8850883678872758268947173651427717/9863131426569152559723668466240000, v189)), Eq(v96, v122), Eq(v97, 35*v108/324 + 125*v175/1296 + 125*v74/1296 + 143*v91/648 + 13*v94/162 + 155/648), Eq(v98, v35), Eq(v99, v152), Eq(v100, v46), Eq(v101, v46), Eq(v102, Max(220238025723283192366108519240454879220580984891661788157350013/282361652451162721818510415100150396932752441247124690042880000, v154)), Eq(v103, v147), Eq(v104, Max(226542105817153888727102225905133917442673378023193639810034013/282361652451162721818510415100150396932752441247124690042880000, v25)), Eq(v105, v14), Eq(v106, 7*v103/162 + 13*v13/162 + 143*v138/648 + 5*v165/81 + 125*v188/1296 + 575/1296), Eq(v107, Max(435458653906733290600307/519997654931407615526400, v106)), Eq(v108, Max(253802678761938692393280149165758641400287017359935136027687/326807468114771668771424091551099996449944955147135057920000, v111)), Eq(v109, v87), Eq(v110, v35), Eq(v111, 5*v109/81 + 143*v194/648 + 125*v20/1296 + 13*v41/162 + 35*v43/324 + 7*v58/162 + 145/432), Eq(v112, v160), Eq(v113, Max(1258342213853621367032641137686368091858037627209580572330220504142983507324061383727708018843022293004047445501873262810166895090483369/1697773340450584222326787919746492518856403606479931590531402171972287877510455635787179828965942765628524996676494796077216101826560000, v3)), Eq(v114, v46), Eq(v115, v46), Eq(v116, Max(8153526396622056593187661712693/9765476659969457979924424224000, v153)), Eq(v117, v152), Eq(v118, v151), Eq(v119, v46), Eq(v120, 125*v105/1296 + 143*v133/648 + 5*v187/81 + 13*v209/162 + 575/1296), Eq(v121, 13*v191/162 + 125*v33/1296 + 143*v34/648 + 5*v51/81 + 575/1296), Eq(v122, 143*v110/648 + 35*v142/324 + 125*v195/1296 + 125*v227/1296 + 13*v80/162 + 103/648), Eq(v123, v156), Eq(v124, Max(2394733571844680723947486046892229995824719990932886027/3293635682949831743806860531669548694292913046945792000, v213)), Eq(v125, v39), Eq(v126, v223), Eq(v127, v151), Eq(v128, v46), Eq(v129, Max(16191931/22535817, v71)), Eq(v130, v35), Eq(v131, v152), Eq(v132, v151), Eq(v133, v9), Eq(v134, v46), Eq(v135, Max(11269204546127/14346413781285, v120)), Eq(v136, v32), Eq(v137, Max(7275710867969158635412531309174807/8437371834213611694654702529536000, v180)), Eq(v138, v3), Eq(v139, v3), Eq(v140, v207), Eq(v141, v221), Eq(v142, Max(2358584054497148792710148400261053499897881214516493192068924801916534408118507179932080824564044765382257088703238080044388691345123909081442699/3220276366745982851101520402065818184565427546942455835141829851831642397140248686055476743948705271055426049839780489592245378301574093209600000, v21)), Eq(v143, Max(226542105817153888727102225905133917442673378023193639810034013/282361652451162721818510415100150396932752441247124690042880000, v25)), Eq(v144, v46), Eq(v145, 5*v100/81 + 125*v137/1296 + 35*v143/324 + 143*v162/648 + 125*v49/1296 + 13*v52/162 + 155/648), Eq(v146, 125*v102/1296 + 143*v131/648 + 125*v174/1296 + 13*v179/162 + 35*v192/324 + 13*v212/162 + 103/648), Eq(v147, 125*v129/1296 + 143*v134/648 + 125*v204/1296 + 35*v26/324 + 155/648), Eq(v148, v156), Eq(v149, Max(279093664883899921882792279729481055727296597875364409018874143/330075542795919385459138332466610996414444404698606408499200000, v201)), Eq(v150, v147), Eq(v151, 7*v115/162 + 125*v167/1296 + 35*v171/324 + 125*v190/1296 + 5*v65/81 + 13*v66/162 + 13*v76/162 + 143*v89/648 + 103/648), Eq(v152, 35*v113/324 + 5*v114/81 + 125*v176/1296 + 13*v200/162 + 143*v234/648 + 125*v57/1296 + 13*v77/162 + 103/648), Eq(v153, 5*v109/81 + 143*v194/648 + 125*v20/1296 + 13*v41/162 + 7*v58/162 + 575/1296), Eq(v154, 143*v1/648 + 7*v161/162 + 5*v2/81 + 35*v22/324 + 125*v29/1296 + 13*v81/162 + 145/432), Eq(v155, v97), Eq(v156, 35*v203/324 + 143*v50/648 + 125*v83/1296 + 13*v85/162 + 103/648), Eq(v157, v8), Eq(v158, v152), Eq(v159, Max(253802678761938692393280149165758641400287017359935136027687/326807468114771668771424091551099996449944955147135057920000, v111)), Eq(v160, 125*v126/1296 + 125*v135/1296 + 13*v206/162 + 35*v47/324 + 143*v98/648 + 155/648), Eq(v161, v147), Eq(v162, v8), Eq(v163, Max(7275710867969158635412531309174807/8437371834213611694654702529536000, v180)), Eq(v164, Max(220238025723283192366108519240454879220580984891661788157350013/282361652451162721818510415100150396932752441247124690042880000, v154)), Eq(v165, v193), Eq(v166, v39), Eq(v167, Max(226542105817153888727102225905133917442673378023193639810034013/282361652451162721818510415100150396932752441247124690042880000, v25)), Eq(v168, v9), Eq(v169, Max(134845448633481464957323497536579332465563743749277/172297325954688833637102978220838496248844582912000, v6)), Eq(v170, v8), Eq(v171, Max(2105200769737816505172382620601987033203554288183000146797721880918236299783478163617023469627750000504871656527720958235556227417449073238133508201/2782318780868529183351713627384866911464529400558281841562540991982539031129174864751931906771681354191888107061570343007700006852560016533094400000, v62)), Eq(v172, v9), Eq(v173, Max(11269204546127/14346413781285, v120)), Eq(v174, v122), Eq(v175, Max(8153526396622056593187661712693/9765476659969457979924424224000, v153)), Eq(v176, v146), Eq(v177, Max(279093664883899921882792279729481055727296597875364409018874143/330075542795919385459138332466610996414444404698606408499200000, v201)), Eq(v178, v160), Eq(v179, v46), Eq(v180, 5*v168/81 + 7*v202/162 + 13*v36/162 + 35*v68/1296 + 143*v7/648 + 125*v78/1296 + 575/1296), Eq(v181, Max(17774053714151175033989655690633134825729039/24642541563136281050949853769599826662195200, v184)), Eq(v182, Max(8153526396622056593187661712693/9765476659969457979924424224000, v153)), Eq(v183, 13*v170/162 + 35*v177/324 + 125*v215/1296 + 7*v217/162 + 5*v237/81 + 35*v4/1296 + 125*v53/1296 + 143*v73/648 + 155/648), Eq(v184, 125*v105/1296 + 143*v133/648 + 5*v187/81 + 13*v209/162 + 35*v38/324 + 145/432), Eq(v185, Max(7839007286424827564645117347/9940275171668787978402662400, v121)), Eq(v186, Max(2105200769737816505172382620601987033203554288183000146797721880918236299783478163617023469627750000504871656527720958235556227417449073238133508201/2782318780868529183351713627384866911464529400558281841562540991982539031129174864751931906771681354191888107061570343007700006852560016533094400000, v62)), Eq(v187, v147), Eq(v188, v145), Eq(v189, 5*v15/81 + 35*v17/1296 + 7*v172/162 + 5*v220/324 + 13*v225/162 + 125*v235/1296 + 143*v61/648 + 575/1296), Eq(v190, v39), Eq(v191, v193), Eq(v192, Max(2045078160365334499396949348874892561929655603986003104366740054193468546098285839528846751540950480006451342487403164371269925361888435218161080201/2782318780868529183351713627384866911464529400558281841562540991982539031129174864751931906771681354191888107061570343007700006852560016533094400000, v216)), Eq(v193, 143*v110/648 + 125*v116/1296 + 125*v227/1296 + 35*v72/324 + 155/648), Eq(v194, v183), Eq(v195, Max(253802678761938692393280149165758641400287017359935136027687/326807468114771668771424091551099996449944955147135057920000, v111)), Eq(v196, Max(540603656426932533425269180020889323149790647994708176949439/653614936229543337542848183102199992899889910294270115840000, v88)), Eq(v197, Max(7839007286424827564645117347/9940275171668787978402662400, v121)), Eq(v198, v147), Eq(v199, Max(8850883678872758268947173651427717/9863131426569152559723668466240000, v189)), Eq(v200, Max(435458653906733290600307/519997654931407615526400, v106)), Eq(v201, 5*v15/81 + 35*v17/1296 + 7*v172/162 + 5*v220/324 + 13*v225/162 + 125*v235/1296 + 143*v61/648 + 35*v79/324 + 145/432), Eq(v202, v87), Eq(v203, Max(33555258956212545159265306236629059868117420803200021272171253837071413093415619201625551702480684238035240126327735818814535013747339/50980804700015500968416644345982125478843256558573011969114731928624900945625419111863289811598969388796638160975604215819955863552000, v208)), Eq(v204, v46), Eq(v205, Max(8850883678872758268947173651427717/9863131426569152559723668466240000, v189)), Eq(v206, v46), Eq(v207, 7*v10/162 + 5*v128/324 + 125*v132/1296 + 143*v140/648 + 35*v148/1296 + 35*v226/324 + 13*v228/162 + 5*v37/81 + 125*v54/1296 + 13*v92/162 + 103/648), Eq(v208, 143*v0/648 + 5*v101/81 + 35*v108/324 + 125*v117/1296 + 125*v175/1296 + 13*v30/162 + 155/648), Eq(v209, v160), Eq(v210, v122), Eq(v211, Max(164303259789428574676434357424189686228606168309404001224495213032351595510640234280468301450122934132127132665995958212872858503/242821238664149429242953838715430791222961708193172782203145156648304855136533203360117025851618319372030931646164857757106176000, v9)), Eq(v212, Max(7082014151066409483273977986654807/8437371834213611694654702529536000, v218)), Eq(v213, 13*v191/162 + 35*v28/324 + 125*v33/1296 + 143*v34/648 + 5*v51/81 + 145/432), Eq(v214, Max(3636449405556959850095431375189982793388139699327/5173750886263588479209023698635022807381206630400, v11)), Eq(v215, Max(8850883678872758268947173651427717/9863131426569152559723668466240000, v189)), Eq(v216, 143*v127/648 + 5*v130/81 + 125*v137/1296 + 35*v143/324 + 13*v166/162 + 7*v219/162 + 125*v90/1296 + 155/648), Eq(v217, v35), Eq(v218, 143*v1/648 + 7*v161/162 + 5*v2/81 + 125*v29/1296 + 13*v81/162 + 575/1296), Eq(v219, v46), Eq(v220, v147), Eq(v221, 35*v124/324 + 125*v197/1296 + 125*v233/1296 + 143*v27/648 + 155/648), Eq(v222, Max(11269204546127/14346413781285, v120)), Eq(v223, 35*v16/324 + 13*v185/162 + 125*v233/1296 + 143*v27/648 + 125*v42/1296 + 103/648), Eq(v224, Max(7082014151066409483273977986654807/8437371834213611694654702529536000, v218)), Eq(v225, v183), Eq(v226, Max(2614555075632615128064129937148739196765976986107121040972735374362471663685097832269563604123454373834497709933777759023191743438472863355594852611/3252479130413442679612535606086476366411081822411880393493248150349958821111651172916031511388192323765980310338178294488167832084589834141696000000, v229)), Eq(v227, v46), Eq(v228, v236), Eq(v229, 7*v10/162 + 5*v128/324 + 125*v132/1296 + 143*v140/648 + 35*v148/1296 + 35*v149/324 + 125*v199/1296 + 13*v228/162 + 5*v37/81 + 155/648), Eq(v230, Max(17162041845107476095157951754621/19530953319938915959848848448000, v84)), Eq(v231, 125*v141/1296 + 143*v178/648 + 35*v222/324 + 13*v45/162 + 145/432), Eq(v232, 5*v119/81 + 35*v164/324 + 13*v210/162 + 125*v224/1296 + 143*v31/648 + 125*v99/1296 + 155/648), Eq(v233, v46), Eq(v234, v152), Eq(v235, v21), Eq(v236, 13*v170/162 + 7*v217/162 + 5*v237/81 + 35*v24/324 + 35*v4/1296 + 125*v53/1296 + 125*v55/1296 + 143*v73/648 + 13*v95/162 + 103/648), Eq(v237, v39), Eq(v238, Max(540603656426932533425269180020889323149790647994708176949439/653614936229543337542848183102199992899889910294270115840000, v88))]''')

import sympy as sp
from sympy import Eq, Expr, Max, Tuple, symbols, linsolve, Matrix, S, Symbol

@dataclass(frozen=True)
class MaxSpec:
m: Symbol # the new epigraph variable (m_i)
args: _T[Expr, ...] # the Max arguments (already xreplaced in eqs2)
raw: Max # original Max expression (before replace)

def build_max_epigraph(eqs: Iterable[Eq]):

eqs = list(eqs)
maxes = Tuple(*eqs).atoms(Max)
syms_max = symbols(f'm0:{len(maxes)}', real=True)
rep_m = dict(zip(maxes, syms_max))

eqs2 = [e.xreplace(rep_m) for e in eqs] # linear (affine) now

# Record each Max’s argument list after replacement
max_specs: List[MaxSpec] = []
for raw, m in zip(maxes, syms_max):
args = Tuple(*[a.xreplace(rep_m) for a in raw.args])
max_specs.append(MaxSpec(m=m, args=args, raw=raw))

# Symbols present
syms_all = Tuple(*eqs2).free_symbols
return eqs2, max_specs, syms_max, syms_all, rep_m

def lower_bound_propagation(eqs2: List[Eq], max_specs: List[MaxSpec]) -> Dict[Symbol, sRational]:
# Start with bounds from Max(const, x, …) where a pure numeric const is present
L: Dict[Symbol, Rational] = {}
def upd(v, val):
if v not in L or val > L[v]:
L[v] = val

for ms in max_specs:
consts = [a for a in ms.args if a.is_Number]
if consts:
c = max(consts) # “epigraph” lower bound
upd(ms.m, nsimplify(c))

# Push through obvious affine equalities v = sum(a_i * w_i) + b with a_i >= 0
# (Very conservative: only pick up clearly nonnegative coefficients.)
# Repeat to a fixpoint (few iterations usually suffice)
unknowns = list(Tuple(*eqs2).free_symbols)
for _ in range(3): # small, safe number of passes; increase if you like
pushed = False
for e in eqs2:
lhs, rhs = e.lhs, e.rhs
if lhs.is_Symbol:
# collect nonnegative linear part and constant
rhs_lin = expand(rhs)
b = rhs_lin.as_independent(*unknowns, as_Add=True)[0]
rest = rhs_lin - b
lb = nsimplify(b) if b.is_Number else 0
ok = True
contrib = 0
for term in rest.as_ordered_terms():
coeff, sym = term.as_independent(*unknowns, as_Add=False)
if not sym.is_Symbol:
ok = False; break
coeff = nsimplify(coeff)
if coeff.is_Number and coeff >= 0 and sym in L:
contrib += coeff * L[sym]
elif coeff == 0:
continue
else:
ok = False; break
if ok:
newL = lb + contrib
if lhs not in L or newL > L[lhs]:
L[lhs] = newL; pushed = True
if not pushed:
break
return L

def solve_linear(eqs: List[Eq], unknowns: List[Symbol]) -> Dict[Symbol, Expr]:
# SymPy returns a FiniteSet with a Tuple solution (possibly parametric).
solset = linsolve(eqs, unknowns)
if not solset: # no solution
return {}
(sol_tuple,) = tuple(solset)
return {u: s for u, s in zip(unknowns, sol_tuple)}

def initial_policy(max_specs: List[MaxSpec], L: Dict[Symbol, Expr] | None = None):
policy = []
for ms in max_specs:
# Prefer a numeric const if present (freezes many Maxes immediately)
const_ix = None
best_val = None
for j, a in enumerate(ms.args):
if a.is_Number:
if best_val is None or a > best_val:
best_val, const_ix = a, j
if const_ix is not None:
policy.append(const_ix)
else:
policy.append(0)
return policy

def policy_iteration(eqs2: List[Eq], max_specs: List[MaxSpec],
syms_max: List[Symbol],
syms_all: Iterable[Symbol],
seed_numeric: Dict[Symbol, float] | None = None,
max_iters: int = 20,
verbose: bool = True):
# Unknown ordering: m’s first, then “don’t care”, then “specials”
# (Your earlier ordering works well; we’ll keep it simple: all symbols.)
unknowns = list(syms_max) + [s for s in syms_all if s not in syms_max]

# Precompute bound-based initial policy
L = lower_bound_propagation(eqs2, max_specs)
policy = initial_policy(max_specs, L)

for it in range(max_iters):
# Enforce chosen branches: add equalities m_i = args[policy[i]]
eqs_pol = list(eqs2)
for ms, j in zip(max_specs, policy):
eqs_pol.append(Eq(ms.m, ms.args[j]))

sol = solve_linear(eqs_pol, unknowns)
if not sol:
if verbose: print(f"[policy {it}] infeasible under current branch choices")
return None # infeasible policy (should be rare if system is well-posed)

# Evaluate to decide if any Max branch is violated
changed = False
for i, ms in enumerate(max_specs):
# Compute numeric-ish value of each arg under current solution
# Substitute symbolic solution first; then optionally a numeric seed
def val(expr):
v = expr.xreplace(sol)
if seed_numeric:
v = v.xreplace({k: nsimplify(vv) for k, vv in seed_numeric.items()})
try:
return N(v, 40)
except Exception:
return v
vals = [val(a) for a in ms.args]
mval = val(ms.m.xreplace(sol))
# pick the arg with the (numerically) largest value
# if ties, keep the current policy choice (stable)
try:
best_ix = max(range(len(vals)), key=lambda j: vals[j])
violated = (vals[best_ix] > mval + 1e-30) # tiny slack to avoid noise
except TypeError:
# Non-numeric comparison; fall back to symbolic max check
violated = False
best_ix = policy[i]
for j, a in enumerate(ms.args):
if simplify((a - ms.args[policy[i]]).xreplace(sol)) > 0:
best_ix = j; violated = True; break

if violated and best_ix != policy[i]:
policy[i] = best_ix
changed = True

if verbose:
print(f"[policy {it}] {'switch' if changed else 'fixed'}")
if not changed:
# Fixed point: return the full solved map (you can also project)
return sol

if verbose:
print("[policy] exceeded iteration cap.")
return None

# eqs = [...] # your big list of Eq(...) from the prompt

eqs2, max_specs, syms_max, syms_all, rep_m = build_max_epigraph(eqs)

# Optional: seed a few core symbols with plausible numerics to break ties consistently.
seed = {} # e.g., {Symbol('v21'): 0.73, Symbol('v111'): 0.78}

solution = policy_iteration(eqs2, max_specs, list(syms_max), syms_all, seed_numeric=seed, verbose=True)

if solution is None:
print("No fixed policy found (under current setup).")
else:
# You now have values/affine forms for ALL vars (including m_i).
# If you want only the original v-variables:
sol_orig = {s: e for s, e in solution.items() if s.name.startswith('v')}
print("Solved entries:", len(sol_orig))
```

/c

Gábor Horváth

unread,

Aug 27, 2025, 4:32:26 PMAug 27

to sympy

Hi Chris,

Hmmmmmm, very interesting. First I thought, "No way that an AI can write a
code that's faster than anyone else's here, including sympy's authors." I
even checked the solution with checksol, and the code below finds the
solution significantly faster than checksol verifies it. :-D
(BTW, that is because I guess checksol uses subs, and here xreplace
suffices for verification.)

It was nagging me all day (especially because I'm quite a non-believer in
AI), so I spent the significant part of the evening to understand in
detail what the code does. Surprisingly, half of the code is not even
used. The lower_bound_propagation is not used at all, which I think is
similar to my original idea, though I used both lower and upper bounds.
Further, initial_policy doesn't really matter I think, it could be set to
[0, 0, ... , 0] by default. This policy is supposed to keep track of which
argument is selected as currently being "maximal".

The key part is of course the policy_iteration. If I understand correctly,
what this does is that evaluates all m_i as per the current policy, then
solves the remaining linear equation, then evaluates all Max()_i
expressions, and checks if m_i == Max()_i or m_i < Max()_i. In the latter
case modifies the policy according to which argument Max()_i takes, and
iterates with the new policy.
(And btw, it does this with always evaluating the symbolic rationals to
float, which I'm not quite sure why on Earth to do, but anyway. :-)
Is it that much faster to calculate with floats instead of with
rationals?)

This feels like a simplex algorithm to me (so AI may have stolen your
linear programming idea, as well :-) ), and I'm very much intrigued by the
concept. Even though I'm not fully convinced that this always finds the
solution. I'm not even fully sure that there is always only one solution.
But I think in the case where each variable is a nonnegative linear
combination of the rest of the variables, then it may work, because it is
(likely) increasing the Max()_i expressions. I'm not fully convinced,
because if, say, m_i = Max(r_i, v_j), and v_j > r_i for the solution where
m_i = r_i is chosen doesn't necessarily mean that v_j > r_i will hold when
m_i = v_j is chosen, and we solve the system all over again. Or at least I
don't see a clear argument for it.

In any case, this is very useful and interesting insight (both into my
problem and into AI's "thinking", as well), thank you very much for this!

Thanks,
Gábor

> To view this discussion visithttps://groups.google.com/d/msgid/sympy/1c1fd924-ceda-45c9-8519-9a8d7c6f375
> 3n%40googlegroups.com.
>
>

Chris Smith

unread,

Aug 27, 2025, 7:18:36 PMAug 27

to sympy

I didn't review the AI solution very closely since we had a working LP solution. (And, BTW, I have found the quality of AI depends a lot on the details of the request and subsequent prompts and also how close a problem is to an already known type of problem.)

As far as " Is it that much faster to calculate with floats instead of with rationals?" It depends what you are doing. If computation with Rationals needs to keep making the denominator larger then floats can be much faster. When trying to compare expressions internally in SymPy we will use a low precision Float, if possible, to make an unambiguous decision about expression equality rather than maintain the exact Rationals.

One thing I did notice, however, was the use of nsimplify which made me a little nervous.

I asked the AI if the method is just simplex in disguise and it said,

```AI

Short answer: no, it isn’t simplex—it’s a policy-iteration / active-set method for a max-affine system.

What you ran there does this each iteration:

pick an “active” branch for every Max (a policy),
solve the resulting linear equalities (with linsolve / SCC),
check which Max constraints are violated and flip those branches,
repeat until no flips are needed.

That’s essentially Howard-style policy iteration (very close to an active-set method). It doesn’t keep a simplex tableau or walk polytope edges with entering/leaving variables; it repeatedly solves the linear system implied by a chosen active set and updates that set.

How it relates to simplex/LP:

If you add the LP objective min Σ m_i with epigraph constraints m_i ≥ each arg, then your problem becomes an LP. Your lpmin run solved exactly this (SymPy uses a simplex-type solver under the hood).
The policy loop can be viewed as an active-set solver for that LP’s KKT conditions: a policy = a guess of which inequalities are tight; the linear solve enforces equalities; the “flip on violation” is an active-set update. In nondegenerate cases it often lands at the same solution as the LP—but it’s not literally simplex, and without the objective you’re just finding a consistent fixed point (when unique, they coincide).

So: policy iteration ≠ simplex, but it’s simplex-adjacent (an active-set scheme). Your final lpmin approach is the canonical LP formulation; the earlier method is a good objective-free solver/certifier that often converges fast on monotone, max-affine systems.

```

/c

Chris Smith

unread,

Aug 28, 2025, 9:08:30 AMAug 28

to sympy

Regarding “lower_bound_propagation”: it was offered as a pre-pruner. Per AI, "It’s most useful when many nodes are Max(const, v) or when you want to freeze obvious branches before the loop. In your example it wasn’t necessary. Keep it in your toolbox for larger, noisier systems." It also observed that "“initial_policy” only affects how many flips you need and [prefers] the “constant branch first” because it freezes a lot of Max(const, ·) immediately, but [0,0,…] is fine."

/c

Oscar Benjamin

unread,

Aug 28, 2025, 10:05:34 AMAug 28

to sy...@googlegroups.com

On Wed, 27 Aug 2025 at 21:32, Gábor Horváth <hungabo...@gmail.com> wrote:
>
> (And btw, it does this with always evaluating the symbolic rationals to
> float, which I'm not quite sure why on Earth to do, but anyway. :-)
> Is it that much faster to calculate with floats instead of with
> rationals?)

It is much faster to use floats, especially machine precision floats
rather than sympy's arbitrary precision floats. This is precisely why
most numeric computing uses floats for non-integer numbers rather than
rationals even in situations where irrational numbers are not needed.

Whether or not using floats is acceptable depends on a number of
things. Usually floats are fine for solving non-degenerate systems for
N equations in N unknowns. What floats are not good for is
distinguishing whether or not a degenerate case holds or whether an
overdetermined linear system has solutions. Of course if you need the
outputs to be exact for further use then using floats is not
acceptable.

In your case you also have the Max conditions which are continuous in
all the variables so I think if every possible combination of
piecewise linear branches has a unique solution for the linear
equations then it is probably okay to use floats to find the
solutions.

If you are not sure whether or not solutions exist or if there can be
open sets of non-unique solutions then you would need to use exact
arithmetic to determine the nature of the solutions.

SymPy's nsolve function solves this system using floats and seems to
give the same answers as the AI-generated code:

In [80]: [solution[s].n() for s in syms][:10]
Out[80]:
[0.615451987540876, 0.744935782479303, 0.669554028724062,
0.50734056790375, 0.50734056790375, 0.779985609984258,
0.540175431864777, 0.802311871497269, 0.646942367048034,
0.792756582043413]

In [81]: list(nsolve(eqs, syms, [0.5]*len(syms)))[:10]
Out[81]:
[0.615451987540876, 0.744935782479303, 0.669554028724062,
0.50734056790375, 0.50734056790375, 0.779985609984258,
0.540175431864777, 0.802311871497269, 0.646942367048034,
0.792756582043413]

Here nsolve is quite slow though because mpmath's LU solver is slow
and the matrix here is 239x239.

If nsolve finds a solution then it means it is likely a locally unique
solution. You can do this yourself using Newton's method and with
numpy and machine precision floats instead of mpmath which is faster:

In [92]: F = Matrix([e.lhs - e.rhs for e in eqs])

In [93]: J = F.jacobian(syms)

In [94]: x0 = [0.5]*len(syms)

In [95]: xn = np.array(x0)

In [96]: Fl = lambdify([syms], F, 'numpy')

In [97]: Jl = lambdify([syms], J, 'numpy')

In [98]: xn1 = xn - np.linalg.solve(Jl(xn.flat), Fl(xn.flat)).flat

In [99]: print(np.sum(np.abs(xn1 - xn))); xn = xn1
46.006302733485825

In [100]: xn1 = xn - np.linalg.solve(Jl(xn.flat), Fl(xn.flat)).flat

In [101]: print(np.sum(np.abs(xn1 - xn))); xn = xn1
0.5634913620344619

In [102]: xn1 = xn - np.linalg.solve(Jl(xn.flat), Fl(xn.flat)).flat

In [103]: print(np.sum(np.abs(xn1 - xn))); xn = xn1
8.548717289613705e-15

In [104]: xn[:10]
Out[104]:
array([0.61545199, 0.74493578, 0.66955403, 0.50734057, 0.50734057,
0.77998561, 0.54017543, 0.80231187, 0.64694237, 0.79275658])

So Newton iteration has converged exactly after two iterations in this
case although a third iteration is used to confirm convergence. I
think what happened is that the initial conditions did not have the
right branches for all Maxes but after one iteration it has the right
branches and a single step of Newton solves the linear equations
exactly.

The slowest part above is computing the Jacobian but you can probably
generate that directly at the same time as generating the equations.
Of course if you are going to use numpy like this then it might be
better just to use numpy directly rather than having sympy set up the
inputs.

Once you have an approximate float solution like this you can use it
to check which branch of each Max holds. Then you can compute the
exact solution in rational numbers and verify that it holds exactly if
that is what is needed.

--
Oscar

Oscar Benjamin

unread,

Aug 28, 2025, 10:26:42 AMAug 28

to sy...@googlegroups.com

On Thu, 28 Aug 2025 at 15:05, Oscar Benjamin <oscar.j....@gmail.com> wrote:
>
> So Newton iteration has converged exactly after two iterations in this
> case although a third iteration is used to confirm convergence. I
> think what happened is that the initial conditions did not have the
> right branches for all Maxes but after one iteration it has the right
> branches and a single step of Newton solves the linear equations
> exactly.
>

> ...

>
> Once you have an approximate float solution like this you can use it
> to check which branch of each Max holds. Then you can compute the
> exact solution in rational numbers and verify that it holds exactly if
> that is what is needed.

Actually not even that is needed. Since Newton iteration converges
exactly for this case you can just do Newton with rational numbers:

In [2]: syms = list(ordered(Tuple(*eqs).atoms(Symbol)))

In [3]: F = Matrix([e.lhs - e.rhs for e in eqs])

In [4]: J = F.jacobian(syms)

In [5]: x0 = [Rational(1, 2)]*len(syms)

In [6]: xn = Matrix(x0)

In [7]: rep = dict(zip(syms, xn.flat()))

In [8]: xn1 = xn - J.xreplace(rep).solve(F.xreplace(rep))

In [9]: print(sum(map(abs, xn - xn1)).evalf(3)); xn = xn1
46.0

In [10]: rep = dict(zip(syms, xn.flat()))

In [11]: xn1 = xn - J.xreplace(rep).solve(F.xreplace(rep))

In [12]: print(sum(map(abs, xn - xn1)).evalf(3)); xn = xn1
0.563

In [13]: rep = dict(zip(syms, xn.flat()))

In [14]: xn1 = xn - J.xreplace(rep).solve(F.xreplace(rep))

In [15]: print(sum(map(abs, xn - xn1)).evalf(3)); xn = xn1
0

In [16]: print([N(e) for e in xn[:10]])
[0.615451987540876, 0.744935782479303, 0.669554028724062,
0.507340567903750, 0.507340567903750, 0.779985609984258,
0.540175431864777, 0.802311871497269, 0.646942367048034,
0.792756582043413]

In [17]: print(xn[0])
29111091937059880243204574405054181649424284485394821464265991381502076096979324693736677856978271553874982581817081789333748397161675933583/47300345967485318925827718364372480084013111896344922847317119820923097431308394914046771654831300577107337882031715520231586587481997312000

--
Oscar

Chris Smith

unread,

Aug 28, 2025, 1:40:40 PMAug 28

to sympy

That's a surprise to me, Oscar. I would not have expected Newton to work with a system involving Max like this. It is interesting, too, that SymPy knows the derivative of Max.

I think you have to watch out for (1) infeasible systems which would cycle and (2) singluar systems. For (1) AI gives as an example,

import sympy as sp
x, m = sp.symbols('x m', real=True)

def step(policy): # policy: 0 -> use m=x, 1 -> use m=-x
eqs = [sp.Eq(m, x if policy == 0 else -x), sp.Eq(x, 2*m + 1)]
sol = sp.solve(eqs, [x, m], dict=True)[0]
next_pol = 0 if sp.Max(sol[x], -sol[x]) == sol[x] else 1
return sol, next_pol

policy = 0
for _ in range(4):
sol, policy = step(policy)
print(f"sol={sol}, next_policy={policy}, m={sol[m]}, Max(x,-x)={sp.Max(sol[x], -sol[x])}")

which oscillates between (x,m) = (-1,-1) and (1/3,-1/3) whereas LP would indicate this as being infeasible.

As an example of (2) consider `x=max(y,0);y=max(x,0)` which has solution (t,t) with t>=0. LP would give a canonical solution like (0,0).

/c

Gábor Horváth

unread,

Sep 6, 2025, 3:17:54 PMSep 6

to sympy

Dear Oscar, Chris,

Thank you very much for all your help and contribution to this.
Finally, I found reasonably good introductory course slides on policy
iteration here:
https://web.stanford.edu/class/cme241/lecture_slides/Tour-DP.pdf

This refers to Markov Decision Processes (MDP) which are essentially what
I'm studying in the somewhat special case of gamma == 1 and no imediate
rewards. My system of equations are essentially the Bell equations for
the strongly connected components of the underlying digraph of an MDP.

Up until now I didn't even realize the similarity of my problem and MDPs.
Without your help I probably wouldn't have. Same with the policy
iteration. I cleaned up the AI's code to my liking and coding style. :-)
It is super-fast, even calling it on a system of 3358 equations (having
2401 Max equations, with 256 Max atoms) it stabilized the policy in 4
iterations in less than a minute on a Raspberry Pi 5!
That said, the linear programming method hung for a system of 2337
equations, having 719 Max equations, with only 45 different Max atoms! So
it seems policy iteration really provided the breakthrough here!

The only nagging thing was that I couldn't prove uniqueness of the
solution, and without that the algorithm is not complete. I was pretty
sure there is only one unique solution, but I couldn't prove it. Nor
could I find a counterexample in my special case. The literature I found
almost exclusively deals with the gamma < 1 situation, because all of the
proofs use the fact that the policy iteration is a gamma-contraction on an
appropriate Banach space, and thus Banach's fixed point theorem can be
used. This doesn't work at all in my case.

I also couldn't prove that the iteration really stops, and won't move
back-and-forth between two policies, but this didn't bother me much,
because in such a case I could simply throw out those policies and start
afresh. (And in practice it didn't seem to happen, anyway.)

Finally, after spending the significant time of the weekend thinking about
this, I finally managed to prove that in my case uniqueness of solutions
also holds for those kinds of systems which I'm interested in! This
completes the policy iteration solution, I don't need to assume
uniqueness, anymore, it is there.

Thank you again for all your help and contributions to this project!

One last question: checksol seems to check solutions by using subs instead
of xreplace. This makes sense of course, as in general this is how a
solution is supposed to be checked. Though would it not make sense to try
to do xreplace first? And if that satisfies the system already, then
return True, and only if not, then do the subs? I doubt the initial
xreplace would have a significant performance impact, especially compared
to subs, but in some (many?) cases it would be able to return the positive
answer significantly faster.

Let me know if you think it would make sense to open an issue about this.

Thanks for all the help!
Gábor

> --
> 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 sympy+un...@googlegroups.com.

> To view this discussion visithttps://groups.google.com/d/msgid/sympy/ba489ce3-cd0c-4281-b0d6-36a671c2a0a
> fn%40googlegroups.com.
>
>

Chris Smith

unread,

Sep 8, 2025, 11:31:15 AMSep 8

to sympy

I'm only observing from afar (so to speak) since I don't have a deep understanding of these systems, but in discussing with AI your observations, the following analysis was suggested which you might find useful to compare to the proofs that you have generated for yourself:

```python

Awesome — here’s a tiny, self-contained SymPy helper that scans your Bellman-style system, builds the **worst-case matrix** $\bar A$ (entrywise max of all action matrices for each row), and prints:

* max row sum $\| \bar A \|_{\infty}$ (a contraction certificate if < 1),
* max column sum $\| \bar A \|_{1}$ (a second bound),
* a **right Perron vector** $w \succ 0$ of $\bar A$ (via power iteration),
* the **weighted∞ operator norm** $\| \bar A \|_{\infty,w}=\max_i \frac{(\bar A w)_i}{w_i}$ (equals the spectral radius when $w$ is Perron),
* a one-line **certificate** that uniqueness & policy-iteration termination hold when any of those bounds is < 1.

You can paste this straight into your notebook.

```python
import sympy as sp
from sympy import Eq, Max, Matrix, Tuple

# ---------- parsing helpers ----------

def _ordered_vars(eqs):
# keep only "v*" unknowns; tweak if you use different names
vars_all = list(sp.ordered(Tuple(*eqs).atoms(sp.Symbol)))
return [v for v in vars_all if v.name.startswith('v')]

def _affine_coeffs(expr, vars_list):
"""Return (coeffs list, const) for an affine expr in vars_list."""
lin = sp.expand(expr)
coeffs = [sp.nsimplify(lin.diff(v)) for v in vars_list]
const = sp.nsimplify(lin - sum(c*v for c, v in zip(coeffs, vars_list)))
return coeffs, const

def build_barA(eqs, vars_list=None):
"""
Build worst-case matrix \bar A: for each equation v_i = Max_k{ a_{ik}^T v + b_{ik} }
(or plain affine), take entrywise maximum over actions k for row i.
Returns (\bar A, vars_list).
"""
if vars_list is None:
vars_list = _ordered_vars(eqs)
n = len(vars_list)
idx = {v:i for i, v in enumerate(vars_list)}

# collect per-row action coefficient lists
row_actions = [[] for _ in range(n)]

for E in eqs:
# normalize to v_i = RHS (ignore equations whose LHS isn't one of our vars)
if E.lhs in idx:
i = idx[E.lhs]
rhs = E.rhs
elif E.rhs in idx:
i = idx[E.rhs]
rhs = E.lhs
else:
continue

# treat Max arms as separate actions; plain affine is a single action
arms = rhs.args if isinstance(rhs, Max) else (rhs,)
for a in arms:
coeffs, _ = _affine_coeffs(a, vars_list)
row_actions[i].append(coeffs)

# build \bar A as entrywise max across actions for each row
Abar = sp.zeros(n, n)
for i in range(n):
if not row_actions[i]:
continue
for j in range(n):
# coeffs are rationals here, so Python max is exact
Abar[i, j] = max(sp.nsimplify(cj[j]) for cj in row_actions[i])
return Abar, vars_list

# ---------- diagnostics & certificate ----------

def analyze_bellman(eqs):
Abar, vars_list = build_barA(eqs)
n = Abar.shape[0]

# sanity: nonnegativity (warn if any negative coefficient sneaks in)
neg_entries = [(i, j, Abar[i, j]) for i in range(n) for j in range(n) if Abar[i, j] < 0]
if neg_entries:
print(f"WARNING: {len(neg_entries)} negative entries found in \\bar A; uniqueness test below may not apply.")

# row/column sum bounds (exact rationals)
row_sums = [sp.nsimplify(sum(Abar[i, j] for j in range(n))) for i in range(n)]
col_sums = [sp.nsimplify(sum(Abar[i, j] for i in range(n))) for j in range(n)]
rmax = max(row_sums) if row_sums else sp.Integer(0)
cmax = max(col_sums) if col_sums else sp.Integer(0)

# right Perron vector via power method (float, but only for a nicer bound/vector)
prec = 80
Af = sp.N(Abar, prec)
v = Matrix([1] * n)
for _ in range(100):
v1 = Af * v
s = max(abs(z) for z in v1) if n else 1
if s == 0:
break
v1 = v1 / s
if (v1 - v).norm() < 1e-30:
v = v1
break
v = v1

# weighted infinity norm with weight vector v (equals rho if v is Perron)
if n and all(vi != 0 for vi in v):
Av = Af * v
beta = max(float(Av[i] / v[i]) for i in range(n))
w = (v / sum(v)).applyfunc(lambda z: sp.N(z, 30)) # normalize to sum 1 for display
else:
beta = float(sp.N(rmax))
w = Matrix([sp.Rational(1, n) for _ in range(n)]) if n else Matrix([])

# report
print("=== Bellman-MDP structural scan ===")
print(f"vars: {n}, ||\\bar A||_∞ (max row sum): {rmax}")
print(f" ||\\bar A||_1 (max col sum): {cmax}")
print(f"weighted-∞ operator norm with Perron-like w: ~ {beta:.12f}")
if n <= 12:
print("w (weights, sum=1):", list(w))

# certificate
cert = None
if rmax.is_Rational and rmax < 1:
cert = f"Contraction in ∞-norm: ||\\bar A||_∞ = {rmax} < 1 ⇒ unique fixed point & policy iteration terminates."
elif cmax.is_Rational and cmax < 1:
cert = f"Contraction in 1-norm: ||\\bar A||_1 = {cmax} < 1 ⇒ unique fixed point & policy iteration terminates."
elif beta < 1 - 1e-12:
cert = f"Spectral certificate: ρ(\\bar A) ≤ {beta:.12f} < 1 via weighted ∞-norm with w ⇒ unique fixed point & PI termination."
else:
cert = "No contraction certificate found (<1) — system may still be fine, but uniqueness/termination not certified by these bounds."
print("Certificate:", cert)

return {
"Abar": Abar,
"row_sums": row_sums,
"col_sums": col_sums,
"rmax": rmax,
"cmax": cmax,
"beta_weighted_inf": beta,
"w": w,
}
```

### How to use

```python
# eqs = [...] # your original list of Eq(...) with Max on the RHS (or LHS)
report = analyze_bellman(eqs)

# Quick check: if report["rmax"] < 1 (or beta_weighted_inf < 1)
# you’ve got a clean uniqueness + termination certificate.
```

### What this is certifying

* $\|\bar A\|_\infty = \max_i \sum_j \bar A_{ij} < 1$ ⇒ $T(v)=\max_a\{A^{(a)}v+b^{(a)}\}$ is a contraction in the $\infty$-norm, hence a **unique** fixed point and **policy iteration** cannot cycle.
* Even if the raw max row sum isn’t <1, the **weighted**∞ norm with a positive weight vector $w$ can drop the operator norm below 1. Choosing $w$ as (an approximation to) the **right Perron vector** of $\bar A$ makes

$$
\| \bar A \|_{\infty,w}=\max_i \frac{(\bar A w)_i}{w_i}\approx \rho(\bar A),
$$

so if that value is <1 you again have a clean certificate.

/c

Gábor Horváth

unread,

Sep 8, 2025, 4:40:01 PMSep 8

to sympy

Hi Chris,

This is the output for that example-equation:

=== Bellman-MDP structural scan ===

vars: 239, ||\bar A||_∞ (max row sum): 1
||\bar A||_1 (max col sum): 18
weighted-∞ operator norm with Perron-like w: ~ 1.000000000000
Certificate: No contraction certificate found (<1) — system may still be

fine, but uniqueness/termination not certified by these bounds.

When I thought about the proof, I thought long and hard if there is
another norm for which my system could have <1 norm, which would make
Banach's fixed point theorem usable. But it looked to me that any natural
norm would just end up using that gamma thing (which is 1 for me), so I'll
need a different way. Apparently, AI's certificate corroborates it, as
well. :-D

One thing I managed to prove much earlier, is that for each policy when
one solves the corresponding system of linear equations, that solution is
unique. That is because the linear system (that comes from these Markov
Decision Processes in my case) for each policy has an underlying matrix
that is weakly chained diagonally dominant, and such non-singular. This
for some weird systems (like [Eq(x,Max(1,y)), Eq(y,Max(2,x))]) needs an
auxiliary information on the domain (like x,y in Interval(1,2)), which is
coming from how my systems are generated.
(Side-note: for this system for policy where we choose the y and x
branch in both Max gives the system x=y and y=x, which is a singular
system without any auxiliary bound on the domain. However, applying x,y<=2
we have Max(2,x) == 2, thus the second equation can actually be rewritten
to y=2, and now the system becomes non-sungular no matter which branch we
choose for the first equation.)

Then it was also reasonably easy to prove that there is at least one
solution to the overall system, using those techniques with the operators
from here:
https://web.stanford.edu/class/cme241/lecture_slides/Tour-DP.pdf
I even managed to prove that there is one solution that is maximal in the
sense that in all the variables it is higher than any other solution.

What I struggled to prove was that if there are multiple policies
providing solutions, then those solutions actually have to be the same.
That needed a combinatorial way of playing around with these operators,
but didn't need any Banach's fixpoint theorem, only monotonicity, and
boundedness of the possible solutions I'm looking for. The latter is not
encoded in the system itself, but is an auxiliary information I am
passing to the equations (the domain).

I'm not surprised that none of the norms worked out. These cases don't
seem like contractions to me, but they had a feel of monotonicity, which
I ended up using in my proof.

That said, these observations made me think of one thing: as AI probably
knows all this stuff from the internet, it would be super helpful if it
gave the exact references to all this stuff and where it gets them. Like
it could have been helpful if it gave a reference to policy iteration or
active sets, or it could really give a few links on where it gets all
these norms from. These are actually natural norms, but the way it creates
the Abar matrix from the arguments of the Maxes is interesting, I wonder
where it got that. It might just be a simple differentiation output,
though, that's how it got the coefficients of the linear part, as well
(instead of doing something like expr.collect(v).coeff(v), which I think
would make more sense...)

Thanks,
Gábor

> To view this discussion visithttps://groups.google.com/d/msgid/sympy/44b88397-79ef-45ef-a715-65e181f4e2f
> 8n%40googlegroups.com.
>
>

Reply all

Reply to author

Forward