BUG: scipy.optimize.minimize with method="newton-cg" is stuck on a local maximum

[legandro-github]

Preface

I either have a usage problem or there is something wrong with the source code. Because I'm unsure which is the case I both opened a stackoverflow question and this Issue. I will just copy the text I've already written on stackoverflow to here.

tl;dr

I want to find the local minimum for a function that depends on 2 variables. For that my plan was to use the scipy.optimize.minimize function with the "newton-cg" method because I can calculate the jacobian and hessian analytically. However, when my starting guesses are on a local maximum, the function terminates successfully in the first iteration step on top of the local maximum, even though the hessian is negative.

Reproducing code

The reproducing code and the wrongly successfull termination message can be seen below.

The same outcome occurs if I use hess=None, hess='cs', hess='2-point' or hess='3-point' instead of my custom hess function. Also if I use other methods like 'dogleg', 'trust-ncg', 'trust-krylov', 'trust-exact' or 'trust-constr' I basically get the same outcome except nhev = 1 but the outcome is still wrong with x = [0,0].

Either I am doing something terribly wrong here (definitely very likely) or there is a major problem with the minimize function, specifically the "newton-cg" method (very unlikely?).

Regarding the latter case I also checked the source code to see if something's wrong there and stumpled upon something kinda weird(?). However, I don't completely understand the whole code, so I am a bit unsure if my worries are legitimate:

Taking a look at the source code

When the minimize function is called with method="newton-cg" it jumps into the _minimize_newtoncg function (see source code here). I want to go into detail what I believe happens here (so you can tell where I am potentially wrong):

On line 2168 A = sf.hess(xk) the hessian is first calculated in dependence of xk which is at first the start guess x0. For my test case the hessian is of course

A = [[f_xx, f_xy], [f_xy, f_yy]]

with f_ij being the derivatives of f after i and j. In my case f_xy = f_yx is also true.

Next on line 2183 Ap = A.dot(psupi) the product of the hessian A and psupi is calculated. psupi is basically equal to b which is the negative gradient of f at xk. So Ap = A.dot(psupi) results in

Ap = [f_xx f_x + f_xy f_y, f_xy f_x + f_yy f_y].

Now to the (possible) problem

Next, the curvature curv is calculated on line 2186 by np.dot(psupi, Ap). As explained above, psupis the negative gradient of f so this results in

curv = f_xx f_x^2 + 2 f_xy f_x f_y + fyy_f_y^2.

However, all of these derivatives are at xk which is equal to the starting parameters x0 at first. If the starting parameters are exactly at a local maximum, the derivatives f_x and f_y are equal to 0. Because of this, curv = 0. This results in a for loop break on the next line, thus, skipping to update xsupi, psupi and all the other parameters. Therefore, pk becomes [0,0] and _line_search_wolfe12 is called with basically all start parameters. This is where my understanding of the source code stops, however, I feel like things already went wrong after curv = 0 and breaking the for loop.

My concluding question is: What am I doing wrong leading to the minimize function being stuck on a local maximum?

Small update: I now understand that this case of starting at a local maximum and finding a local minimum from there might not have been implemented. So I am unsure if that would be doable to add in the future. If it were, I'd greatly appreciate it.

Reproducing Code Example

import numpy as np
import scipy.optimize as o

def get_f_df(var):
    x = var[0]
    y = var[1]

    f = np.cos(x) + np.cos(y)
    df_dx = -np.sin(x)
    df_dy = -np.sin(y)

    return f, (df_dx, df_dy)

def hess(var):
    x = var[0]
    y = var[1]

    f_hess = np.zeros((2,2))
    f_hess[0,0] = -np.cos(x)
    f_hess[1,1] = -np.cos(y)
    return f_hess

min = o.minimize(get_f_df, (0, 0), jac=True, hess=hess, method="newton-cg")
print(min)

Error message (not an error but the OptimizedResult)

message: Optimization terminated successfully.
 success: True
  status: 0
     fun: 2.0
       x: [ 0.000e+00  0.000e+00]
     nit: 1
     jac: [-0.000e+00 -0.000e+00]
    nfev: 1
    njev: 1
    nhev: 1

SciPy/NumPy/Python version and system information

1.12.0 1.26.4 sys.version_info(major=3, minor=11, micro=0, releaselevel='final', serial=0)
Build Dependencies:
  blas:
    detection method: pkgconfig
    found: true
    include directory: /c/opt/64/include
    lib directory: /c/opt/64/lib
    name: openblas
    openblas configuration: USE_64BITINT= DYNAMIC_ARCH=1 DYNAMIC_OLDER= NO_CBLAS=
      NO_LAPACK= NO_LAPACKE= NO_AFFINITY=1 USE_OPENMP= SKYLAKEX MAX_THREADS=2    
    pc file directory: c:/opt/64/lib/pkgconfig
    version: 0.3.21.dev
  lapack:
    detection method: pkgconfig
    found: true
    include directory: /c/opt/64/include
    lib directory: /c/opt/64/lib
    name: openblas
    openblas configuration: USE_64BITINT= DYNAMIC_ARCH=1 DYNAMIC_OLDER= NO_CBLAS=
      NO_LAPACK= NO_LAPACKE= NO_AFFINITY=1 USE_OPENMP= SKYLAKEX MAX_THREADS=2    
    pc file directory: c:/opt/64/lib/pkgconfig
    version: 0.3.21.dev
  pybind11:
    detection method: config-tool
    include directory: unknown
    name: pybind11
    version: 2.11.1
Compilers:
  c:
    commands: cc
    linker: ld.bfd
    name: gcc
    version: 10.3.0
  c++:
    commands: c++
    linker: ld.bfd
    name: gcc
    version: 10.3.0
  cython:
    commands: cython
    linker: cython
    name: cython
    version: 3.0.8
  fortran:
    commands: gfortran
    linker: ld.bfd
    name: gcc
    version: 10.3.0
  pythran:
    include directory: C:\Users\runneradmin\AppData\Local\Temp\pip-build-env-d02me32y\overlay\Lib\site-packages/pythran
    version: 0.15.0
Machine Information:
  build:
    cpu: x86_64
    endian: little
    family: x86_64
    system: windows
  cross-compiled: false
  host:
    cpu: x86_64
    endian: little
    family: x86_64
    system: windows
Python Information:
  path: C:\Users\runneradmin\AppData\Local\Temp\cibw-run-aug4qibh\cp311-win_amd64\build\venv\Scripts\python.exe
  version: '3.11'