Re: [SciPy-dev] scipy.optimize.nonlin.broyden2
![](https://secure.gravatar.com/avatar/f65ea1dfe7ef7b9f889ed6877fcc69b8.jpg?s=120&d=mm&r=g)
On Wed, Mar 19, 2008 at 9:30 PM, argriffi <argriffi@ncsu.edu> wrote:
Hi Ondrej,
Running the attached script should demonstrate the problem I mentioned. I removed dependencies on other modules I had written, but let me know if you want it reduced even further.
Thanks, that is enough. I can reproduce the problem now: ------------------ test_solver_20 (__main__.TestCodonInversion) ... ok test_solver_40 (__main__.TestCodonInversion) ... ERROR ====================================================================== ERROR: test_solver_40 (__main__.TestCodonInversion) ---------------------------------------------------------------------- Traceback (most recent call last): File "CodonInversion.py", line 181, in test_solver_40 self._test_solver_helper(40) File "CodonInversion.py", line 173, in _test_solver_helper raise RuntimeError('failed on iteration %d of %d: %s' % (len(objective_function.guesses), iterations, str(e))) RuntimeError: failed on iteration 22 of 40: detected possibly invalid input to the objective function: (nan, nan, nan) ---------------------------------------------------------------------- Ran 2 tests in 0.029s FAILED (errors=1) ------------------------- And the reason is simple. If you look at this code: for n in range(iter): deltaxm=-Gm*Fxm xm=xm+deltaxm Fxm1=myF(F,xm) deltaFxm=Fxm1-Fxm Fxm=Fxm1 Gm=Gm+(deltaxm-Gm*deltaFxm)*deltaFxm.T/norm(deltaFxm)**2 And if you debug it, as I just did (for example using winpdb), you'll notice, that deltaFxm gets zero, after a lot of iterations, because you converged to a solution already. Then norm(deltaFxm) is 0 and after dividing by 0, you get a NaN in Gm, thus you get NaN in deltaxm and xm and you get NaN as an input into your function and here we are. The solution is to check that deltaFxm is near zero and stop iterating. Do you think you could please provide such a patch? :) Also with a test case. Thanks, Ondrej P.S. I CCed scipy-dev as I think you forgot to forward your email to it. I hope it's ok.
![](https://secure.gravatar.com/avatar/fc2baa22f0b0e9fc8dc93733cb94103c.jpg?s=120&d=mm&r=g)
Hi Ondrej, As requested I have attached a modification and a test case. The modification checks a value for equality to zero, but maybe it should more appropriately handle values near zero. Also, other functions in nonlin.py may benefit from similar modifications. Alex Ondrej Certik wrote:
On Wed, Mar 19, 2008 at 9:30 PM, argriffi <argriffi@ncsu.edu> wrote:
Hi Ondrej,
Running the attached script should demonstrate the problem I mentioned. I removed dependencies on other modules I had written, but let me know if you want it reduced even further.
Thanks, that is enough. I can reproduce the problem now:
------------------
test_solver_20 (__main__.TestCodonInversion) ... ok test_solver_40 (__main__.TestCodonInversion) ... ERROR
====================================================================== ERROR: test_solver_40 (__main__.TestCodonInversion) ---------------------------------------------------------------------- Traceback (most recent call last): File "CodonInversion.py", line 181, in test_solver_40 self._test_solver_helper(40) File "CodonInversion.py", line 173, in _test_solver_helper raise RuntimeError('failed on iteration %d of %d: %s' % (len(objective_function.guesses), iterations, str(e))) RuntimeError: failed on iteration 22 of 40: detected possibly invalid input to the objective function: (nan, nan, nan)
---------------------------------------------------------------------- Ran 2 tests in 0.029s
FAILED (errors=1)
-------------------------
And the reason is simple. If you look at this code:
for n in range(iter): deltaxm=-Gm*Fxm xm=xm+deltaxm Fxm1=myF(F,xm) deltaFxm=Fxm1-Fxm Fxm=Fxm1 Gm=Gm+(deltaxm-Gm*deltaFxm)*deltaFxm.T/norm(deltaFxm)**2
And if you debug it, as I just did (for example using winpdb), you'll notice, that deltaFxm gets zero, after a lot of iterations, because you converged to a solution already. Then norm(deltaFxm) is 0 and after dividing by 0, you get a NaN in Gm, thus you get NaN in deltaxm and xm and you get NaN as an input into your function and here we are.
The solution is to check that deltaFxm is near zero and stop iterating. Do you think you could please provide such a patch? :)
Also with a test case.
Thanks, Ondrej
P.S. I CCed scipy-dev as I think you forgot to forward your email to it. I hope it's ok.
import unittest import scipy.optimize import numpy import math def easy_function(X): return tuple(x*x - 1 for x in X) class MyTest(unittest.TestCase): def helper(self, solve): guess = (0, 0, 0) iterations = 20 result = solve(easy_function, guess, iterations) for x in result: self.failUnless(x in (-1.0, 1.0), x) def test_current(self): "The current version of broyden2 fails." self.helper(scipy.optimize.nonlin.broyden2) def test_modified(self): "A modified version of broyden2 does not fail so hard." def norm(v): return math.sqrt(numpy.sum((numpy.array(v)**2).flat)) def myF(F,xm): return numpy.matrix(F(tuple(xm.flat))).T def broyden2(F, xin, iter=10, alpha=0.4, verbose = False): xm=numpy.matrix(xin).T Fxm=myF(F,xm) Gm=-alpha*numpy.matrix(numpy.identity(len(xin))) for n in range(iter): deltaxm=-Gm*Fxm xm=xm+deltaxm Fxm1=myF(F,xm) # begin modification deltaFxm=Fxm1-Fxm denominator = norm(deltaFxm)**2 if not denominator: break Fxm=Fxm1 Gm=Gm+(deltaxm-Gm*deltaFxm)*deltaFxm.T/denominator # end modification if verbose: print "%d: |F(x)|=%.3f"%(n+1, norm(Fxm)) return xm.flat self.helper(broyden2) unittest.main()
![](https://secure.gravatar.com/avatar/f65ea1dfe7ef7b9f889ed6877fcc69b8.jpg?s=120&d=mm&r=g)
On Thu, Mar 20, 2008 at 4:11 AM, argriffi <argriffi@ncsu.edu> wrote:
Hi Ondrej,
As requested I have attached a modification and a test case. The modification checks a value for equality to zero, but maybe it should more appropriately handle values near zero. Also, other functions in nonlin.py may benefit from similar modifications.
Excellent, thanks a lot for the work. Here are my comments: denominator = norm(deltaFxm)**2 if not denominator: break use abs(denominator) < eps and define eps as a default parameter of the function, for example eps=1e-10. Also could you please apply this into the scipy svn, and send use the output of "svn di"? So that someone with a svn access can just take it and apply it. As to the tests, please add them to scipy/optimize/tests/test_nonlin.py. Simply do all of this in your local svn copy and send use the result of "svn di". Thanks very much, Ondrej
![](https://secure.gravatar.com/avatar/f65ea1dfe7ef7b9f889ed6877fcc69b8.jpg?s=120&d=mm&r=g)
On Thu, Mar 20, 2008 at 1:45 PM, Ondrej Certik <ondrej@certik.cz> wrote:
On Thu, Mar 20, 2008 at 4:11 AM, argriffi <argriffi@ncsu.edu> wrote:
Hi Ondrej,
As requested I have attached a modification and a test case. The modification checks a value for equality to zero, but maybe it should more appropriately handle values near zero. Also, other functions in nonlin.py may benefit from similar modifications.
Excellent, thanks a lot for the work. Here are my comments:
denominator = norm(deltaFxm)**2 if not denominator: break
use abs(denominator) < eps
and define eps as a default parameter of the function, for example eps=1e-10.
Also could you please apply this into the scipy svn, and send use the
use -> us
output of "svn di"? So that someone with a svn access can just take it and apply it.
As to the tests, please add them to scipy/optimize/tests/test_nonlin.py. Simply do all of this in your local svn copy and send use the result of "svn di".
use -> us. O.
![](https://secure.gravatar.com/avatar/fc2baa22f0b0e9fc8dc93733cb94103c.jpg?s=120&d=mm&r=g)
I've been using a packaged scipy release, and I just checked out the svn scipy trunk but haven't used or modified it yet. What is the best way to manage the two installations so that by default my scripts use the official release, but so that I can alternatively use the svn version by specifying it somehow? I also have a question about the proposed eps argument for broyden2. Could this be named more descriptively? For example, I have seen small values named rel_tol and abs_tol for relative and absolute tolerance in other contexts. If eps has a meaningful interpretation in broyden2 then maybe it should be named more appropriately. Alex Ondrej Certik wrote:
On Thu, Mar 20, 2008 at 4:11 AM, argriffi <argriffi@ncsu.edu> wrote:
Hi Ondrej,
As requested I have attached a modification and a test case. The modification checks a value for equality to zero, but maybe it should more appropriately handle values near zero. Also, other functions in nonlin.py may benefit from similar modifications.
Excellent, thanks a lot for the work. Here are my comments:
denominator = norm(deltaFxm)**2 if not denominator: break
use abs(denominator) < eps
and define eps as a default parameter of the function, for example eps=1e-10.
Also could you please apply this into the scipy svn, and send use the output of "svn di"? So that someone with a svn access can just take it and apply it.
As to the tests, please add them to scipy/optimize/tests/test_nonlin.py. Simply do all of this in your local svn copy and send use the result of "svn di".
Thanks very much, Ondrej _______________________________________________ Scipy-dev mailing list Scipy-dev@scipy.org http://projects.scipy.org/mailman/listinfo/scipy-dev
![](https://secure.gravatar.com/avatar/f65ea1dfe7ef7b9f889ed6877fcc69b8.jpg?s=120&d=mm&r=g)
On Thu, Mar 20, 2008 at 4:12 PM, argriffi <argriffi@ncsu.edu> wrote:
I've been using a packaged scipy release, and I just checked out the svn scipy trunk but haven't used or modified it yet. What is the best way to manage the two installations so that by default my scripts use the official release, but so that I can alternatively use the svn version by specifying it somehow?
Maybe there are other ways how to deal with this, but I do: ./setup install --home=/some/path everytime I change something - it takes just a while to do it if you only modify python files. And then you add /some/path into your PYTHONPATH - just temporarily.
I also have a question about the proposed eps argument for broyden2. Could this be named more descriptively? For example, I have seen small values named rel_tol and abs_tol for relative and absolute tolerance in other contexts. If eps has a meaningful interpretation in broyden2 then maybe it should be named more appropriately.
Absolutely, name it more appropriate. That was just a general suggestion from me how to do things. Ondrej
participants (2)
-
argriffi
-
Ondrej Certik