Jacobian matrix produced by least_squares contains inverse signs

Question

I have been working with scipy least_squares and working to be able to produce the different outputs by hand. When I get to the jacobian matrix, I am able to produce it by hand, but the values I get are positive while the matrix produced by python contains negative values. Otherwise, the values are an exact match. Does anyone know what causes this?

import numpy as np
from scipy.optimize import least_squares
import pandas as pd
from sympy import Matrix, pprint
import matplotlib.pyplot as plt

# Define your model function
def model_func(x, params):
 return params[0] * x ** 2 + params[1] * x + params[2]
    
# Define the function to compute the residuals
def residuals(params, x_data, y_data):
 return y_data - model_func(x_data, params)
    
# Sample data
x_data = np.array([1, 2, 3, 4, 5])
y_data = np.array([2, 4, 9, 16, 25])

# Initial guess for parameters
initial_params = [1, 1, 1]

# Perform the least squares optimization
result = least_squares(residuals, initial_params, args=(x_data, y_data))

# Results
x = result.x # Optimized parameters
resnorm = result.cost # Sum of residuals squared
residual = result.fun # Array of residuals
exitflag = result.status # Status of the solver
output = result.message # Description of the cause of the termination
jacobian = result.jac # Jacobian matrix at the solution

print("Optimized parameters:", x)
print("Residual norm:", resnorm)
print("Residuals:", residual)
print("Exit flag:", exitflag)
print("Output message:", output)
print("Jacobian:", jacobian)

And returns the following:

When I compute it by hand I get:

I have looked everywhere it feels like to try and understand why this is the case, anyone have any ideas?

Answer 1

Provided the function you are analyzing is (which is not the function in your code):

def func(x, p):
    return p[0]  * x **2 + p[1] * x + p[2]

The result of the optimization is:

     message: `xtol` termination condition is satisfied.
     success: True
      status: 3
         fun: [ 1.143e-01 -2.571e-01  8.571e-02  1.429e-01 -8.571e-02]
           x: [ 1.143e+00 -1.057e+00  1.800e+00]
        cost: 0.0571428571428573
         jac: [[-1.000e+00 -1.000e+00 -1.000e+00]
               [-4.000e+00 -2.000e+00 -1.000e+00]
               ...
               [-1.600e+01 -4.000e+00 -1.000e+00]
               [-2.500e+01 -5.000e+00 -1.000e+00]]
        grad: [ 8.941e-09  1.611e-09  1.608e-08]
  optimality: 1.6083786736995895e-08
 active_mask: [ 0.000e+00  0.000e+00  0.000e+00]
        nfev: 5
        njev: 4

Where the matrix called jac effectively have been multiplied by a negative sign. If your read the definition of jac in the documentation, you will find that this is not exactly the Jacobian but a modified version of it:

Modified Jacobian matrix at the solution, in the sense that J^T J is a Gauss-Newton approximation of the Hessian of the cost function. The type is the same as the one used by the algorithm.

That is, the name is a bit misleading and hide an implementation detail.