[SciPy-User] possible bug in scipy 0.10.1 newton_krylon solver

[Wu Degang]

Hi,

I'm using scipy on my macbook pro with mac os x lion through macports 
(py27-scipy). I'm not sure whether it is a platform-specific problem. 
Here is a program that should reproduces the bug:


import numpy as np
import math
from scipy.optimize import newton_krylov
from numpy import zeros_like, zeros
from math import log
import pickle

with open('bug_data','rb') as pkl_file:
     k_seq = pickle.load(pkl_file)

N = len(k_seq)

def func(ws):
     new_lhs = zeros_like(ws)

     for i in xrange(N):
         new_lhs[i] = k_seq[i] - (N-1) + np.sum(1/(1+np.exp(ws[i]+ws))) 
- 1/(1+np.exp(2*ws[i]))

     return new_lhs

guess = zeros(N, float)
sol = newton_krylov(func,guess,verbose=1)


On my MBP, it produces the following output:

0:  |F(x)| = 1308.14; step 1; tol 8.46824e-05
bug.py:17: RuntimeWarning: overflow encountered in exp
   new_lhs[i] = k_seq[i] - (N-1) + np.sum(1/(1+np.exp(ws[i]+ws))) - 
1/(1+np.exp(2*ws[i]))
1:  |F(x)| = 181529; step 1; tol 0.9999
2:  |F(x)| = 128432; step 1; tol 0.89982
3:  |F(x)| = 128000; step 1; tol 0.893968
4:  |F(x)| = 89085; step 1; tol 0.719261
/opt/local/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/site-packages/scipy/sparse/linalg/isolve/lgmres.py:264: 
RuntimeWarning: invalid value encountered in divide
   outer_v.append((dx/nx, None))
Traceback (most recent call last):
   File "bug.py", line 22, in <module>
     sol = newton_krylov(func,guess,verbose=1)
   File "<string>", line 8, in newton_krylov
   File 
"/opt/local/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/site-packages/scipy/optimize/nonlin.py", 
line 297, in nonlin_solve
     raise ValueError("Jacobian inversion yielded zero vector. "
ValueError: Jacobian inversion yielded zero vector. This indicates a bug 
in the Jacobian approximation.


To data file (a pickle file actually) can be found here (I don't know 
whether it is possible to upload file, perhaps attachment?): 
https://dl.dropbox.com/u/7383429/bug_data. It contains a numpy array.


I also have a five-year old laptop running xubuntu 12.04 with scipy 0.9. 
The same program and the same data gives a very different output:

0:  |F(x)| = 32140.1; step 1; tol 0.0511187
1:  |F(x)| = 11196.2; step 1; tol 0.109217
2:  |F(x)| = 4105.06; step 1; tol 0.120988
3:  |F(x)| = 1787.21; step 1; tol 0.170589
4:  |F(x)| = 574.511; step 1; tol 0.0930013
5:  |F(x)| = 209.088; step 1; tol 0.119207
6:  |F(x)| = 68.6437; step 1; tol 0.0970035
7:  |F(x)| = 20.731; step 1; tol 0.0820881
8:  |F(x)| = 3.55224; step 1; tol 0.0264245
9:  |F(x)| = 0.109991; step 1; tol 0.000862892
10:  |F(x)| = 7.99967e-05; step 1; tol 4.76067e-07
11:  |F(x)| = 4.1262e-08; step 1; tol 2.39441e-07

It runs successfully.


Thanks.

Wu Degang