Number of function calls and jacobian calls in scipy.root

Question

Just as an exercise, I am numerically solving the following system of equations:

$$ \begin{equation} \begin{cases} x^2 + y^2 = 32 \\ 3x + 7y = 15 \end{cases} \end{equation} $$

I use the following code:

from scipy.optimize import root

nfev = 0
def fun(x):
    global nfev
    nfev += 1
    return x[0]**2 + x[1]**2 - 32, 3 * x[0] + 7 * x[1] - 15

njev = 0
def jac(x):
    global njev
    njev += 1
    jacobian = [[2 * x[0], 2 * x[1]], [3, 7]]
    return jacobian

x0 = [0, 0]
result = root(fun, x0, jac=jac)

print(result)
print(nfev, njev)

I get the following output:

 message: The solution converged.
 success: True
  status: 1
     fun: [-1.269e-08 -7.468e-10]
       x: [ 1.500e+00 -5.454e+00]
    nfev: 13
    njev: 1
    fjac: [[-9.147e-02 -9.958e-01]
           [-9.958e-01  9.147e-02]]
       r: [-9.918e+00  9.915e-01  1.082e+01]
     qtf: [-2.960e-08  4.167e-07]
15 2

There are 15 calls to fun and 2 calls to jac which are different from 13 and 1 for respectively nfev and njev result fields.

Removing jac=jac, I get the following output:

 message: The solution converged.
 success: True
  status: 1
     fun: [-3.553e-15  0.000e+00]
       x: [ 1.500e+00  5.454e+00]
    nfev: 17
    fjac: [[-7.626e-02 -9.971e-01]
           [ 9.971e-01 -7.626e-02]]
       r: [-1.003e+01 -8.319e-01  1.088e+01]
     qtf: [ 6.485e-10 -8.479e-09]
19 0

There are 19 calls to fun which is different than nfev = 17 shown in the result.

What could explain the discrepancies? Would providing the Jacobian accelerate the calculations or improve the accuracy? If so in which circumstances?

Adding two print statements at the end of the jac function as so:

print(x)
print(jacobian)

I get the following extra output:

[0 0]
[[0, 0], [3, 7]]
[0. 0.]
[[0.0, 0.0], [3, 7]]

This shows that jac is called twice with the same x values. This does not make much sense. Any clue?

Answer 1

At the point where you print out the Jacobian, adding

traceback.print_stack()

reveals that the first evaluation comes from _check_func(). This function validates the shape of the results of the user-supplied functions before calling the minpack nonlinear system solver.

Answer 2

Studying your code and comparing it with the provided system of nonlinear equations, there is an error in the following code snippet.

def fun(x):
    global nfev
    nfev += 1
    return x[0]**2 + x[1]**2 - 32, 3 * x[0] + 7 * x[0] - 15

Correction

def fun(x):
    global nfev
    nfev += 1
                                                 #x[0] - 15
    return x[0]**2 + x[1]**2 - 32, 3 * x[0] + 7 * x[1] - 15

Nonlinear system solution

With the correction made, your code generates the correct approximate solutions for the given problem. Since we have a nonlinear system formed by a line and a circle centered at the origin with radius $r = \sqrt{32}$, it admits two real solutions that satisfy the problem. For this problem, I wrote a code in JavaScript applying the Newton-Raphson algorithm. An important observation regarding this numerical method lies in the fact that we calculate the inverse of the Jacobian matrix; consequently, it cannot be null.

"use strict";

// This function performs the subtraction of two vectors.
function subtractionVector(A,B) {
    if (A.length !== B.length) {
        console.log('Problem! The inserted matrices are not the same size.');
        process.exit(0);
    } else {
        let matrixOut = [];
        for (let i = 0; i < A.length; i++) {
            matrixOut[i] = A[i] - B[i];
        }
        return matrixOut;
    }
}


// Multiplies an nxm matrix by a mx1 vector
function multiplyMatrix(A,B) {
    if (A[0].length !== B.length) {
        console.log("Problem! The inserted matrices have different orders.");
        process.exit(0);
    } else {
        let matrixOut = [];
        for (let i = 0; i < A.length; i++) {
            matrixOut[i] = 0;
            for (let j = 0; j < A[i].length; j++) {
                matrixOut[i] += A[i][j] * B[j];
            }
        }
        return matrixOut;
    }
}

// Inverse matrix
function inverseMatrix(A) {
    const n = A.length;
    for (let k = 0; k < n; k++) {
        const con = A[k][k];
        A[k][k] = 1;
        for (let j = 0; j < n; j++) {
            A[k][j] = A[k][j] / con;
        }
        for (let i = 0; i < n; i++) {
            if (i !== k) {
                const con = A[i][k];
                A[i][k] = 0.0;
                for (let j = 0; j < n; j++) {
                    A[i][j] = A[i][j] - A[k][j] * con;
                }
            }
        }
    }
    return A;
}

// Norm the vector
function norm(vector) {
    let sumSquare = 0;
    for (let i = 0; i < vector.length; i++) {
        sumSquare += vector[i]**2;
    }
    return Math.sqrt(sumSquare);
}

// Norm
function normInfinity(vector) {
    function valorMaximo(A){
        let n = A.length;
        let max = A[0];
        for (let i = 1; i<=n-1; i++) {
            if(max<A[i]){
                max = A[i];
            }
        }
        return max;
    }

    let matrixA = [];
    for (let i = 0; i < vector.length; i++) {
        matrixA[i] = Math.abs(vector[i]);
    }
    return valorMaximo(matrixA);
}

// System equation f(r) = [f1(x,y,...),f2(x,y,...)]
function f(r) {
    let x = r[0];
    let y = r[1];

    // System of equations
    let f1 = x**2+y**2-32;
    let f2 = 3*x+7*y-15;

    return [f1,f2];
}

// Jacobian system
// Attention! The Jacobian cannot be null, because in the
// iterative Newton-Rapshon method we calculate the inverse
// of this, which is a square matrix.
function Jacobian(r) {
    let x = r[0];
    let y = r[1];

    // Partial derivative of f1(x,y) with respect to x and y
    let f1_x = 2*x; // Partial derivative of x
    let f1_y = 2*y; // Partial derivative of y

    // Partial derivative of f2(x,y) with respect to x and y
    let f2_x = 3; // Partial derivative of x
    let f2_y = 7; // Partial derivative of y

    return [[f1_x,f1_y],[f2_x,f2_y]];
}

function NewtonRapshon(r0,tol) {
    let counter = 0, x = r0, s = r0;
    do {
        s = x; // previous
        x = subtractionVector(x,multiplyMatrix(inverseMatrix(Jacobian(x)),f(x)));
        if(counter>1000){
            console.log('It didn\'t converge!');
            break;
        }
        //console.log(normInfinity(subtractionVector(x,s)))
        counter++;
    } while (norm(subtractionVector(x,s))>tol);
    for (let i = 0; i < x.length; i++) {
        x[i] = parseFloat(x[i].toFixed(3));
    }
    console.log(`\nNumber of iterations: ${counter}`);
    return x;
}

// Tolerance
let tol = 1.0e-15; // Tolerance

// First root
let r0_1 = [4,5]; // First kick
let sol_1 = NewtonRapshon(r0_1,tol);
console.log('Solution: ',sol_1);

// Second root
let r0_2 = [-5,-3]; // First kick
let sol_2 = NewtonRapshon(r0_2,tol);
console.log('Solution: ',sol_2);

Final result