Strange fsolve behavior when fvec array is reused

[tadeu]

fsolve seems to be messing with the returned fvec array, because when it is a "reused" array, fsolve have a strange behavior (after some function evaluations it passes weird x values, and diverges).

Reproducing code example:

import numpy as np
from scipy.optimize import fsolve

y = np.full(1, np.nan)

def f(x):
    y[:] = x * x - np.log(x + 1)
    print(f'x: {x}, y: {y}')
    return y

x = fsolve(f, x0 = 1.0)

print(f'\nfinal x: {x}')

Output:

x: [1.], y: [0.30685282]                                                                                               
x: [1.], y: [0.30685282]                                                                                               
x: [1.], y: [0.30685282]                                                                                               
x: [1.00000001], y: [0.30685284]                                                                                       
test_scipy.py:7: RuntimeWarning: invalid value encountered in log                                                      
  y[:] = x * x - np.log(x + 1)                                                                                         
x: [-99.], y: [nan]                                                                                                    
x: [nan], y: [nan]                                                                                                     
x: [1.00000001], y: [0.30685284]                                                                                       
x: [-99.], y: [nan]                                                                                                    
x: [nan], y: [nan]                                                                                                     
x: [nan], y: [nan]                                                                                                     
x: [nan], y: [nan]                                                                                                     
x: [nan], y: [nan]                                                                                                     
x: [nan], y: [nan]                                                                                                     
x: [nan], y: [nan]                                                                                                     
x: [nan], y: [nan]                                                                                                     
(...)\site-packages\scipy\optimize\minpack.py:162: RuntimeWarning: The iteration is not making good progress, as measured by the improvement from the last ten iterations.
  warnings.warn(msg, RuntimeWarning)                                                                                   
                                                                                                                       
final x: [1.]

It passed a -99. at the fifth function evaluation above. Compare with the results if we change y[:] = x * x - np.log(x + 1) for y = x * x - np.log(x + 1) (so, creating a new fvec array every time):

x: [1.], y: [0.30685282]            
x: [1.], y: [0.30685282]            
x: [1.], y: [0.30685282]            
x: [1.00000001], y: [0.30685284]    
x: [0.79543146], y: [0.04746584]    
x: [0.75799697], y: [0.01038433]    
x: [0.74751379], y: [0.00058278]    
x: [0.74689048], y: [8.05201074e-06]
x: [0.74688175], y: [6.42351716e-09]
x: [0.74688174], y: [7.09432513e-14]
x: [0.74688174], y: [1.11022302e-16]
                                    
final x: [0.74688174]

Scipy/Numpy/Python version information:

1.2.0 1.15.4 sys.version_info(major=3, minor=7, micro=2, releaselevel='final', serial=0)

[tadeu]

This is probably because of internal behavior from minpack's HYBRJ, e.g.:

         IF IFLAG = 1 CALCULATE THE FUNCTIONS AT X AND
         RETURN THIS VECTOR IN FVEC.  DO NOT ALTER FJAC.
         IF IFLAG = 2 CALCULATE THE JACOBIAN AT X AND
         RETURN THIS MATRIX IN FJAC.  DO NOT ALTER FVEC.

but if it's too hard to circumvent (detect that fvec still uses the same buffer and making a copy in that case, etc.), perhaps a note could be added to the documentation.

[volpatto]

@tadeu, I think you would like to know:

Just as a note, I tried to reproduce this issue locally. Instead using fsolve, I used root, which is stated as the up-to-date interface for root-finding methods in scipy.optimize. The same behavior is observed with method="hybr". However, when you choose method="krylov", the strange behavior is not observed. Below I provide the MWE:

import numpy as np
from scipy.optimize import root

y = np.full(1, np.nan)

def f(x):
    y[:] = x * x - np.log(x + 1)
    print(f'x: {x}, y: {y}')
    return y

x = root(f, x0 = 1.0, method='krylov')
# x = root(f, x0 = 1.0, method='hybr')  # uncomment to check hybr behavior

print(f"\noutput:\n{x}")