How to stop the iteration when the Jacobian reached to an arbitrary (small) value in Newton-CG method?

Question

How to put a stopping condition on jacobian (or gradient) for Newton-CG methode? I want the algorithme to stop when the jacobian reaches to 1e-2, is it possible to do with Newton-CG ??

input:

scipy.optimize.minimize(f, [5.0,1.0,2.0,5.0], args=Data, method='Newton-CG',jac=Jacf)

output:

     jac: array([7.64265411e-08, 1.74985718e-08, 4.12408407e-07, 5.02972841e-08])
 message: 'Optimization terminated successfully.'
    nfev: 12
    nhev: 0
     nit: 11
    njev: 68
  status: 0
 success: True
       x: array([0.22545395, 0.3480084 , 1.06811724, 1.64873479])

in BFGS method, which is symilar to Newton-CG, there is a gtol option, it allows to stop the iteration when the gradient reaches to some value. But in Newton-CG theres no that type of option.

Does anyone know how to stop the iteration when the jacobien reaches to 1e-2.

Here are some details to reproduce my code:

def convert_line2matrix(a):
    n = len(a)
    if (np.sqrt(n) % 1 == 0) :
        d = int(np.sqrt(n))
        Mat = np.zeros((d,d))
        for i in range(d):
            for j in range(d):
                Mat[i,j] = a[j+d*i] 
    else:
        raise  ValueError(f"{a} cant be converted into a (n x n) matrix. The array has {len(a)} elements, \n\t    thus impossible to build a square matrix with {len(a)} elements.")
    return Mat

def convert_matrix2line(Matrix):
    result = []
    dim = len(Matrix) 
    for i in range(dim):
        for j in range(dim):
            result.append(Matrix[i,j])
    return np.array(result)


my_data = np.array([[0.21530249, 0.32450331, 0        ],
       [0.1930605 , 0.31788079, 0        ],
       [0.17793594, 0.31788079, 0        ],
       [0.16459075, 0.31125828, 1        ],
       [0.24822064, 0.31125828, 0        ],
       [0.28647687, 0.32450331, 0        ],
       [0.32829181, 0.31788079, 0        ],
       [0.38879004, 0.32450331, 0        ],
       [0.42882562, 0.32450331, 0        ],
       [0.47419929, 0.32450331, 0        ],
       [0.5044484 , 0.32450331, 0        ],
       [0.1797153 , 0.31125828, 0        ],
       [0.16548043, 0.31125828, 1        ],
       [0.17793594, 0.29801325, 1        ],
       [0.1930605 , 0.31788079, 0        ]])

Data = pd.DataFrame(my_data, columns=['X_1','X_2', 'Allum'])



def logLB(params,Data):
    B = convert_line2matrix(params)
    X = np.array(Data.iloc[:,:len(B)]) 
    Y = np.array(Data.iloc[:,len(B)])

    result = 0
    n = len(Data)
    BB = np.transpose(B) @ B
    for i in range(n):
        if(1-np.exp(-X[i].T @ BB @ X[i]) > 0):
            result += Y[i]*(-np.transpose(X[i]) @ BB @ X[i]) + (1 - Y[i])*np.log(1-np.exp(-X[i].T @ BB @ X[i]))

    return result

def f(params, Data):
    return -logLB(params, Data)



def dlogLB(params, Data):
    B = convert_line2matrix(params)
    X = np.array(Data.iloc[:,:len(B)]) 
    Y = np.array(Data.iloc[:,len(B)])
 
    BB = B.T @ B
    N = len(Data)
    M = len(B)
    Jacobian = np.zeros(np.shape(B))
    for n in range(len(B)):
        for m in range(len(B)):
            result = 0
            for c in range(N):
                som = 0
                for i in range(M):
                    som += X[c,m]*B[n,i]*X[c,i]
                if (1 - np.exp(-X[c].T @ BB @ X[c]) > 0):
                    result += -2*Y[c]*som + (1-Y[c])*np.exp(-X[c].T @ BB @ X[c])*(2*som)/(1 - np.exp(-X[c].T @ BB @ X[c]))

                Jacobian[n,m] = result 

    return convert_matrix2line(Jacobian)

def Jacf(params, Data):
    return -dlogLB(params, Data)

Answer 1

I assume that you want to stop the optimizer as soon as the euclidian norm of the gradient reaches a specific value, which is exactly the meaning of the BFGS method's gtol option. Otherwise, it doesn't make any sense mathematically, since the evaluated gradient is a vector and thus can't be compared to a scalar value.

The Newton-CG method doesn't provide a similar option. However, you could use a simple callback that is called after each iteration and terminates the algorithm when the callback returns True. Unfortunately, you can only terminate the optimizer by a callback with the trust-constr method. For all other methods, the callback's return value is ignored, so it's very limited.

A possible hacky and ugly way to terminate the optimizer by the callback anyway would be raising an exception:

import numpy as np
from scipy.optimize import minimize

class Callback:
    def __init__(self, eps, args, jac):
        self.eps = eps
        self.args = args
        self.jac = jac
        self.x = None
        self.gtol = None

    def __call__(self, xk):
        self.x = xk
        self.gtol = np.linalg.norm(self.jac(xk, *self.args))
        if self.gtol <= self.eps:
            raise Exception("Gradient norm is below threshold")

Here, xk is the current iterate, eps your desired tolerance, args a tuple containing your optional objective und gradient arguments and jac the gradient. Then, you can use it like this:

from scipy.optimize import minimize

cb = Callback(1.0e-1, (Data,), Jacf)

try:
    res = minimize(f, [5.0,1.0,2.0,5.0], args=Data, method='Newton-CG', 
    jac=Jacf, callback=cb)
except:
    x = cb.x
    gtol = cb.gtol
    print(f"gtol = {gtol:E}, x = {x}")

which yields

gtol = 5.515263E-02, x = [14.43322108 -5.18163542  0.22582261 -0.04859385]