[SciPy-Dev] SciPy-Dev Digest, Vol 113, Issue 21

The Helmbolds helmrp@yahoo....
Mon Mar 18 11:53:34 CDT 2013


Regarding fsolve, 

*Cautionary Note*: According to [2], the `hybrd` Fortran programs
computes an initial value of the "approximate Jacobian", but update
that value only when "the rank-1 method is not producing satisfactory
progress", where the reference is to Broyden's rank-1 method [4]_
(used in an `fsolve` inner loop). Examining the number of Jacobian 
evaluations used (see `njev`) will give some idea as to how often
`fsolve` updated the Jacobian. Because the program's QR-related
outputs (`fjac`, `r`, and `qtf`) are based on the program's
"approximate Jacobian", they should not be used in subsequent
analyses unless their validity is confirmed by independent
computations. (Example 2 illustrates this point.)

References
----------
.. [1] Powell, M. J. D., "A Hybrid Method for Nonlinear Equations,"
   in `Numerical Methods for Nonlinear Algebraic Equations`, P. 
   Rabinowitz, ed., Gordon and Breach, 1970.
.. [2] User Guide for MINPACK-1, Mor'e, Jorge J., et al, 
   Argonne National Laboratory, Report ANL-80-74, August 1980.
   (Out of print. For contents, see
   http://www.mcs.anl.gov/~more/ANL8074a.pdf, and 
   http://www.mcs.anl.gov/~more/ANL8074b.pdf.)
.. [3] A copy of the Fortran code is in SciPy module 
   "optimize/MINPACK".
.. [4] Broyden, C. G., "A Class of Methods for Solving Nonlinear
   Simultaneous Equations", Mathematics of Computation 
   (American Mathematical Society), 19 (92), October 1965, 
   pp 577–593. doi:10.2307/2003941. JSTOR 2003941.   

Examples
--------
`Example 1`: Supppose we wish to solve the following system of three
simultaneous nonlinear equations in the three variables (u, v, w):

* 2 * a*u + b*v + d - w*v = 0
* b*u + 2 * c*v + e - w*u = 0
* -u*v + g  = 0

For this system, the Jacobian matrix with derivatives 
across columns is:

    +-------+-------+----+
    |  2*a  | b - w | -v |
    +-------+-------+----+
    | b - w |  2*c  | -u |
    +-------+-------+----+
    |   -v  |   -u  | 0  |
    +-------+-------+----+

It is symmetrical, and could be considered as "banded", 
though we make no use of that in these examples. 
We may proceed as follows.

>>> params = (2, 3, 7, 8, 9, 10, 2)
>>> def Myfunc(z, *params):
>>> ... (u, v, w) = z
>>> ... (a, b, c, d, e, f, g) = params
>>> ... # Make a list of the LHS values.
>>> ... temp = [2*a*u + b*v + d - w*v]
>>> ... temp.append(b*u + 2*c*v + e - w*u)
>>> ... temp.append(-u*v + g)
>>> ... return temp
>>>
>>> def Myfprime(z, *params):
>>> ... (u, v, w) = z
>>> ... (a, b, c, d, e, f, g) = params
>>> ... J = np.ndarray([2*a, b - w, -v], 
            [b-w, 2*c, -u], 
            [-v, -u, 0])
>>> ... J.reshape((3,3))
>>> ... return J
>>>
>>> z0 = np.ndarray([1, 1, 5])         # Initial guess.
>>> # Accept the default values for all parameters.
>>> from scipy import optimize
>>> Example_1 = optimize.fsolve(Myfunc,
            z0,
            args = params)
>>> print "Example_1 = ", Example_1

Prints the solution point value::

[  1.79838825    1.11210691    16.66195357]


`Example 2`: To solve this problem using the new 
(as of Version 0.11.0) `root` method we need to (a) convert
the options to an "option dictionary" using keys acceptable 
to the `root` method for `fsolve` (as given in the
following code snippet),
(b) call `root` with `method = 'hybr'` (not `'hybrd'`),
and (c) take account of the fact that the returned value 
will be a `Result` dictionary as defined in SciPy's 
"optimize/optimize.py" module.
We may proceed as follows (where we note some changes 
in terminology from those used in `fsolve`):

>>> Myopts = {
          'col_deriv' : 0,  # Default, must be `0`, not `None`.
            'xtol'    : 1.49012e-8,   # Default.
            'maxfev'  : 0,    # Default.
            'band'    : None, # Default.
            'eps'     : 0.0,  # Default, formerly called `epsfcn`.
            'factor'  : 100,  # Default.
            'diag'    : None  # Default.
              }
>>>              
>>> from scipy import optimize
>>> Example_2 = optimize.root(Myfunc,
                z0,
                args = params,
                method = 'hybr',  # Not 'hybrd'.
                jac = Myprime,    # Formerly called 'fprime'.
                options = Myopts)
>>>
>>> print "Example_2 = ", Example_2
        
This prints the following (in a slightly different format)::

Example_2 =
    status: 1
    success: True
    qtf: array([ 1.32699176e-07, -8.49465473e-08, -3.09277588e-08])
    nfev: 12
    r: array([ -7.77806716, 30.02199802, -0.819055,
              -10.74878184, 2.00090268,   1.02706198])
    fun: array([ -2.30691910e-10, 4.58971172e-09, 5.55402391e-11])
    x: array([  1.79838825,   1.11210691,  16.66195357])
    message: 'The solution converged.'
    fjac: array([[-0.64093238,  0.75748326,  0.1241966 ],
              [-0.62403598, -0.60841098,  0.4903215 ],
              [-0.44697291, -0.23675978, -0.8626471 ]])
    njev: 1

As the `status` code, `success` value, and `message` indicate, the 
run is considered a success. The value of `fun` confirms this by 
showing that the equations are satisfied to at least 8 significant
figures. The solution point `x` is the same as in Example_1. 
The value of `nfev` shows that this run used 12 function 
evaluations.

Note that `njev`, the number of Jacobian evaluations, is 1. So the 
initial "approximate Jacobian" was used throughout. Thus the
reported values of `fjac`, `r`, and `qtf` may not be those at the
solution point. Therefore subsequent analyses should not use them
as the values at the solution point until their validity is
confirmed by independent computations.
 
Bob and Paula Helmbold
"Many of life's failures are people who did not realize how close they were to success when they gave up." (Thomas Edison)



>________________________________
> From: "scipy-dev-request@scipy.org" <scipy-dev-request@scipy.org>
>To: scipy-dev@scipy.org 
>Sent: Monday, March 18, 2013 8:36 AM
>Subject: SciPy-Dev Digest, Vol 113, Issue 21
> 
>Send SciPy-Dev mailing list submissions to
>    scipy-dev@scipy.org
>
>To subscribe or unsubscribe via the World Wide Web, visit
>    http://mail.scipy.org/mailman/listinfo/scipy-dev
>or, via email, send a message with subject or body 'help' to
>    scipy-dev-request@scipy.org
>
>You can reach the person managing the list at
>    scipy-dev-owner@scipy.org
>
>When replying, please edit your Subject line so it is more specific
>than "Re: Contents of SciPy-Dev digest..."
>
>
>Today's Topics:
>
>   1. Re: Lomb-Scargle Periodogram: Press & Rybicki Algorithm
>      (Christian Geier)
>   2. Re: optimize.fsolve endless loop with nan (josef.pktd@gmail.com)
>   3. Re: optimize.fsolve endless loop with nan (josef.pktd@gmail.com)
>   4. Re: optimize.fsolve endless loop with nan (Skipper Seabold)
>
>
>----------------------------------------------------------------------
>
>Message: 1
>Date: Mon, 18 Mar 2013 15:19:22 +0000
>From: Christian Geier <geier@lostpackets.de>
>Subject: Re: [SciPy-Dev] Lomb-Scargle Periodogram: Press & Rybicki
>    Algorithm
>To: SciPy Developers List <scipy-dev@scipy.org>
>Message-ID: <20130318151922.GA20372@brutus.lostpackets.de>
>Content-Type: text/plain; charset=us-ascii
>
>On Fri, Mar 08, 2013 at 05:18:56PM -0800, Jacob Vanderplas wrote:
>> Hi,
>> We have a cython version of this Lomb-Scargle algorithm in astroML [1], as
>> well as a generalized version that does not depend on the sample mean being
>> a good approximation of the true mean.  We've not yet implemented the FFT
>> trick shown in Press & Rybicki, but it's fairly fast as-is for problems of
>> reasonable size.
>> 
>> For some examples of it in use on astronomical data, see [2-3] below (the
>> examples are figures from our upcoming astro/statistics textbook). This
>> code is BSD-licensed, so if it seems generally useful enough to include in
>> scipy, it would be no problem to port it.
>> 
>> Also, if you have implemented an FFT-based version, it would get a fair bit
>> of use in the astronomy community if you were willing to contribute it to
>> astroML.
>> 
>> Thanks,
>>    Jake
>> 
>> [1]
>> https://github.com/astroML/astroML/blob/master/astroML_addons/periodogram.pyx
>> [2] http://astroml.github.com/book_figures/chapter10/fig_LS_example.html
>> [3]
>> http://astroml.github.com/book_figures/chapter10/fig_LS_sg_comparison.html
>> 
>> 
>> On Fri, Mar 8, 2013 at 5:26 PM, Christian Geier <geier@lostpackets.de>wrote:
>> 
>> > Hello everyone!
>> >
>> > Would you in general be considering to include the Lombscargle
>> > Periodogram by Press & Rybicki [1] into scipy in addition to the already
>> > present one? I find the included algorithm by Townend rather slow and
>> > had recently some "interesting" results returned by it.
>> >
>> > I've recently translated the original FORTRAN code (which is actually
>> > the description of the algorithm [1]) to (pure) python [2] and would like
>> > to know what the legal situation is in this case: can I release this
>> > translated code under the BSD license?
>> >
>> > In this case I would translate the code further to cython and supply
>> > tests and more documentation.
>> >
>> > Greetings
>> >
>> > Christian Geier
>> >
>> > [1] http://adsabs.harvard.edu/full/1989ApJ...338..277P
>> > [2]
>> > https://github.com/geier/scipy/commit/710bf4ca514d223df39891eb20401aba2edd8cdb
>
>Hi,
>thanks for your answer. When time allows for it (hopefully next week)
>I will improve the code some more and then contact astroML.
>
>Regards
>Christian
>
>
>------------------------------
>
>Message: 2
>Date: Mon, 18 Mar 2013 11:14:41 -0400
>From: josef.pktd@gmail.com
>Subject: Re: [SciPy-Dev] optimize.fsolve endless loop with nan
>To: SciPy Developers List <scipy-dev@scipy.org>
>Message-ID:
>    <CAMMTP+Bt9cNpGC0cB1ovF=XmEZV3kPHf0aoffbK7pkDxEt6EYA@mail.gmail.com>
>Content-Type: text/plain; charset=ISO-8859-1
>
>On Wed, Mar 13, 2013 at 1:24 PM,  <josef.pktd@gmail.com> wrote:
>> On Wed, Mar 13, 2013 at 12:59 PM,  <josef.pktd@gmail.com> wrote:
>>> On Wed, Mar 13, 2013 at 12:17 PM,  <josef.pktd@gmail.com> wrote:
>>>> On Wed, Mar 13, 2013 at 10:51 AM,  <josef.pktd@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)
>
>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
>
>
>------------------------------
>
>Message: 3
>Date: Mon, 18 Mar 2013 11:29:17 -0400
>From: josef.pktd@gmail.com
>Subject: Re: [SciPy-Dev] optimize.fsolve endless loop with nan
>To: SciPy Developers List <scipy-dev@scipy.org>
>Message-ID:
>    <CAMMTP+CWKGq4J+WbYBJouMrq+PynsLM2wZDc_Jy+8a=F04Ck9w@mail.gmail.com>
>Content-Type: text/plain; charset=ISO-8859-1
>
>On Mon, Mar 18, 2013 at 11:14 AM,  <josef.pktd@gmail.com> wrote:
>> On Wed, Mar 13, 2013 at 1:24 PM,  <josef.pktd@gmail.com> wrote:
>>> On Wed, Mar 13, 2013 at 12:59 PM,  <josef.pktd@gmail.com> wrote:
>>>> On Wed, Mar 13, 2013 at 12:17 PM,  <josef.pktd@gmail.com> wrote:
>>>>> On Wed, Mar 13, 2013 at 10:51 AM,  <josef.pktd@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
>
>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
>
>
>------------------------------
>
>Message: 4
>Date: Mon, 18 Mar 2013 11:40:49 -0400
>From: Skipper Seabold <jsseabold@gmail.com>
>Subject: Re: [SciPy-Dev] optimize.fsolve endless loop with nan
>To: SciPy Developers List <scipy-dev@scipy.org>
>Message-ID:
>    <CAKF=DjtzcPFUFJZvOONQsnSsy4mzQFWhTxyMP7Sw8nJ5cjmRbw@mail.gmail.com>
>Content-Type: text/plain; charset="iso-8859-1"
>
>On Mon, Mar 18, 2013 at 11:29 AM, <josef.pktd@gmail.com> wrote:
>
>> On Mon, Mar 18, 2013 at 11:14 AM,  <josef.pktd@gmail.com> wrote:
>> > On Wed, Mar 13, 2013 at 1:24 PM,  <josef.pktd@gmail.com> wrote:
>> >> On Wed, Mar 13, 2013 at 12:59 PM,  <josef.pktd@gmail.com> wrote:
>> >>> On Wed, Mar 13, 2013 at 12:17 PM,  <josef.pktd@gmail.com> wrote:
>> >>>> On Wed, Mar 13, 2013 at 10:51 AM,  <josef.pktd@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@scipy.org
>> http://mail.scipy.org/mailman/listinfo/scipy-dev
>>
>-------------- next part --------------
>An HTML attachment was scrubbed...
>URL: http://mail.scipy.org/pipermail/scipy-dev/attachments/20130318/674f08fc/attachment.html 
>
>------------------------------
>
>_______________________________________________
>SciPy-Dev mailing list
>SciPy-Dev@scipy.org
>http://mail.scipy.org/mailman/listinfo/scipy-dev
>
>
>End of SciPy-Dev Digest, Vol 113, Issue 21
>******************************************
>
>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://mail.scipy.org/pipermail/scipy-dev/attachments/20130318/8aee912d/attachment-0001.html 


More information about the SciPy-Dev mailing list