Nonlinear Least Squares

In the quest to understand the world, scientists and engineers constantly seek to connect theoretical models with experimental facts. This often involves a fundamental challenge: finding the specific set of parameters that makes a mathematical model best describe a collection of data points. While simple linear relationships can be solved directly, many natural and engineered systems are governed by complex, nonlinear equations. This creates a significant knowledge gap: how do we reliably fit these intricate models to reality?

This article introduces Nonlinear Least Squares (NLLS), a powerful and versatile method designed to solve precisely this problem. It is the quantitative engine that allows us to test our theories against data, turning scattered measurements into concrete knowledge. Across the following chapters, you will gain a comprehensive understanding of this essential technique. First, we will explore the "Principles and Mechanisms," dissecting what makes a problem nonlinear and examining the clever iterative algorithms, like Gauss-Newton and Levenberg-Marquardt, that navigate the complex error landscape to find a solution. Following that, in "Applications and Interdisciplinary Connections," we will journey through diverse fields—from astronomy to biology and cryptography—to witness how this single mathematical idea provides a common language for scientific discovery.

Imagine you are a detective, and you have a scattering of clues—a set of data points. Your goal is to uncover the underlying law, the mathematical model, that connects these clues. The method of least squares is one of the most powerful tools in your detective kit. Its principle is beautifully simple: the best model is the one that minimizes the total "error," defined as the sum of the squared differences between your observed data and the model's predictions.

If we have a set of data points $(x_i, y_i)$ and a model function $f(x; \mathbf{p})$ with parameters $\mathbf{p}$, we are trying to find the specific parameters $\mathbf{p}$ that make the model "hug" the data as closely as possible. We quantify this "closeness" with the sum-of-squares error function, $S$:

For instance, if we were modeling the decaying oscillation of a mechanical system, our model might be $f(t; p_1, p_2, p_3) = p_1 \exp(-p_2 t) \cos(p_3 t)$. Our task would be to find the initial amplitude $p_1$, damping factor $p_2$, and frequency $p_3$ that minimize the sum of squared differences between this function and our measured displacements over time. This minimization problem is the heart of all least squares fitting. But how we solve it depends critically on the nature of the function $f$.

You might think that the dividing line between "easy" and "hard" fitting problems is the shape of the model's curve. A straight line is easy, a curve is hard. Nature, however, has a more subtle distinction. The crucial property is not whether the model is linear in its variable (like $x$), but whether it is ​​linear in its parameters​​ (the coefficients we are trying to find).

A model is linear in its parameters if it can be written as a simple sum of those parameters, each multiplied by a function of the independent variable $x$. That is, $f(x; \mathbf{c}) = c_1 g_1(x) + c_2 g_2(x) + \dots + c_k g_k(x)$. The functions $g_j(x)$ can be as wildly nonlinear as you like—$\sin(x)$, $\ln(x)$, $x^2$, anything—as long as they don't contain the parameters $c_j$ we are solving for.

Consider these two models:

1. $f(x; c_1, c_2) = c_1 x + c_2 x^2$
2. $f(x; c_1, c_2) = c_1 (x - c_2)^2$

The first model, a simple quadratic, is a ​​linear​​ least squares problem. It might look curved when you plot it against $x$, but it is a straightforward combination of the parameters $c_1$ and $c_2$. The second model, which also produces a parabola, is a ​​nonlinear​​ least squares problem. Why? Because if you expand it, you get $c_1 x^2 - 2c_1c_2 x + c_1c_2^2$. The parameters are now tangled up—multiplied by each other and squared. They are not in the simple additive form required.

This distinction is everything. For linear problems, the error surface $S(\mathbf{c})$ is a perfect, predictable parabolic bowl. Finding its bottom is a simple matter of calculus, leading to a set of linear equations (the "normal equations") that a computer can solve in one direct step. For nonlinear problems, the error surface can be a wild landscape of winding valleys, ridges, and plateaus. Finding the lowest point is no longer a one-shot calculation; it requires an iterative search, a journey into the unknown.

Solving a nonlinear least squares problem is like being a hiker dropped into a vast, foggy mountain range, tasked with finding the absolute lowest point. The altitude at any location $(\mathbf{p})$ is given by the error function $S(\mathbf{p})$. You can only see your immediate surroundings, and from your current position, you must decide which way to step next to get lower. This is the essence of iterative algorithms like Gauss-Newton and Levenberg-Marquardt.

At each step, we stand at a point $\mathbf{p}_k$ on the error surface. We can't see the whole landscape, but we can approximate it. The core idea is to pretend, just for a moment, that our complex nonlinear model is a simple linear one in the neighborhood of our current guess. This simplification transforms the terrifying, craggy landscape into a nice, simple parabolic bowl that we do know how to solve. We find the bottom of this local, approximate bowl and jump there. That's our next guess, $\mathbf{p}_{k+1}$. Then we repeat the process: create a new approximation, find its bottom, and jump again. Hopefully, each jump takes us to a lower altitude, and we eventually converge on the true minimum of the real landscape.

The ​​Gauss-Newton algorithm​​ is the classic implementation of this strategy. To create its local approximation, it uses the ​​Jacobian matrix​​, $J$, which is the collection of all the first partial derivatives of the model's residuals with respect to each parameter. The Jacobian tells us how sensitive the model's output is to a small wiggle in each parameter.

The update step in the Gauss-Newton method is found by solving the system:

Here, $\mathbf{r}$ is the vector of current residuals (the differences $y_i - f(x_i; \mathbf{p}_k)$) and $\boldsymbol{\delta}$ is the step we should take to get to our next, better guess $\mathbf{p}_{k+1} = \mathbf{p}_k + \boldsymbol{\delta}$.

The matrix $J^T J$ is a masterpiece of approximation. It is, in fact, an approximation of the true Hessian matrix (the matrix of second derivatives) of the error function $S$. The full Hessian is $\nabla^2 S = J^T J + \sum_{i=1}^m r_i(\mathbf{x}) \nabla^2 r_i(\mathbf{x})$. The Gauss-Newton method bravely ignores the second term. This is a brilliant simplification because this second term depends on the residuals $r_i$. Near the solution, the model fits the data well, the residuals are small, and the ignored term vanishes! The approximation becomes increasingly accurate as we get closer to our goal.

However, this boldness can be the algorithm's downfall. Far from the solution, where the residuals are large, the approximation can be poor. The algorithm might confidently leap across the valley and land higher up on the other side. In some pathological cases, it can even become trapped in an oscillation, jumping back and forth between two points, never reaching the minimum that lies between them.

This is where the venerable ​​Levenberg-Marquardt (LM) algorithm​​ comes to the rescue. It is the seasoned explorer to Gauss-Newton's reckless youth. The LM algorithm introduces a "damping" parameter, $\lambda$, that adaptively controls the size and direction of each step. The update equation is modified to:

The magic lies in how $\lambda$ is adjusted:

• ​​When $\lambda$ is very small ($\lambda \to 0$):​​ The term $\lambda I$ vanishes, and the algorithm becomes the pure, bold Gauss-Newton method. This is used when previous steps have been successful, indicating we are in a well-behaved, "bowl-like" region of the landscape. We take big, confident strides.
• ​​When $\lambda$ is very large:​​ The $\lambda I$ term dominates the $J^T J$ matrix. In this limit, the update step becomes a very small step in the direction of steepest descent—the most cautious move possible, simply heading straight "downhill" from the current position. This is used when a Gauss-Newton step would have taken us uphill.

The LM algorithm is a brilliant hybrid. It "interpolates" between the fast convergence of Gauss-Newton and the guaranteed (but slow) descent of the steepest descent method. By rewarding successful steps with a smaller $\lambda$ (more boldness) and penalizing failed steps with a larger $\lambda$ (more caution), it robustly and efficiently navigates almost any landscape to find the minimum.

Faced with a nonlinear problem, a tempting shortcut is to try to transform the data to make the model linear. For an exponential model like $y = a e^{bx}$, why not just take the logarithm? Then we have $\ln(y) = \ln(a) + bx$, which looks like a simple line. Problem solved, right?

Wrong. This shortcut is a dangerous trap. When we transform our data, we also transform the experimental error. Suppose our original data has a simple additive error: $y_i = a e^{bx_i} + \varepsilon_i$, where the error $\varepsilon_i$ has a constant variance. When we take the logarithm, the new error term becomes $\ln(a e^{bx_i} + \varepsilon_i) - (\ln a + bx_i)$, which is a complicated beast. Its mean is no longer zero, and its variance now depends on $x_i$. Applying standard linear regression to this transformed data violates the core statistical assumptions and will produce biased and incorrect estimates for the parameters. The direct NLLS fit on the original, untransformed data is the statistically sound method, as it honors the original error structure of the measurement.

This brings us to the crucial idea of ​​weighting​​. The basic least-squares formulation implicitly assumes that every data point is equally trustworthy—that the variance of the measurement error is the same for all points. But what if this isn't true? In many experiments, the measurement error depends on the magnitude of the signal itself.

Consider fitting impedance data that spans many orders of magnitude. A common (but often poor) choice is "modulus weighting," which weights each point by the inverse square of its impedance magnitude, $1/|Z_i|^2$. If you have two points, one at $50 \, \Omega$ and one at $200,000 \, \Omega$, and your model misses both by the exact same absolute amount, this weighting scheme will penalize the error at the high-frequency ($50 \, \Omega$) point nearly 16 million times more than the error at the low-frequency point. The fit will be completely dominated by the low-impedance data, effectively ignoring what happens at high impedances. A proper fit requires weights that are inversely proportional to the actual error variance at each point, $w_i \propto 1/\sigma_i^2$, ensuring that each point contributes to the fit in proportion to its true reliability.

After all this searching, the NLLS algorithm gives us the "best-fit" parameters—the coordinates of the lowest point it could find. But science demands more than just a single number. It demands to know: how sure are we? If we repeated the experiment, how much would these parameters change?

Herein lies the final, beautiful piece of the puzzle. The very same matrix we used to navigate the error surface, $J^T J$, holds the key to this question. The shape of the error valley at the minimum tells us about the uncertainty of our parameters. A narrow, steep valley means the parameter is well-determined; even a small change in its value leads to a large increase in error. A wide, flat valley means the parameter is poorly determined; it can be changed quite a lot without making the fit much worse.

This "shape" is mathematically captured by the inverse of the $J^T J$ matrix. Under standard statistical assumptions, the covariance matrix of the estimated parameters is given by:

Here, $s^2$ is our estimate of the measurement error variance, calculated from the final sum of squared residuals. The diagonal elements of this covariance matrix give us the variance for each parameter, and the square root of the variance gives us the coveted ​​standard error​​—a measure of the statistical uncertainty in our best-fit values.

Thus, the journey of nonlinear least squares comes full circle. The machinery used to iteratively step toward the solution simultaneously provides the tool to assess the quality of that solution. It not only finds the bottom of the valley but also tells us how deep and narrow that valley is, providing not just an answer, but a measure of our confidence in that answer.

So, we have learned about this marvelous intellectual machine called nonlinear least squares. You give it a theory—a mathematical model of how you think some part of the world works—and you give it the facts—your experimental data. The machine churns and grinds, and it hands you back the best possible version of your theory, the specific parameters that make your model cling most tightly to reality.

This is a powerful tool, to be sure. But to treat it as a mere "curve-fitting" utility is to miss the point entirely. It is not just about drawing a line through a set of points. It is a universal lens, a way of asking precise questions of nature and getting quantitative answers. To truly appreciate its power, we must go on a journey and see it in action. We will see that this single mathematical idea provides a common thread, weaving together the vast tapestry of scientific inquiry, from the dance of distant galaxies to the secret whispers of a microchip.

Let's start by looking up at the night sky. We see a faint, fuzzy patch of light—a distant galaxy. What is its shape? Is it a perfect sphere, or is it flattened into an ellipse like a cosmic discus? To answer this, an astronomer measures the positions of many stars or glowing gas clouds within it. Each position is a single point of light. We can then propose a model, the equation of an ellipse with semi-axes $a$ and $b$. For any proposed ellipse, some data points will lie inside it, and some outside. Nonlinear least squares provides a principled way to find the one and only ellipse that best "fits" this cloud of points by minimizing the sum of squared "errors" or distances from each point to the curve. By finding the optimal $a$ and $b$, we are no longer just guessing; we have quantitatively characterized the galaxy's shape. This same principle applies to tracking the orbits of planets and asteroids, turning a collection of nightly observations into a predictable trajectory.

Now, let's come down from the heavens and into the chemist's laboratory. A flask is gently heated, and the reaction inside begins to bubble faster. We know that temperature dramatically affects the rate of chemical reactions, and Svante Arrhenius gave us a beautiful law for it: the rate $k$ depends exponentially on the inverse of the temperature $T$. This relationship, however, contains a crucial parameter: the activation energy, $E_a$, which represents the energy barrier that molecules must overcome to react. To find this fundamental value, a chemist measures the reaction rate at several different temperatures. The resulting data points don't form a straight line; they trace an exponential curve. Nonlinear least squares is the perfect tool for the job. By fitting the Arrhenius equation directly to the data, we can extract a precise estimate of $E_a$. This number is not just an abstract parameter; it is a key piece of knowledge that allows us to control chemical processes, from manufacturing pharmaceuticals to modeling the complex chemistry of our planet's atmosphere.

Let's get even more tangible and look at the world of materials and electronics. How does a battery store and release energy, or how does a metal surface corrode? These processes occur at a complex, chaotic interface. To study it, electrochemists use a technique called Electrochemical Impedance Spectroscopy (EIS). They apply a small, oscillating voltage at different frequencies and measure the oscillating current that flows in response. The relationship between this voltage and current is the impedance, a complex number that changes with frequency. To make sense of this, scientists build a simplified "map" of the interface using an equivalent circuit composed of resistors, capacitors, and other elements like the Warburg impedance, which describes diffusion. The total impedance of this model circuit is a highly nonlinear function of its component values. By using nonlinear least squares—extended to the domain of complex numbers—we can adjust the parameters of our circuit map, such as the charge-transfer resistance $R_{\text{ct}}$, until its predicted impedance spectrum perfectly matches the measured data. In doing so, we turn a bewildering set of measurements into a clear diagnosis of the hidden processes governing the battery's health or the material's decay.

The same principles that describe galaxies and batteries turn out to be indispensable for understanding the machinery of life. At the heart of biology are enzymes, the tiny protein machines that catalyze the chemical reactions in our cells. To understand an enzyme, we must know its "specifications": its maximum speed, $V_{\text{max}}$, and its affinity for its target molecule, $K_M$. The relationship between the reaction rate and the concentration of the target molecule is described by the famous Michaelis-Menten equation, another nonlinear model.

For a long time, scientists used a clever mathematical trick called the Lineweaver-Burk plot to turn this curve into a straight line, making the parameters easy to estimate with a ruler. But this trick came at a cost: it distorted the experimental error, giving undue weight to measurements at low concentrations, which are often the least reliable. Nonlinear least squares is the modern, honest approach. It fits the true Michaelis-Menten curve directly to the untransformed data, respecting the error in each measurement equally. This provides a far more accurate and reliable picture of how these vital cellular machines truly function.

From a single enzyme, we can scale up to the response of an entire cell or tissue. How does a cell respond to a hormone or a drug? Often, the response is "sigmoidal"—it's weak at low doses, rises steeply in a narrow range, and then plateaus at a maximum effect. The Hill function is a beautiful nonlinear model that captures this switch-like behavior. Fitting this model using nonlinear least squares allows pharmacologists to determine two critical parameters: the $EC_{50}$, which is the concentration needed for a half-maximal response (a measure of potency), and the Hill coefficient $n$, which describes the steepness or "cooperativity" of the response. These numbers are the bedrock of drug development and quantitative biology.

Life is not static; it is a story of populations in flux. Imagine the population of memory T cells in your body after a booster vaccine. They multiply rapidly at first, then the expansion slows as they compete for limited resources like space and signaling molecules called cytokines. This pattern of growth—fast at first, then leveling off—is ubiquitous, describing everything from yeast in a vat to fish in a pond. It is captured by the logistic growth model. The model's key parameter is the "carrying capacity" $K$, the maximum sustainable population size. By counting the cells over time and fitting the logistic model with nonlinear least squares, immunologists can estimate the carrying capacity of the "immune niche," providing a quantitative understanding of how our body regulates its defenses.

And what about the grandest scale of biology—evolution? Natural selection acts on populations over generations. Population genetics provides us with precise, but often nonlinear, equations that predict how the frequency of a beneficial gene will change over time. If we observe a population evolving in the lab—say, bacteria adapting to an antibiotic—we can track the frequency of the resistance gene over many generations. We can then use nonlinear least squares to fit the theoretical model of selection to this time-series data. The result is astonishing: we get a direct measurement of the selection coefficient $s$, the very force of evolution in action.

Our journey concludes with a surprising twist, moving from the natural world to the world of human-made information. Can a mathematical fitting procedure be used to steal a secret? In the fascinating field of cryptography, a "side-channel attack" does just that. Every time your computer performs a calculation, it consumes a tiny amount of electrical power. This power consumption is not constant; it fluctuates depending on the exact operations being performed and the data being processed.

Imagine a cryptographic algorithm that uses a secret key, $k$, to encrypt a message. It turns out that the power drawn by the processor at specific moments can be modeled as a subtle, nonlinear function of the secret key and the (known) input data. An attacker can carefully measure these power fluctuations—the physical "leakage" of the computation. They now have a set of data points. Using their knowledge of the hardware, they write down a nonlinear model that predicts the power leakage given a key. Then, they turn to nonlinear least squares. The algorithm finds the value of the key $k$ that makes the model's predictions best match the observed power traces. The model, born from physics and statistics, becomes a skeleton key to unlock a digital secret.

As we have seen, the applications of nonlinear least squares are as diverse as science itself. Yet, a unifying theme emerges. The method is more than just a tool for finding parameters. It is the engine of the scientific method, the quantitative link between theory and experiment.

We start with an idea, a story about how the world works, and we formalize it as a mathematical model. We then go out and collect the facts. Nonlinear least squares provides the stage for the confrontation between our idea and the facts. It tells us, in no uncertain terms, what the best version of our story is.

But perhaps the most profound part of this process is what happens when the fit is not good. When our best model still leaves a large, systematic pattern in the residuals—the leftover differences between the model and the data—we have not failed. We have discovered something wonderful. Those residuals are a message from nature, a clue that our story is incomplete. They are a treasure map pointing toward a deeper, more beautiful, and more accurate law. This iterative process of modeling, fitting, and analyzing the residuals is the very essence of scientific discovery—a perpetual and fruitful dialogue with the universe.