SciPy - Non-Linear Least Squares

---

---

The Non-Linear least squares (NLLS) is a method for fitting a model to data where the model's parameters are non-linear. It minimizes the sum of squared residuals between the observed values and the model's predictions.

Mathematical equation of Non-Linear Least Squares for a set of residuals r(x) given a model f(x, t) is given as follows −

r(x)22

Where r(x) is the vector of residuals r_i(x) = y_i - f(x, t_i) and x is the vector of parameters to be estimated.

Implementation of Non-Linear Least Squares

SciPy provides the function scipy.optimize.least-squares() to solve non-linear least-squares problems efficiently. This function is flexible and supports various optimization methods and robust loss functions.

How Non-Linear Least Squares Method Works

Here are the steps how the Non-Linear least squares works in scipy −

• Residual Function: Defines the difference between the observed and predicted values for the model.
• Optimization Algorithms: SciPy uses methods like Trust Region Reflective ('trf'), Levenberg-Marquardt ('lm') and Dogleg ('dogbox') to find the optimal parameters.
• Robust Loss Functions: Functions such as 'soft_l1', 'huber' and others are used to reduce the influence of outliers in the fitting process.

Syntax

Following is the syntax for using the function scipy.optimize.least_squares() −

scipy.optimize.least_squares(
   fun, 
   x0, 
   jac='2-point', 
   bounds=(-inf, inf), 
   method='trf', 
   loss='linear', ...
)

Parameters

Here are the parameters of the function scipy.optimize.least_squares() −

• fun: The residual function that computes the difference between the observed data and the model.
• x0: Initial guess for the parameters to be optimized.
• jac (optional): Jacobian matrix or a function to compute it.
• bounds(optional): Constraints on the parameters given as lower_bounds, upper_bounds.
• method: Optimization algorithm ('trf', 'dogbox', 'lm').
• loss: Robust loss function to handle outliers ('linear', 'soft_l1', 'huber').

Fitting a Nonlinear Function

Lets solve a simple example where we fit a non-linear model y=aexp(bx) to some noisy data.

import numpy as np
from scipy.optimize import least_squares

# Generate synthetic data
x_data = np.linspace(0, 1, 10)
true_params = [2.5, -1.3]
y_data = true_params[0] * np.exp(true_params[1] * x_data) + 0.1 * np.random.randn(len(x_data))

# Define the residual function
def residuals(params, x, y):
    a, b = params
    return y - a * np.exp(b * x)

# Initial guess
x0 = [1.0, -1.0]

# Solve the least-squares problem
result = least_squares(residuals, x0, args=(x_data, y_data))

# Output results
print("Optimal parameters:", result.x)
print("Cost:", result.cost)
print("Residuals:", result.fun)

Following is the output of the above example which shows the optimal parameters, cost and residuals after fitting the model to the data −

Optimal parameters: [ 2.46335211 -1.22610415]
Cost: 0.0631726759601759
Residuals: [-0.14739481  0.13317732  0.12567623 -0.0355416  -0.12790527  0.19364081
 -0.04241043 -0.11512075 -0.02253795  0.02025557]

Constrained Nonlinear Fitting

Here is the example in which we can solve the same problem as mentioned above with bounds on the parameters such as a0, b-2 −

import numpy as np
from scipy.optimize import least_squares

# Generate synthetic data
x_data = np.linspace(0, 1, 10)
true_params = [2.5, -1.3]
y_data = true_params[0] * np.exp(true_params[1] * x_data) + 0.1 * np.random.randn(len(x_data))

# Define the residual function
def residuals(params, x, y):
    a, b = params
    return y - a * np.exp(b * x)

# Initial guess
x0 = [1.0, -1.0]

# Solve with bounds on the parameters
result = least_squares(residuals, x0, bounds=([0, -2], [5, 0]), args=(x_data, y_data))

# Output results
print("Optimal parameters with bounds:", result.x)

Following is the output of solving the non linear least squares with the bound parameter of the function scipy.optimize.least-squares() −

Optimal parameters with bounds: [ 2.4754014  -1.26138182]

Robust Nonlinear Fitting

Sometimes the data may contain outliers and we can use a robust loss function to reduce their impact. Lets solve the above example but introduce an outlier and use the 'soft_l1' loss function in it −

import numpy as np
from scipy.optimize import least_squares

# Generate synthetic data
x_data = np.linspace(0, 1, 10)
true_params = [2.5, -1.3]
y_data = true_params[0] * np.exp(true_params[1] * x_data) + 0.1 * np.random.randn(len(x_data))

# Define the residual function
def residuals(params, x, y):
    a, b = params
    return y - a * np.exp(b * x)

# Initial guess
x0 = [1.0, -1.0]

# Add an outlier to the data
y_data_with_outliers = y_data.copy()
y_data_with_outliers[2] += 1  # Introduce an outlier

# Solve using a robust loss function
result = least_squares(residuals, x0, args=(x_data, y_data_with_outliers), loss='soft_l1')

# Output results
print("Optimal parameters (robust):", result.x)

Here is the output of the robust Non-Linear least square fitting −

Optimal parameters (robust): [ 2.63552355 -1.16439492]