Least squares (least_squares) in SciPy - Deep Dive

Overview - Least squares (least_squares)

What is it?

Least squares is a method to find the best fit solution to a system of equations that may not have an exact answer. It works by minimizing the sum of the squares of the differences between observed values and the values predicted by a model. The scipy library provides a function called least_squares to solve these problems efficiently. This method is widely used in data fitting, regression, and optimization.

Why it matters

Without least squares, we would struggle to find good approximations when data is noisy or when exact solutions don't exist. It helps us make sense of imperfect data by finding the closest possible match. This is crucial in fields like science, engineering, and economics where measurements have errors. Least squares turns messy real-world data into useful insights and predictions.

Where it fits

Before learning least squares, you should understand basic algebra, functions, and simple optimization concepts. After mastering least squares, you can explore advanced regression techniques, machine learning models, and nonlinear optimization methods.

Mental Model

Core Idea

Least squares finds the best solution by making the total squared difference between predicted and actual values as small as possible.

Think of it like...

Imagine trying to draw a straight line through a scatter of points on paper. Least squares is like adjusting the line so that the total squared distance from all points to the line is as small as possible, giving the best overall fit.

Observed points: •  •   •  •  •
Fitted line:  ─────────────
Differences: | |  |  |  |
Sum of squares: Σ(differences²) minimized

Build-Up - 7 Steps

1

Foundation

Understanding residuals and errors

Concept: Introduce the idea of residuals as differences between observed and predicted values.

When we try to predict data points using a model, the difference between the actual data and the model's prediction is called a residual. Residual = observed value - predicted value. These residuals show how far off our model is for each point.

Result

You can calculate residuals for any model and data, which measure prediction errors.

Understanding residuals is key because least squares works by minimizing these errors to improve the model.

2

Foundation

Why square residuals?

Concept: Explain why residuals are squared before summing.

If we just add residuals, positive and negative errors cancel out. Squaring residuals makes all errors positive and emphasizes larger errors more. This helps find a solution that balances all errors fairly.

Result

Sum of squared residuals is always positive and highlights big mistakes.

Squaring residuals prevents error cancellation and focuses on reducing large mistakes, which improves model accuracy.

3

Intermediate

Linear least squares with scipy

🤔Before reading on: do you think least_squares can solve linear problems directly or only nonlinear ones? Commit to your answer.

Concept: Learn how to use scipy's least_squares function for simple linear models.

You can define a function that calculates residuals for a linear model y = mx + b. Then use scipy.optimize.least_squares to find m and b that minimize residuals. Example: import numpy as np from scipy.optimize import least_squares def residuals(params, x, y): m, b = params return m * x + b - y x = np.array([1, 2, 3, 4]) y = np.array([2, 4, 6, 8]) result = least_squares(residuals, x0=[0, 0], args=(x, y)) print(result.x) # slope and intercept

Result

Output shows the best slope and intercept fitting the data.

Knowing how to set up residual functions and call least_squares unlocks flexible fitting for many models.

4

Intermediate

Handling nonlinear models

🤔Before reading on: do you think least_squares can handle models where parameters appear inside nonlinear functions? Commit to yes or no.

Concept: Extend least squares to nonlinear models where parameters affect the model in complex ways.

Least squares can fit models like y = a * exp(b * x) by defining residuals accordingly. For example: def residuals(params, x, y): a, b = params return a * np.exp(b * x) - y This flexibility allows fitting curves, growth models, and more.

Result

The function finds parameters a and b that best fit the nonlinear curve to data.

Understanding that least_squares works beyond lines opens doors to modeling real-world complex relationships.

5

Intermediate

Using bounds and constraints

🤔Before reading on: do you think least_squares can restrict parameter values within limits? Commit to yes or no.

Concept: Learn how to limit parameter values to realistic ranges during fitting.

least_squares supports bounds to keep parameters within specified intervals. For example: result = least_squares(residuals, x0=[1, 1], bounds=([0, 0], [10, 10]), args=(x, y)) This prevents nonsensical values like negative rates or impossible constants.

Result

Fitted parameters respect the bounds, improving model realism.

Knowing how to apply bounds helps avoid unrealistic fits and improves model trustworthiness.

6

Advanced

Jacobian and optimization speed

🤔Before reading on: do you think providing the Jacobian matrix speeds up least_squares? Commit to yes or no.

Concept: Learn about the Jacobian matrix and how supplying it can make fitting faster and more accurate.

The Jacobian is a matrix of partial derivatives of residuals with respect to parameters. If you provide a function that calculates this matrix, least_squares uses it to optimize more efficiently. Example: def jacobian(params, x, y): m, b = params return np.vstack((x, np.ones_like(x))).T result = least_squares(residuals, x0=[0, 0], jac=jacobian, args=(x, y))

Result

Optimization converges faster and more reliably.

Understanding Jacobians reveals how optimization algorithms use gradient information to improve performance.

7

Expert

Robust loss functions and outliers

🤔Before reading on: do you think least_squares always uses simple squared errors? Commit to yes or no.

Concept: Explore how least_squares can use different loss functions to reduce the effect of outliers.

By default, least_squares minimizes squared residuals, which can be sensitive to outliers. You can specify loss='soft_l1' or 'huber' to reduce outlier impact. For example: result = least_squares(residuals, x0=[0, 0], args=(x, y), loss='huber') This makes fitting more robust when data has errors or extreme points.

Result

Fitted parameters are less influenced by outliers, improving model reliability.

Knowing about robust loss functions helps build models that work well with real, messy data.

Under the Hood

least_squares uses iterative optimization algorithms like the Levenberg-Marquardt or Trust Region Reflective methods. It starts with initial guesses for parameters and repeatedly updates them to reduce the sum of squared residuals. Internally, it calculates residuals and optionally Jacobians to guide the search. The process stops when improvements become very small or a maximum number of iterations is reached.

Why designed this way?

The design balances flexibility and efficiency. Iterative methods handle both linear and nonlinear problems. Providing Jacobians speeds convergence but is optional for ease of use. Bounds and loss functions add robustness for real-world data. Alternatives like direct matrix inversion are limited to linear problems and less stable with noise.

Initial guess
    ↓
Calculate residuals and Jacobian
    ↓
Update parameters (optimization step)
    ↓
Check convergence
   ┌─────────────┐
   │ Not converged? ├─Yes─> Repeat
   └─────────────┘
    ↓ No
Return best parameters

Myth Busters - 4 Common Misconceptions

Quick: Does least_squares always find the global best solution? Commit yes or no.

Common Belief: Least squares always finds the perfect best fit solution.

Reality: least_squares finds a local minimum near the initial guess, which may not be the global best solution, especially for nonlinear problems.

Why it matters: Relying on a single run can lead to suboptimal fits; multiple initial guesses or methods may be needed.

Quick: Is least_squares only for linear models? Commit yes or no.

Common Belief: Least squares only works for straight-line or linear models.

Reality: least_squares can fit nonlinear models by minimizing residuals defined by any function.

Why it matters: Limiting least squares to linear models misses its power for complex real-world data.

Quick: Does squaring residuals always improve model accuracy? Commit yes or no.

Common Belief: Squaring residuals always leads to the best model fit.

Reality: Squaring residuals can overly penalize outliers, sometimes harming model robustness.

Why it matters: Ignoring this can cause poor fits when data contains errors or extreme values.

Quick: Can you ignore bounds when fitting parameters? Commit yes or no.

Common Belief: Bounds are optional and rarely affect results.

Reality: Bounds can be crucial to keep parameters realistic and prevent nonsensical fits.

Why it matters: Ignoring bounds can produce invalid models that mislead decisions.

Expert Zone

1

Providing an accurate Jacobian can drastically reduce computation time and improve convergence stability.

2

Choosing the right loss function is critical for handling outliers and noisy data effectively.

3

Initial parameter guesses strongly influence the solution in nonlinear least squares, requiring domain knowledge or heuristics.

When NOT to use

least_squares is not ideal for very large datasets where stochastic or batch optimization methods like gradient descent are more efficient. Also, if the model is not differentiable or residuals are not smooth, alternative optimization methods should be considered.

Production Patterns

In production, least_squares is often wrapped in pipelines that preprocess data, validate parameter bounds, and run multiple fits with different initial guesses. It is combined with robust loss functions and automated diagnostics to ensure reliable model deployment.

Connections

Linear Regression

least_squares generalizes the core idea of linear regression by allowing nonlinear models and constraints.

Understanding least squares deepens comprehension of regression as an optimization problem, not just a formula.

Gradient Descent Optimization

least_squares uses gradient-based iterative methods similar to gradient descent to minimize errors.

Knowing least squares helps grasp how gradient information guides parameter updates in many machine learning algorithms.

Physics: Experimental Data Fitting

least_squares is widely used in physics to fit models to experimental measurements with noise.

Recognizing least squares in physics experiments shows how math tools translate raw data into scientific laws.

Common Pitfalls

#1 Ignoring initial parameter guesses for nonlinear problems.

Wrong approach:

result = least_squares(residuals, x0=[0, 0], args=(x, y)) # no thought to starting values

Correct approach:

result = least_squares(residuals, x0=[1, 0.5], args=(x, y)) # informed initial guess

Root cause:

Nonlinear optimization depends on starting points; poor guesses lead to bad local minima.

#2 Not using bounds when parameters must be positive.

Wrong approach:

result = least_squares(residuals, x0=[1, 1], args=(x, y)) # no bounds

Correct approach:

result = least_squares(residuals, x0=[1, 1], bounds=([0, 0], [np.inf, np.inf]), args=(x, y))

Root cause:

Without bounds, optimization can produce negative or invalid parameter values.

#3 Assuming least_squares handles outliers well by default.

Wrong approach:

result = least_squares(residuals, x0=[0, 0], args=(x, y)) # default loss

Correct approach:

result = least_squares(residuals, x0=[0, 0], args=(x, y), loss='huber') # robust loss

Root cause:

Default squared loss is sensitive to outliers, which can skew results.

Key Takeaways

Least squares finds the best fit by minimizing the sum of squared differences between predicted and observed data.

It works for both linear and nonlinear models by defining residual functions and using iterative optimization.

Providing bounds and robust loss functions improves model realism and resilience to outliers.

Jacobian matrices speed up optimization but are optional for ease of use.

Initial guesses and understanding optimization limits are crucial for reliable results.