Optimization1805foundational9 min read

Nouvelles méthodes pour la détermination des orbites des comètes (Method of Least Squares)

طُرُق جديدة لتحديد مدارات المذنّبات (طريقة المربعات الصغرى)

Legendre, A.-M. — Firmin Didot, Paris

The problem

Astronomers tracking comets took measurements of position and time, but every observation carried — wobbling telescopes, atmospheric distortion, human imprecision. With more equations than unknowns (each measurement is one equation, but the orbit has only a few parameters), the system was inconsistent: no single set of orbital parameters satisfied all observations exactly. Before Legendre, there was no principled, general-purpose method for extracting the "best" parameters from an of noisy equations.

The contribution

Legendre proposed minimizing the : square each observation's (predicted minus observed), add them all up, and choose the parameters that make this total as small as possible. Setting the partial derivatives to zero yields the "" — a linear system that can be solved directly. This replaced ad-hoc judgment with an algebraic procedure any scientist could follow.

The impact

Least squares is the foundation of modern regression and a direct ancestor of the functions used to train every . The that minimizes in 2026 is the same sum-of-squared-residuals idea Legendre published in 1805 to fit comet orbits. Every time you hear "minimize the loss," you are hearing an echo of this paper.

Imagine you are a tailor making a suit for a customer you have never met. Five friends each guess the customer's chest measurement: 98 cm, 101 cm, 99 cm, 103 cm, 100 cm. No two agree. You cannot ask the customer directly.

Least squares says: pick the number that makes the total of all squared miss-distances as small as possible. That number turns out to be the average — 100.2 cm — and the suit fits better than trusting any single guess.

Legendre did the same thing, but instead of a chest measurement he was fitting an entire comet's orbit, and instead of five friends he had dozens of telescope readings.

The problem: too many measurements, no exact solution

By the late 1700s, astronomers could observe the position of a comet on several different nights. Each observation gave one equation relating the orbital parameters (the unknowns) to the measured angle and time. But there were always more observations than unknowns — the system was over-determined.

Worse, every observation contained measurement error. So the system was also inconsistent: no choice of parameters satisfied all the equations at once. The question was: when a perfect solution does not exist, what is the best compromise?

Before Legendre, astronomers relied on subjective judgment — picking two or three "trustworthy" observations, solving exactly, and ignoring the rest. Different astronomers got different answers from the same data. Science needed an objective, repeatable rule.

Open in Lab

Six observations (orange dots) cannot all lie on one line. Drag the line to see how the squared residuals (red squares) change. The least-squares line makes their total area smallest.

The demo wakes as you arrive…

The idea: minimize the sum of squared errors

Legendre's insight was elegant: define the "best" parameters as those that minimize a single number — the sum of squared residuals. A residual is the gap between what your predicts and what was actually observed. Square each gap (so that positive and negative errors don't cancel), add them all up, and choose the parameters that make this total as small as possible.

Why square the errors rather than just take their absolute values? Three practical reasons:

Squaring is differentiable everywhere — you can use calculus to find the minimum directly, whereas absolute values have a sharp corner at zero that complicates .
Squaring penalizes large errors more heavily — a single wildly wrong measurement dominates the sum, pushing the solution away from outliers. This reflects the intuition that a big error is not merely twice as bad as a small one, but disproportionately worse.
The resulting system of equations — the normal equations — is linear, so it can be solved with basic algebra, even by hand.

The formula: from intuition to algebra

Suppose you have $n$ observations. For each observation $i$ , the model predicts a value $\hat{y}_i$ that depends on parameters you can adjust, and the actual measured value is $y_i$ . The residual (error) is $e_i = y_i - \hat{y}_i$ . The objective is to minimize the total :

S = \sum_{i=1}^{n} e_i^2 = \sum_{i=1}^{n} (y_i - \hat{y}_i)^2

The least-squares objective — the number we want to minimize — Least squares measures how well a model's predictions match the observed data. Every prediction error contributes to the objective, but larger errors receive disproportionately greater penalties because they are squared. The goal is to find the model parameters that make the overall discrepancy between predictions and observations as small as possible. These parameter values define the best-fitting model under the least-squares criterion.

For a straight line $\hat{y}_i = a + bx_i$ , the objective becomes $S(a, b) = \sum (y_i - a - bx_i)^2$ . To find the minimum, take the with respect to each and set it to zero:

\frac{\partial S}{\partial a} = 0 \quad \Longrightarrow \quad \sum y_i = n \cdot a + b \sum x_i

Normal equation 1 — solving for the intercept

\frac{\partial S}{\partial b} = 0 \quad \Longrightarrow \quad \sum x_i y_i = a \sum x_i + b \sum x_i^2

Normal equation 2 — solving for the slope

These two equations are linear in $a$ and $b$ — they can be solved with pencil and paper. That was the practical breakthrough: a recipe any scientist could apply mechanically, regardless of the data.

Open in Lab

Enter data points and watch the normal equations form and solve in real time. Drag points to see how the optimal line changes.

The demo wakes as you arrive…

The modern form: matrices make it compact

In modern notation, stack all observations into a matrix equation $\mathbf{y} = \mathbf{X}\boldsymbol{\beta} + \mathbf{e}$ , where $\mathbf{X}$ is the (observations as rows, features as columns), $\boldsymbol{\beta}$ is the parameter vector, and $\mathbf{e}$ is the error vector. The normal equations become a single matrix equation:

\hat{\boldsymbol{\beta}} = (\mathbf{X}^\top \mathbf{X})^{-1} \mathbf{X}^\top \mathbf{y}

The closed-form least-squares solution — the "normal equation" — This equation provides a direct way to compute the parameters of a linear regression model without iterative optimization. It combines information about the relationships among the input features with information about how those features relate to the target values. The resulting parameter vector is the one that minimizes the total squared prediction error, making it the optimal solution under the least-squares criterion.

Read this formula as: "project the target $\mathbf{y}$ onto the of $\mathbf{X}$ ." Geometrically, the predicted values $\hat{\mathbf{y}} = \mathbf{X}\hat{\boldsymbol{\beta}}$ form the closest point in the column space to the observation vector — the . The residual vector $\mathbf{e} = \mathbf{y} - \hat{\mathbf{y}}$ is perpendicular to every column of $\mathbf{X}$ . That perpendicularity condition is the normal equations.

Open in Lab

The observation vector y (orange) is projected onto the column space of X (the green plane). The residual vector (red) is perpendicular — that's what the normal equations enforce.

The demo wakes as you arrive…

The idea in code

Least squares — both the closed-form and iterative approachespython

Simplified to show the idea — not the real implementation.

import numpy as np

# === Closed-form solution (the normal equation) ===
def least_squares_closed(X, y):
    """Solve β = (XᵀX)⁻¹Xᵀy directly."""
    return np.linalg.inv(X.T @ X) @ X.T @ y

# === Gradient descent approach (preview of modern ML) ===
def least_squares_gd(X, y, lr=0.01, steps=1000):
    """Minimize ‖y − Xβ‖² by gradient descent."""
    beta = np.zeros(X.shape[1])
    for _ in range(steps):
        residuals = y - X @ beta          # the errors
        gradient = -2 * X.T @ residuals   # direction of steepest ascent
        beta -= lr * gradient              # step opposite to gradient
    return beta

# --- Demo: fit a line y = a + b·x to noisy data ---
np.random.seed(1805)  # Legendre's year
x = np.linspace(0, 10, 20)
y_true = 2.0 + 0.5 * x
y_noisy = y_true + np.random.normal(0, 0.5, size=len(x))

# Design matrix: column of ones (intercept) + column of x values
X = np.column_stack([np.ones_like(x), x])

beta_closed = least_squares_closed(X, y_noisy)
beta_gd     = least_squares_gd(X, y_noisy)

# Both converge to the same answer:
# beta ≈ [2.0, 0.5] — the true intercept and slope
print(f"Closed-form: a={beta_closed[0]:.3f}, b={beta_closed[1]:.3f}")
print(f"Grad descent: a={beta_gd[0]:.3f}, b={beta_gd[1]:.3f}")

The bridge to machine learning: from comet orbits to neural networks

The core idea of least squares — define an error measure, then minimize it — is the template for all . When you train a neural network:

The (usually mean squared error for regression) is Legendre's sum of squared residuals, divided by $n$ .
is the iterative cousin of the normal equations: instead of solving $\nabla S = 0$ in one shot, it takes small steps downhill. The normal equations work only for linear models; gradient descent works for any differentiable model, including deep networks.
Backpropagation is the applied to compute $\nabla S$ efficiently through layers — but $S$ is still the same squared-error sum that Legendre invented.

The intellectual lineage is direct: Legendre (1805) → Gauss (1809, probabilistic justification) → gradient descent (Cauchy, 1847) → stochastic gradient descent (Robbins & Monro, 1951) → backpropagation (Rumelhart, Hinton, Williams, 1986) → every modern neural network.

Open in Lab

See the same squared-error idea evolve from Legendre's pen-and-paper solution to modern gradient descent. Toggle between the closed-form and iterative approaches.

The demo wakes as you arrive…

The priority dispute: Legendre vs. Gauss

Legendre published the method of least squares in 1805 in a nine-page appendix to this very book on comet orbits. Four years later, Carl Friedrich Gauss published his own work on celestial mechanics and claimed he had been using the method privately since 1795 — a full decade before Legendre's publication — but offered no written proof from that period.

Legendre was furious. In a letter to Gauss he wrote that a claim of priority without prior publication was, at the very least, inappropriate. Modern historians generally credit Legendre with the first publication and Gauss with the probabilistic justification: Gauss showed in 1809 that if measurement errors follow a Gaussian (normal) distribution, then least squares gives the most probable values of the unknowns — the maximum likelihood estimate. Laplace later strengthened this connection in 1810.

The irony: the bell-shaped curve we now call the "" is itself named after Gauss partly because of this least-squares connection, even though the distribution was studied earlier by de Moivre. Names in mathematics often reward the person who reveals the deepest connection, not the first discoverer.

1805
Legendre publishes: first description of least squares
A nine-page appendix in "Nouvelles méthodes pour la détermination des orbites des comètes" presents the method and applies it to the shape of the Earth.
1809
Gauss: probabilistic justification (maximum likelihood)
In "Theoria motus," Gauss shows that least squares gives the maximum likelihood estimate when errors are normally distributed, and claims prior independent use since 1795.
1810
Laplace: central limit theorem strengthens the case
Laplace proves the central limit theorem, showing that averages of many errors tend to be normally distributed — thus broadening the conditions under which least squares is optimal.
1847
Cauchy: gradient descent
Cauchy introduces the method of steepest descent — an iterative alternative to solving the normal equations, enabling optimization of non-linear objectives.
1951
Robbins & Monro: stochastic approximation
Stochastic approximation provides the theoretical foundation for stochastic gradient descent — using one noisy sample at a time rather than all the data.
1986
Backpropagation: least squares meets deep networks
Rumelhart, Hinton & Williams apply the chain rule to compute gradients of the squared-error loss through multi-layer networks, unlocking deep learning.
2026
Today — every neural network minimizes a loss
GPT, Claude, and every deep model train by minimizing a loss function. The template — define error, square it, minimize — is Legendre's 1805 idea, industrialized.

CitationLegendre, A.-M.. Nouvelles méthodes pour la détermination des orbites des comètes. Firmin Didot, Paris, 1805.

Terms in this paper

Squared errorخطأ تربيعي
Normal Equationsالمعادلات السوية
Over-Determined Systemنظام مُفرَط التحديد
Residualالبُعد المتبقّي
Closed-Form Solutionحل بصيغة مغلقة
Orthogonal Projectionالإسقاط العمودي
Design Matrixمصفوفة التصميم
Sum of Squared Errorsمجموع مربعات الأخطاء