[SciPy-Dev] optimize.fsolve endless loop with nan
Skipper Seabold
jsseabold at gmail.com
Mon Mar 18 11:40:49 EDT 2013
On Mon, Mar 18, 2013 at 11:29 AM, <josef.pktd at gmail.com> wrote:
> On Mon, Mar 18, 2013 at 11:14 AM, <josef.pktd at gmail.com> wrote:
> > On Wed, Mar 13, 2013 at 1:24 PM, <josef.pktd at gmail.com> wrote:
> >> On Wed, Mar 13, 2013 at 12:59 PM, <josef.pktd at gmail.com> wrote:
> >>> On Wed, Mar 13, 2013 at 12:17 PM, <josef.pktd at gmail.com> wrote:
> >>>> On Wed, Mar 13, 2013 at 10:51 AM, <josef.pktd at gmail.com> wrote:
> >>>>> preliminary question, I didn't have time yet to look closely
> >>>>>
> >>>>>>>> scipy.__version__
> >>>>> '0.9.0'
> >>>>>
> >>>>> I have a problem where fsolve goes into a range where the values are
> >>>>> nan. After that it goes into an endless loop, as far as I can tell.
> >>>>>
> >>>>> Something like this has been fixed for optimize.fmin_bfgs. Was there
> a
> >>>>> fix for this also for fsolve, since 0.9.0?
> >>>>>
> >>>>>
> >>>>> (The weirder story: I rearranged some test, and made unfortunately
> >>>>> also some other changes, and now when I run nosetests it never
> >>>>> returns. Ctrl+C kills nosetests, but leaves a python process running.
> >>>>> I have no clue why the test sequence should matter.)
> >>>>
> >>>> I had left the starting values in a module global even after I started
> >>>> to adjust them in one of the cases.
> >>>>
> >>>> The starting value for fsolve was in a range where the curvature is
> >>>> very flat, and fsolve made huge steps into the invalid range. After
> >>>> getting nans, it went AWOL.
> >>>>
> >>>> If I return np.inf as soon as I get a nan, then fsolve seems to stop
> >>>> right away. Is there a way to induce fsolve to stay out of the nan
> >>>> zone, for example returning something else than inf?
> >>>>
> >>>> I don't want to find a very specific solution, because I'm throwing
> >>>> lot's of different cases at the same generic method.
> >>>
> >>> same result with python 2.7, scipy version 0.11.0b1
> >>>
> >>>>"C:\Programs\Python27\python.exe" fsolve_endless_nan.py
> >>> scipy version 0.11.0b1
> >>> args 0.3 [100] 0.05
> >>> args 0.3 [ 100.] 0.05
> >>> args 0.3 [ 100.] 0.05
> >>> args 0.3 [ 100.00000149] 0.05
> >>> args 0.3 [-132.75434239] 0.05
> >>> fsolve_endless_nan.py:36: RuntimeWarning: invalid value encountered in
> sqrt
> >>> pow_ = stats.nct._sf(crit_upp, df, d*np.sqrt(nobs))
> >>> fsolve_endless_nan.py:39: RuntimeWarning: invalid value encountered in
> sqrt
> >>> pow_ += stats.nct._cdf(crit_low, df, d*np.sqrt(nobs))
> >>> args 0.3 [ nan] 0.05
> >>>
> >>>
> >>> standalone test case (from my power branch)
> >>>
> >>> Don't run in an interpreter (session) that you want to keep alive!
> >>> And open TaskManager if you are on Windows :)
> >>>
> >>> ------------
> >>> # -*- coding: utf-8 -*-
> >>> """Warning: endless loop in runaway process, requires hard kill of
> process
> >>>
> >>> Created on Wed Mar 13 12:44:15 2013
> >>>
> >>> Author: Josef Perktold
> >>> """
> >>>
> >>> import numpy as np
> >>> from scipy import stats, optimize
> >>>
> >>> import scipy
> >>> print "scipy version", scipy.__version__
> >>>
> >>>
> >>> def ttest_power(effect_size, nobs, alpha, df=None,
> alternative='two-sided'):
> >>> '''Calculate power of a ttest
> >>> '''
> >>> print 'args', effect_size, nobs, alpha
> >>> d = effect_size
> >>> if df is None:
> >>> df = nobs - 1
> >>>
> >>> if alternative in ['two-sided', '2s']:
> >>> alpha_ = alpha / 2. #no inplace changes, doesn't work
> >>> elif alternative in ['smaller', 'larger']:
> >>> alpha_ = alpha
> >>> else:
> >>> raise ValueError("alternative has to be 'two-sided', 'larger'
> " +
> >>> "or 'smaller'")
> >>>
> >>> pow_ = 0
> >>> if alternative in ['two-sided', '2s', 'larger']:
> >>> crit_upp = stats.t.isf(alpha_, df)
> >>> # use private methods, generic methods return nan with
> negative d
> >>> pow_ = stats.nct._sf(crit_upp, df, d*np.sqrt(nobs))
> >>> if alternative in ['two-sided', '2s', 'smaller']:
> >>> crit_low = stats.t.ppf(alpha_, df)
> >>> pow_ += stats.nct._cdf(crit_low, df, d*np.sqrt(nobs))
> >>> return pow_
> >>>
> >>> func = lambda nobs, *args: ttest_power(args[0], nobs, args[1])
> >>>
> >>> print optimize.fsolve(func, 100, args=(0.3, 0.05))
> >>> ------------
> >>
> >> correction to get a solution that would make sense (last two lines)
> >>
> >> func = lambda nobs, *args: ttest_power(args[0], nobs, args[1]) - args[2]
> >>
> >> print optimize.fsolve(func, 10, args=(0.76638635, 0.1, 0.8))
> >> #converges to 12
> >> print optimize.fsolve(func, 100, args=(0.76638635, 0.1, 0.8))
> #runaway
> >>
> >> Josef
> >> "Lost in Translation"
> >
> > (continuing my monologue)
> >
> > more problems with fsolve, looks like on Ubuntu
> >
> > Skipper uses scipy master, but python-xy daily testing uses
> > python-scipy (0.10.1+dfsg1-4)
> > TravisCI, which also runs Ubuntu doesn't have any problems (using
> > python-scipy amd64 0.9.0+dfsg1-1ubuntu2 )
> > https://travis-ci.org/statsmodels/statsmodels/jobs/5585028
> >
> > (just realized: I've tested on Windows so far only with scipy 0.9)
>
> I should have done the testing across versions:
>
> I'm getting the same failures on Windows with
> >>> import scipy
> >>> scipy.__version__
> '0.11.0b1'
>
> but not with scipy 0.9
>
> something has changed in the recent scipy versions that makes fsolve
> more fragile
>
>
Another data point if it helps. Using 'lm' of the new optimize.root is ok,
but all of the other method also fail on this - even when starting very
close the actual root. Most reach maximum iterations.
[~/statsmodels/statsmodels-skipper/statsmodels/stats/tests/]
[80]: optimize.root(func, .15, method="lm")
[80]:
status: 2
qtf: array([ 0.])
cov_x: array([[ 0.08608474]])
ipvt: array([1], dtype=int32)
success: True
nfev: 16
fun: array([ 0.])
x: array([ 0.05])
message: 'The relative error between two consecutive iterates is at most
0.000000'
fjac: array([[-3.40829292]])
[~/statsmodels/statsmodels-skipper/statsmodels/stats/tests/]
[81]: optimize.root(func, .15, method="hybr")
[81]:
status: 1
success: True
qtf: array([-0.20178728])
nfev: 4
r: array([ inf])
fun: array([ 0.20178728])
x: array([ 0.15])
message: 'The solution converged.'
fjac: array([[-1.]])
> Josef
>
>
> >
> > from Skipper's run (
> https://github.com/statsmodels/statsmodels/issues/710 ):
> >
> > -----------
> > [~/statsmodels/statsmodels-skipper/statsmodels/stats/tests/]
> > [67]: optimize.fsolve(func, .14, full_output=True)
> > [67]:
> > (array([ 0.05]),
> > {'fjac': array([[-1.]]),
> > 'fvec': array([ 0.]),
> > 'nfev': 12,
> > 'qtf': array([ 0.]),
> > 'r': array([-3.408293])},
> > 1,
> > 'The solution converged.')
> >
> > [~/statsmodels/statsmodels-skipper/statsmodels/stats/tests/]
> > [68]: optimize.fsolve(func, .15, full_output=True)
> > [68]:
> > (array([ 0.15]),
> > {'fjac': array([[-1.]]),
> > 'fvec': array([ 0.20178728]),
> > 'nfev': 4,
> > 'qtf': array([-0.20178728]),
> > 'r': array([ inf])},
> > 1,
> > 'The solution converged.')
> >
> > It looks like the QR factorization is failing (r == inf) and then it's
> > reporting convergence still.
> > ---------
> >
> > That it stops with a "solution converged" also doesn't trigger my
> > backup root-finding.
> >
> >
> > I can work around this since I have no idea how to debug this.
> > However, there might be something fishy with fsolve.
> >
> > I never had problems with fsolve before, and scipy.stats.distributions
> > is a heavy user of it.
> >
> > Josef
> >
> >>
> >>>
> >>>>
> >>>> Josef
> >>>>
> >>>>>
> >>>>> Josef
> _______________________________________________
> SciPy-Dev mailing list
> SciPy-Dev at scipy.org
> http://mail.scipy.org/mailman/listinfo/scipy-dev
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mail.python.org/pipermail/scipy-dev/attachments/20130318/674f08fc/attachment.html>
More information about the SciPy-Dev
mailing list