ZeroDivisionError in broyden1 when initial guess is the right answer. #3227

Closed
rmjarvis opened this Issue Jan 20, 2014 · 4 comments

Projects

None yet

4 participants

@rmjarvis

I'm getting a ZeroDivisionError when I use scipy.optimize.broyden1 and the initial guess happens to be the right answer already -- not so uncommon a situation in real life, so I'm kind of surprised that this issue hasn't been reported already. Anyway, here is a code snippet to reproduce the problem:

import scipy.optimize
xsq = [ 1, 4, 9 ]  # Solve x**2 = xsq for x.
xinit = [ 1, 2, 3 ]  # Initial guess.  Also the right answer.
func = lambda x: x**2 - xsq  # The function to pass to broyden1
soln = scipy.optimize.broyden1(func, xinit, x_tol=1.e-6)
print 'soln = ',soln

I'm using scipy version 0.13.2 with Mac OS 10.6.8 (using the fink installation of scipy-py27), and I get the following output:

Traceback (most recent call last):
  File "test_broyden.py", line 8, in <module>
    soln = scipy.optimize.broyden1(func, xinit, x_tol=1.e-6)
  File "<string>", line 8, in broyden1
  File "/sw64/lib/python2.7/site-packages/scipy/optimize/nonlin.py", line 284, in nonlin_solve
    jacobian.setup(x.copy(), Fx, func)
  File "/sw64/lib/python2.7/site-packages/scipy/optimize/nonlin.py", line 965, in setup
    GenericBroyden.setup(self, x, F, func)
  File "/sw64/lib/python2.7/site-packages/scipy/optimize/nonlin.py", line 672, in setup
    self.alpha = 0.5*max(norm(x0), 1) / norm(f0)
ZeroDivisionError: float division by zero

I think the fix is trivial. Just check if norm(f0) == 0, and set alpha to 1 (or whatever) in that case.

The workaround until this is fixed is to make sure to set alpha to something reasonable. Passing alpha=1 to broyden1 works just fine in this case.

@rmjarvis rmjarvis referenced this issue in astropy/astropy Jan 20, 2014
Closed

Problems with wcs.all_world2pix function. #1977

@nicodelpiano

Do you think it's correct to set alpha to 1?
Mathematically speaking, if norm(f0) is 0, then we must set alpha to the maximum value...
I'm new here, and interested to start solving issues :).

Cheers!

@rmjarvis

I believe it shouldn't matter what you set alpha to when f0=0, since we're already at the right answer, so we just needs to get to the next step where it should pass the convergence criterion and return.

@argriffing

@nicodelpiano It looks like alpha is just some initial guess related to the jacobian, so if the initial guess of x0 is correct then we do not care about the jacobian.
https://github.com/scipy/scipy/blob/master/scipy/optimize/nonlin.py#L881
I don't have much familiarity with this code but it seems OK to just set alpha to 1 if alpha is initially None and if the default initialization of alpha tries to divide by zero because the x0 guess is correct. Maybe add a unit test too.

I'm new here, and interested to start solving issues :).

@argriffing

submitted a PR

@ev-br ev-br closed this in #3390 Feb 25, 2014
@pv pv added this to the 0.14.0 milestone Feb 25, 2014
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment