mpmath
mpmath copied to clipboard
mpmath
findroot bisect solver should dynamically set maxsteps
The default maxsteps for the bisect solver is not high enough to find roots for higher precision. Consider the following:
In [37]: import mpmath as mp
In [38]: mp.mp.dps = 15
In [39]: mpmath.findroot(lambda x: mp.cos(x) - x, [0, 1], solver='bisect')
Out[39]: mpf('0.73908513321516064')
Works fine, but
In [40]: mp.mp.dps = 100
In [41]: mpmath.findroot(lambda x: mp.cos(x) - x, [0, 1], solver='bisect')
---------------------------------------------------------------------------
ValueError Traceback (most recent call last)
<ipython-input-41-3f13b0529a9a> in <module>()
----> 1 mpmath.findroot(lambda x: mp.cos(x) - x, [0, 1], solver='bisect')
/Users/aaronmeurer/anaconda3/lib/python3.5/site-packages/mpmath/calculus/optimization.py in findroot(ctx, f, x0, solver, tol, verbose, verify, **kwargs)
973 '(%g > %g)\n'
974 'Try another starting point or tweak arguments.'
--> 975 % (norm(f(*xl))**2, tol))
976 return x
977 finally:
ValueError: Could not find root within given tolerance. (1.00455e-63 > 1.39525e-104)
Try another starting point or tweak arguments.
With a little bisection, I find that the maxsteps needs to be at least 172 for 100 digits of precision:
In [42]: mpmath.findroot(lambda x: mp.cos(x) - x, [0, 1], solver='bisect', maxsteps=171)
---------------------------------------------------------------------------
ValueError Traceback (most recent call last)
<ipython-input-42-196333f2d7bd> in <module>()
----> 1 mpmath.findroot(lambda x: mp.cos(x) - x, [0, 1], solver='bisect', maxsteps=171)
/Users/aaronmeurer/anaconda3/lib/python3.5/site-packages/mpmath/calculus/optimization.py in findroot(ctx, f, x0, solver, tol, verbose, verify, **kwargs)
973 '(%g > %g)\n'
974 'Try another starting point or tweak arguments.'
--> 975 % (norm(f(*xl))**2, tol))
976 return x
977 finally:
ValueError: Could not find root within given tolerance. (5.31455e-104 > 1.39525e-104)
Try another starting point or tweak arguments.
In [43]: mpmath.findroot(lambda x: mp.cos(x) - x, [0, 1], solver='bisect', maxsteps=172)
Out[43]: mpf('0.7390851332151606416553120876738734040134117589007575191876080010538405475879238025538562218703542515087')
For 200 digits I need 337 steps
In [44]: mp.mp.dps = 200
In [45]: mpmath.findroot(lambda x: mp.cos(x) - x, [0, 1], solver='bisect', maxsteps=336)
---------------------------------------------------------------------------
ValueError Traceback (most recent call last)
<ipython-input-45-0e1552ec2280> in <module>()
----> 1 mpmath.findroot(lambda x: mp.cos(x) - x, [0, 1], solver='bisect', maxsteps=336)
/Users/aaronmeurer/anaconda3/lib/python3.5/site-packages/mpmath/calculus/optimization.py in findroot(ctx, f, x0, solver, tol, verbose, verify, **kwargs)
973 '(%g > %g)\n'
974 'Try another starting point or tweak arguments.'
--> 975 % (norm(f(*xl))**2, tol))
976 return x
977 finally:
ValueError: Could not find root within given tolerance. (3.91903e-204 > 1.59475e-204)
Try another starting point or tweak arguments.
In [46]: mpmath.findroot(lambda x: mp.cos(x) - x, [0, 1], solver='bisect', maxsteps=337)
Out[46]: mpf('0.73908513321516064165531208767387340401341175890075746496568063577328465488354759459937610693176653184919819275735830418435568063419760355316796482170968483234993805142330598738012953875703419358306567112')
And 504 for 300 digits:
In [62]: mp.mp.dps = 300
In [63]: mpmath.findroot(lambda x: mp.cos(x) - x, [0, 1], solver='bisect', maxsteps=503)
---------------------------------------------------------------------------
ValueError Traceback (most recent call last)
<ipython-input-63-2f854689d955> in <module>()
----> 1 mpmath.findroot(lambda x: mp.cos(x) - x, [0, 1], solver='bisect', maxsteps=503)
/Users/aaronmeurer/anaconda3/lib/python3.5/site-packages/mpmath/calculus/optimization.py in findroot(ctx, f, x0, solver, tol, verbose, verify, **kwargs)
973 '(%g > %g)\n'
974 'Try another starting point or tweak arguments.'
--> 975 % (norm(f(*xl))**2, tol))
976 return x
977 finally:
ValueError: Could not find root within given tolerance. (2.30543e-304 > 1.82278e-304)
Try another starting point or tweak arguments.
In [64]: mpmath.findroot(lambda x: mp.cos(x) - x, [0, 1], solver='bisect', maxsteps=504)
Out[64]: mpf('0.73908513321516064165531208767387340401341175890075746496568063577328465488354759459937610693176653184980124664398716302771490369130842031578044057462077916282162797578542473636933798179130401805184490463795874555827715563044280142803837427540806357309007481391227727881966551733721964594436217179273782')
It appears to be a roughly linear relationship (bisection should converge exponentially, so this makes sense).
maxsteps=1.7*dps seems to do the trick (maybe it can be a bit tighter; I'm guessing the theoretical value is somewhere around log2(10)/2*dps=1.66*dps).
Related issue (maybe also related to https://github.com/fredrik-johansson/mpmath/issues/338), the error message here:
In [100]: mpmath.mp.dps = 1000
In [101]: mpmath.findroot(lambda x: mp.cos(x) - x, [0, 1], solver='bisect', maxsteps=int(1.6*1000))
---------------------------------------------------------------------------
ValueError Traceback (most recent call last)
<ipython-input-101-f7d2e2402e0a> in <module>()
----> 1 mpmath.findroot(lambda x: mp.cos(x) - x, [0, 1], solver='bisect', maxsteps=int(1.6*1000))
/Users/aaronmeurer/anaconda3/lib/python3.5/site-packages/mpmath/calculus/optimization.py in findroot(ctx, f, x0, solver, tol, verbose, verify, **kwargs)
973 '(%g > %g)\n'
974 'Try another starting point or tweak arguments.'
--> 975 % (norm(f(*xl))**2, tol))
976 return x
977 finally:
ValueError: Could not find root within given tolerance. (0 > 0)
Try another starting point or tweak arguments.
In [102]: mpmath.findroot(lambda x: mp.cos(x) - x, [0, 1], solver='bisect', maxsteps=int(1.7*1000))
Out[102]: mpf('0.7390851332151606416553120876738734040134117589007574649656806357732846548835475945993761069317665318498012466439871630277149036913084203157804405746207786885249038915392894388450952348013356312767722315809563537765724512043734199364335125384097800343406467004794021434780802718018837711361382042066316335037277991696731223230061388658203621770810997897062684240588094898683261860600485898958548725736764015075227608180391459518101628159120096461646067544051326415171064466281109360825848783713839555561751414947236484445905951725761941148124675640177779314940116256134645924226342365252988830797652100338640082594708891018634196239732723884319303214358459663243138212618441901510120562257664282369345426395275701837074195755827548189527226772711639712666269485383630038160237483500476900830750278028872088259473095780628227615380400702052684894648023683076752844039751753554952926602731285027770932877599438100266338506241996659896926245284831688941621543668733559784374419780341206614185258016707072152')
says (0 > 0).
The 0 > 0 is because the error message uses "%g". That's the only place in the code that uses that. Is there a better alternative?
The 0 > 0 is because the error message uses "%g". That's the only place in the code that uses that. Is there a better alternative?
@asmeurer, this seems to be fixed by #341 and is present in 1.3.0.
It appears to be a roughly linear relationship (bisection should converge exponentially, so this makes sense).
- There a other solvers.
- bisection shouldn't converge exponentially, ~~rather the newton method~~.
BTW, 'anderson' solver works fine on your examples:
>>> mp.dps = 100
>>> mpmath.findroot(lambda x: mp.cos(x) - x, [0, 1], solver='anderson')
mpf('0.7390851332151606416553120876738734040134117589007574649656806357732846548835475945993761069317665318498')
>>> mp.dps = 200
>>> mpmath.findroot(lambda x: mp.cos(x) - x, [0, 1], solver='anderson')
mpf('0.73908513321516064165531208767387340401341175890075746496568063577328465488354759459937610693176653184980124664398716302771490369130842031578044057462077868852490389153928943884509523480133563127677223158')
>>> mp.dps = 300
>>> mpmath.findroot(lambda x: mp.cos(x) - x, [0, 1], solver='anderson')
mpf('0.739085133215160641655312087673873404013411758900757464965680635773284654883547594599376106931766531849801246643987163027714903691308420315780440574620778688524903891539289438845095234801335631276772231580956353776572451204373419936433512538409780034340646700479402143478080271801883771136138204206631634')
So, I'm not sure if should make API more complex (i.e. set maxsteps=None and then change it dynamically for some solvers. I think there are better alternatives. bisect is too slow, anyway.
The fact that there are other solvers is irrelevant to this issue. Of course other solvers work for my toy example. In general, bisection may be more appropriate for some types of equations than other solvers as it has nice properties that other solvers do not.
bisection shouldn't converge exponentially
Not sure why you think that. The convergence rate is trivial to compute https://en.wikipedia.org/wiki/Bisection_method#Analysis.
So, I'm not sure if should make API more complex (i.e. set maxsteps=None and then change it dynamically for some solvers. I think there are better alternatives. bisect is too slow, anyway.
mpmath should automatically do the right thing when the precision is set higher than the default.
Not sure why you think that.
I was thinking of https://en.wikipedia.org/wiki/Rate_of_convergence
maxsteps=1.7dps seems to do the trick (maybe it can be a bit tighter; I'm guessing the theoretical value is somewhere around log2(10)/2dps=1.66*dps).
Keep in mind, it will depend on the interval (it's not scaled to [0,1] in the Bisect solver).
I think the correct theoretical value (but only for default norm function) is something like: ceil(log2(|b-a|/(tol*max(1, max(|x| in [a, b])))).
mpmath should automatically do the right thing when the precision is set higher than the default.
Yep, probably you are right in this case. Solvers have different values for maxsteps properties (now class-wide), maxsteps value could be overridden from kwargs. So, you could just add a dynamic property instead. No API break, PR welcomed.