BUG: fsolve + pchip + max, can't find root

[nealockwood]

Describe your issue.

I originally posted this on stack thinking the error might be mine but now I think it's some sort of bug.

I am trying to find the roots of multiple equations over a grid (want the equations to be 0 evaluated at the points in the grid) and am stacking all the function values as a vector. I was able to do this in Matlab,

After playing around with every aspect of the code, I discovered the max functions are giving it trouble. The max terms are squared to make them a smooth function, but when they are squared, fsolve just returns the initial guess. Ironically, when the squared term is removed, fsolve does return something besides the initial guess, but it's not the root. When the max terms (and mu terms) are completely removed, fsolve successfully finds a root.

There is also a pchip being used, which may add some complications, but the same thing is happening when I use a spline

Reproducing Code Example

import numpy as np
from scipy.interpolate import pchip, Akima1DInterpolator
from scipy.interpolate import InterpolatedUnivariateSpline as spline

def ap(x, alpha, beta, delta, r, w, kgrid, zgrid, piz):
    m = len(kgrid)
    pp1 = pchip(kgrid, x[:m])
    pp2 = pchip(kgrid,x[m:2*m])
    #pp1 = spline(kgrid, x[:m],k=3)
    #pp2 = spline(kgrid,x[m:2*m],k=3)
    res = np.zeros(len(x))
    for i in range(m):
        kp1 = x[i]
        c = (1 + r - delta) * kgrid[i] + w * zgrid[0] - kp1
        kpp1 = pp1(kp1)
        cp1 = (1 + r - delta) * kp1 + w * zgrid[0] - kpp1
        kpp2 = pp2(kp1)
        cp2 = (1 + r - delta) * kp1 + w * zgrid[1] - kpp2
        mu1 = x[i + 2*m]
        res[i] = c**(-1) - max(mu1, 0)**2 - beta * (1 + r - delta) * (piz[0, 0] * cp1**(-1) + piz[0, 1] * cp2**(-1))
        kp2 = x[i + m]
        c = (1 + r - delta) * kgrid[i] + w * zgrid[1] - kp2
        kpp1 = pp1(kp2)
        cp1 = (1 + r - delta) * kp2 + w * zgrid[0] - kpp1
        kpp2 = pp2(kp2)
        cp2 = (1 + r - delta) * kp2 + w * zgrid[1] - kpp2
        mu2 = x[i + 3*m]
        res[i + m] = c**(-1) - max(mu2, 0)**2 - beta * (1 + r - delta) * (piz[1, 0] * cp1**(-1) + piz[1, 1] * cp2**(-1))
        res[i + 2*m] = max(-mu1, 0)**2 - kp1
        res[i + 3*m] = max(-mu2, 0)**2 - kp2
    return res

Is the function and the evaluation is

alpha = 0.36
beta = 0.99
delta = 0.025
zgrid = np.array([1.25, 0.75])
piz = np.array([[0.9, 0.1], [0.1, 0.9]])
m = 51
kgrid = np.concatenate((np.linspace(0, 2, 20), np.linspace(2.5, 25, 20), np.linspace(30, 200, 11)))
theta = 1.8  
r=.034
w = (1-alpha)*(r/alpha)**(alpha/(alpha-1))
x0 = np.zeros(4*m)
x0[:m] = kgrid
x0[m:2*m] = kgrid
x0[2*m:3*m] = -np.sqrt(kgrid)
x0[3*m:] = -np.sqrt(kgrid)

def to_solve(x):
    return ap(x, alpha, beta, delta, r, w, kgrid, zgrid, piz)

x, f = fsolve(to_solve, x0, full_output=True,xtol=1e-12)[:2]

Changing the tolerance didn't help. I also tried simplifying and writing like

data = alpha, beta, delta, r, w, kgrid, zgrid, piz

x = fsolve(ap,x0,args=data)

But it changed nothing. I looked a long time for something on this, searched github for files that used this combination of functions, but to no avail.

SciPy/NumPy/Python version and system information

1.10.1 1.24.3 sys.version_info(major=3, minor=11, micro=4, releaselevel='final', serial=0)
lapack_armpl_info:
  NOT AVAILABLE
lapack_mkl_info:
    libraries = ['mkl_rt']
    library_dirs = ['C:/ProgramData/anaconda3\\Library\\lib']
    define_macros = [('SCIPY_MKL_H', None), ('HAVE_CBLAS', None)]
    include_dirs = ['C:/ProgramData/anaconda3\\Library\\include']
lapack_opt_info:
    libraries = ['mkl_rt']
    library_dirs = ['C:/ProgramData/anaconda3\\Library\\lib']
    define_macros = [('SCIPY_MKL_H', None), ('HAVE_CBLAS', None)]
    include_dirs = ['C:/ProgramData/anaconda3\\Library\\include']
blas_armpl_info:
  NOT AVAILABLE
blas_mkl_info:
    libraries = ['mkl_rt']
    library_dirs = ['C:/ProgramData/anaconda3\\Library\\lib']
    define_macros = [('SCIPY_MKL_H', None), ('HAVE_CBLAS', None)]
    include_dirs = ['C:/ProgramData/anaconda3\\Library\\include']
blas_opt_info:
    libraries = ['mkl_rt']
    library_dirs = ['C:/ProgramData/anaconda3\\Library\\lib']
    define_macros = [('SCIPY_MKL_H', None), ('HAVE_CBLAS', None)]
    include_dirs = ['C:/ProgramData/anaconda3\\Library\\include']
Supported SIMD extensions in this NumPy install:
    baseline = SSE,SSE2,SSE3
    found = SSSE3,SSE41,POPCNT,SSE42,AVX,F16C,FMA3,AVX2
    not found = AVX512F,AVX512CD,AVX512_SKX,AVX512_CLX,AVX512_CNL,AVX512_ICL

[nealockwood]

Someone on stack found that if number of iterations was high enough, scipy.optimize.minimize() will return a non initial guess. But if you just run the default it returns the initial guess after about a minute (Matlab finds the root in about 3 seconds on the same machine)

[nealockwood]

@mdhaber is it possible for this to get tagged?

Perhaps I didn't meet some threshold for posting a minimal version of the code. But I don't know how to simplify it. Namely, this system (of equations) requires two pchip interpolation components, so if I wanted to remove one of the max terms, (comment out mu2, res[i+m], and res[i+3m]) I would have more unknowns than equations. Moreover, I know that the code I presented has a root, and even if I found a way to simplify, it may not

[mdhaber]

@nealockwood hmm I'm not seeing evidence of a bug here. It's not surprising that a gradient-based rootfinder is having trouble with a function that involves max, because max is not smooth. You could try a soft maximum function and/or try scipy.optimize.root with another method.

I am trying to find the roots of multiple equations over a grid (want the equations to be 0 evaluated at the points in the grid) and am stacking all the function values as a vector.

How many inputs and outputs does each equation have? That is, for each equation, how many values are varied to find the root, and how many outputs need to become zero?

If each of your equations is scalar-in, scalar-out, you have other options.

[nealockwood]

@mdhaber thank you so much for quickly responding!

To briefly highlight why I think it's a bug: fsolve in matlab is able to find the root pretty fast, the max terms are squared (making them smooth functions), removing the squared term actually seems to improve the performance (returning something other than the initial guess), I have unsuccessfully tried other max alternatives, , and one is able to get a non-trivial solution by treating it as a minimization problem (minimizing normed value) when given a lot of iterations.

Each function has a scalar output, but the goal is to solve the system over a fixed grid. So essentially, there are 4 equations I want to be 0 over a grid of m points. To do this, I stack the residual terms in a vector. So for example, with m=51, there are 204 residual terms, or you could think of it as each grid point corresponding to 4 of the 204 residual terms.

[mdhaber]

The documentation of Matlab's fsolve says it uses a trust-region dogleg algorithm "similar in nature" to the that of MINPACK. This suggests that it is not identical to the one in MINPACK, which scipy.optimize.fsolve wraps. It is not guaranteed for all algorithms to solve all problems from an arbitrary guess, so failure to solve this problem does not necessarily indicate a bug. Often, it is just a shortcoming of the algorithm.

Assuming that max is from math (not NumPy), I see that max(x, 0)**2 is continuous and has continuous first derivative. However, the function may not be sufficiently smooth. The second derivative is discontinuous, and the first derivative is 0 for half the real line, which I would expect to cause problems.

You might try different algorithms (exposed with scipy.optimize.root) or a minimization algorithm - maybe a global minimization algorithm, even. I can't guarantee that any will solve the problem, though, even if there is a solution.

There is an effort to replace the Fortran code of MINPACK (#18566), which we can hope will uncover and fix bugs in the implementation. In the meantime, I don't think there is much we can do to illuminate bugs with this example. Would you agree @ilayn?

[nealockwood]

Thanks for the detail!

[ilayn]

Our past and current role was/is to only facilitate/vendor these external Fortran77 libraries and hope that these are internally consistent and tested over time by actual experts.

As we started looking into these old libraries a bit more critically, it is becoming more apparent that there are a lot of edge cases that such historical code do not handle and do not report back. However we are also bounded by our own expertise in all different subjects hence it would be typical that we would not catch up the current state of art in integration or optimization and so on unfortunately.

But our hope is that as we convert these tools to a more readable format, more experts can look into the code and spot obvious points of improvement, low-hanging-fruit fixes etc. or get more encouraged to teach us how to do things better. Also we would be able to see what happens when we push in functions like max() and see if we can incorporate some tiny bumps here and there. Currently, we are still far from that point, but getting there. All help is welcome.

[nealockwood]

I will see if I can recruit some help!