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
On Thursday, August 28, 2025 at 9:26:42 AM UTC-5 Oscar wrote:
> 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
>
--
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+unsubscr...@googlegroups.com.
To view this discussion visit
https://groups.google.com/d/msgid/sympy/ba489ce3-cd0c-4281-b0d6-36a671c2a0afn%40googlegroups.com.
