I was working on a data analysis project where I needed to fit a complex exponential decay model to experimental measurements. The challenge was that my model was non-linear, making traditional linear regression methods unsuitable. This is where SciPy’s leastsq function came to the rescue.
In this article, I’ll share how to use SciPy’s leastsq function to solve nonlinear least squares problems in Python. Whether you’re fitting complex models to experimental data or optimizing parameters in scientific computing, this useful tool can help you get accurate results with minimal effort.
So let’s start..
SciPy’s Leastsq Function
SciPy’s leastsq function is part of the scipy.optimize module and provides a way to solve nonlinear least squares problems. It’s essentially used to find the parameters that minimize the sum of squares of a set of functions.
In simpler terms, when you have a model with unknown parameters and you want to find the parameter values that best fit your data, leastsq can help you do that, even when your model is non-linear.
Basic Usage of Leastsq
Let’s start with a simple example of how to use the leastsq function:
import numpy as np
from scipy.optimize import leastsq
import matplotlib.pyplot as plt
# Generate some sample data with noise
def true_func(x):
return 2.5 * np.exp(-x * 0.5) + 0.5
# Create noisy data
x_data = np.linspace(0, 5, 50)
y_true = true_func(x_data)
y_noise = 0.2 * np.random.normal(size=len(x_data))
y_data = y_true + y_noise
# Define the function to minimize (residuals)
def residuals(params, x, y):
a, b, c = params
model = a * np.exp(-b * x) + c
return model - y
# Initial guess for parameters
initial_guess = [1.0, 1.0, 1.0]
# Perform the fit
result, success = leastsq(residuals, initial_guess, args=(x_data, y_data))
# Extract the optimized parameters
a_opt, b_opt, c_opt = result
# Plot the results
plt.figure(figsize=(10, 6))
plt.scatter(x_data, y_data, label='Noisy data')
plt.plot(x_data, true_func(x_data), 'r-', label='True function')
plt.plot(x_data, a_opt * np.exp(-b_opt * x_data) + c_opt, 'g--', label='Fitted function')
plt.legend()
plt.xlabel('X')
plt.ylabel('Y')
plt.title('Exponential Decay Fitting with Leastsq')
plt.show()
print(f"Optimized parameters: a={a_opt:.4f}, b={b_opt:.4f}, c={c_opt:.4f}")
print(f"True parameters: a=2.5000, b=0.5000, c=0.5000")You can see the output in the screenshot below.
In this example, we:
- Generate synthetic data following an exponential decay with added noise
- Define a residual function that computes the difference between our model and the data
- Provide initial parameter guesses
- Use leastsq to find the optimal parameters
- Plot the results to visualize the fit
Advanced Usage: Working with Full Jacobian
For more complex problems, you might want to provide the Jacobian matrix of the residual function. This can significantly improve the performance and accuracy of the fitting process.
import numpy as np
from scipy.optimize import leastsq
import matplotlib.pyplot as plt
# Generate data for a damped oscillation
def true_func(x, params):
a, w, phi, c = params
return a * np.exp(-0.1 * x) * np.cos(w * x + phi) + c
# True parameters: amplitude, frequency, phase, offset
true_params = [2.0, 1.5, 0.3, 0.5]
# Create sample data with noise
x_data = np.linspace(0, 10, 100)
y_true = true_func(x_data, true_params)
y_noise = 0.2 * np.random.normal(size=len(x_data))
y_data = y_true + y_noise
# Define the residual function
def residuals(params, x, y):
a, w, phi, c = params
model = a * np.exp(-0.1 * x) * np.cos(w * x + phi) + c
return model - y
# Define the Jacobian (partial derivatives with respect to each parameter)
def jacobian(params, x, y):
a, w, phi, c = params
exp_term = np.exp(-0.1 * x)
cos_term = np.cos(w * x + phi)
sin_term = np.sin(w * x + phi)
# Partial derivatives
da = exp_term * cos_term
dw = -a * exp_term * x * sin_term
dphi = -a * exp_term * sin_term
dc = np.ones_like(x)
return np.vstack([da, dw, dphi, dc]).T
# Initial guess
initial_guess = [1.0, 1.0, 0.0, 0.0]
# Perform the fit with Jacobian
result, success = leastsq(residuals, initial_guess, args=(x_data, y_data),
Dfun=jacobian, col_deriv=0, full_output=1)
# Extract the optimized parameters
params_opt = result[0]
a_opt, w_opt, phi_opt, c_opt = params_opt
# Plot the results
plt.figure(figsize=(12, 6))
plt.scatter(x_data, y_data, alpha=0.6, label='Noisy data')
plt.plot(x_data, true_func(x_data, true_params), 'r-', label='True function')
plt.plot(x_data, true_func(x_data, params_opt), 'g--', label='Fitted function')
plt.legend()
plt.xlabel('Time (s)')
plt.ylabel('Amplitude')
plt.title('Damped Oscillation Fitting with Jacobian')
plt.grid(True)
plt.show()
print("True parameters:", true_params)
print("Fitted parameters:", [round(p, 4) for p in params_opt])You can see the output in the screenshot below.
By providing the Jacobian, we give the algorithm analytical derivatives instead of having it estimate them numerically. This makes the fitting process more accurate and efficient, especially for complex models.
Check out Python SciPy Confidence Interval
Common Challenges Faced While Working with SciPy’s Leastsq Function and Solutions
Now, I will explain to you the common challenges that can be faced while working with SciPy’s Leastsq Function.
1. Convergence Issues
If your fitting process doesn’t converge, try these approaches:
# Improve convergence with better initial guesses
# Instead of:
initial_guess = [1.0, 1.0, 1.0]
# Try a more informed guess based on data characteristics:
amplitude_guess = np.max(y_data) - np.min(y_data)
offset_guess = np.mean(y_data)
decay_guess = 0.5 # Domain knowledge
initial_guess = [amplitude_guess, decay_guess, offset_guess]
# You can also try different method parameters
result, success = leastsq(residuals, initial_guess, args=(x_data, y_data),
ftol=1e-10, xtol=1e-10, maxfev=10000)2. Parameter Bounds
The leastsq function doesn’t directly support parameter bounds. If you need bounded optimization, consider using scipy.optimize.least_squares (a newer alternative):
from scipy.optimize import least_squares
# Define bounds (lower and upper bounds for each parameter)
lower_bounds = [0, 0, 0] # All parameters must be positive
upper_bounds = [10, 5, 2] # Upper limits for each parameter
# Use least_squares instead of leastsq
result = least_squares(residuals, initial_guess,
args=(x_data, y_data),
bounds=(lower_bounds, upper_bounds))
# Extract optimized parameters
params_opt = result.x3. Weight Data Points
Sometimes, not all data points are equally reliable. You can add weights to give more importance to certain measurements:
# Define weights (higher weight = more important)
weights = 1.0 / np.array([0.1, 0.2, 0.1, 0.3, 0.2, 0.1, 0.2, 0.3, 0.2, 0.1])
# Modify the residual function to include weights
def weighted_residuals(params, x, y, weights):
model = some_model_function(x, params)
return weights * (model - y)
# Use the weighted residuals in leastsq
result, success = leastsq(weighted_residuals, initial_guess,
args=(x_data, y_data, weights))Check out Python SciPy Freqz
Real-World Example: Temperature Sensor Calibration
Let’s look at a more practical example. Imagine you’re working for a weather station in New York City, and you need to calibrate temperature sensors. You have measurements from a reference sensor and your new sensor that needs calibration.
import numpy as np
from scipy.optimize import leastsq
import matplotlib.pyplot as plt
# Reference temperature readings (°F) in NYC throughout a day
reference_temps = np.array([55.2, 57.8, 60.5, 65.3, 72.1, 75.6, 78.4,
79.2, 77.5, 74.3, 70.1, 65.8, 60.2])
# Your uncalibrated sensor readings for the same times
uncalibrated_temps = np.array([54.0, 56.3, 58.9, 63.8, 70.5, 73.9, 76.5,
77.0, 75.2, 72.1, 68.0, 63.5, 58.1])
# Define the calibration model: reference = a * uncalibrated + b
def residuals(params, uncalibrated, reference):
a, b = params
return a * uncalibrated + b - reference
# Initial guess for calibration parameters
initial_guess = [1.0, 0.0] # Start with no adjustment
# Perform the calibration
result, success = leastsq(residuals, initial_guess, args=(uncalibrated_temps, reference_temps))
# Extract the optimized parameters
a_opt, b_opt = result
# Apply calibration to the uncalibrated data
calibrated_temps = a_opt * uncalibrated_temps + b_opt
# Plot the results
plt.figure(figsize=(10, 6))
plt.scatter(uncalibrated_temps, reference_temps, label='Uncalibrated vs Reference')
plt.plot([min(uncalibrated_temps), max(uncalibrated_temps)],
[min(uncalibrated_temps), max(uncalibrated_temps)], 'k--', label='Perfect calibration')
plt.plot(uncalibrated_temps, calibrated_temps, 'r-', label='Calibrated')
plt.legend()
plt.xlabel('Sensor Temperature (°F)')
plt.ylabel('Reference Temperature (°F)')
plt.title('Temperature Sensor Calibration')
plt.grid(True)
plt.show()
print(f"Calibration parameters: a={a_opt:.4f}, b={b_opt:.4f}")
print(f"This means: Reference = {a_opt:.4f} × Sensor + {b_opt:.4f}")You can see the output in the screenshot below.
In this example, we calibrated a temperature sensor by finding the linear relationship between the uncalibrated readings and reference values, which is a common task in instrumentation and measurement science.
Read Python SciPy Stats Multivariate_Normal
Leastsq vs. Curve_fit: Which One to Choose?
SciPy provides another function called curve_fit that’s built on top of leastsq. It’s often easier to use for common curve fitting tasks:
from scipy.optimize import curve_fit
# Define your model function
def model_func(x, a, b, c):
return a * np.exp(-b * x) + c
# Fit using curve_fit
params_opt, params_cov = curve_fit(model_func, x_data, y_data)
# params_opt contains the optimized parameters
# params_cov gives the covariance matrix (useful for error estimation)
So, when should you use each?
- Use
curve_fitwhen:- You have a simple model function
- You want parameter error estimates
- You prefer a simpler interface
- Use
leastsqwhen:- You need more control over the optimization process
- Your residual function is complex
- You want to provide an analytical Jacobian
- You’re solving a system of equations, not just fitting a curve
Error Analysis with Leastsq
Understanding the uncertainty in your fitted parameters is crucial for scientific work:
import numpy as np
from scipy.optimize import leastsq
from scipy import stats
# After running leastsq with full_output=True
result = leastsq(residuals, initial_guess, args=(x_data, y_data), full_output=True)
params_opt, cov_x, infodict, mesg, ier = result
# Number of data points and parameters
n = len(x_data)
p = len(initial_guess)
# Residual sum of squares
residual = residuals(params_opt, x_data, y_data)
ss_res = np.sum(residual**2)
# Degrees of freedom
dof = max(0, n - p)
# Compute reduced chi-squared
chi_squared = ss_res / dof
# Scale covariance matrix
if cov_x is not None:
cov_x = cov_x * chi_squared
# Standard errors
param_errors = np.sqrt(np.diag(cov_x))
# Confidence intervals (95%)
t_value = stats.t.ppf(0.975, dof)
conf_intervals = np.array([(params_opt[i] - t_value * param_errors[i],
params_opt[i] + t_value * param_errors[i])
for i in range(p)])
print("Parameter values and 95% confidence intervals:")
for i, (param, error, interval) in enumerate(zip(params_opt, param_errors, conf_intervals)):
print(f"p{i} = {param:.4f} ± {error:.4f} ({interval[0]:.4f}, {interval[1]:.4f})")
This code calculates standard errors and confidence intervals for your fitted parameters, giving you a measure of how reliable your results are.
Best Practices for Using Leastsq
Based on my years of experience working with SciPy’s optimization tools, here are some best practices:
- Start with reasonable initial guesses: The closer your initial guess is to the solution, the faster and more reliably leastsq will converge.
- Normalize your data: If your x and y variables have very different scales, normalize them before fitting to improve numerical stability.
- Provide analytical Jacobians when possible: This improves both speed and accuracy, especially for complex models.
- Check the return status: The second output from leastsq tells you if the optimization was successful.
- Plot your results: Always visualize your fits to ensure they make sense.
- Experiment with different models: Sometimes, a simpler model can give more robust results than a complex one.
- Use bounds when needed: If you need parameter constraints, consider using
least_squaresinstead ofleastsq.
Check out Python SciPy Eigenvalues
Conclusion
SciPy’s leastsq function is a useful tool for solving nonlinear least squares problems in Python. I’ve used it extensively in my work, from fitting complex physical models to calibrating sensors and analyzing experimental data.
While it may take some practice to master, the flexibility and power it offers make it worth the effort, especially for scientific and engineering applications where accurate parameter estimation is crucial.
You may also read other SciPy-related articles.
I am Bijay Kumar, a Microsoft MVP in SharePoint. Apart from SharePoint, I started working on Python, Machine learning, and artificial intelligence for the last 5 years. During this time I got expertise in various Python libraries also like Tkinter, Pandas, NumPy, Turtle, Django, Matplotlib, Tensorflow, Scipy, Scikit-Learn, etc… for various clients in the United States, Canada, the United Kingdom, Australia, New Zealand, etc. Check out my profile.
