Wilkinson's Polynomial

Polynomials are a fundamental tool in science and engineering, used to model everything from the arc of a projectile to the stability of a control system. A core task associated with them is root-finding, a process that seems direct and reliable. We intuitively expect that a small adjustment to a polynomial's coefficients will result in only a small shift in its roots. However, this intuition can be dangerously flawed, leading to catastrophic computational errors. This article addresses this critical knowledge gap by exploring the phenomenon of numerical instability through its most famous and startling example: Wilkinson's polynomial.

This article will guide you through this profound computational lesson. In the "Principles and Mechanisms" chapter, we will dissect the concept of ill-conditioning, uncover the mathematical reason for this sensitivity, and witness the dramatic failure of intuition with Wilkinson's classic experiment. Following that, the "Applications and Interdisciplinary Connections" chapter will demonstrate that this is not just a mathematical curiosity, but a crucial consideration with far-reaching consequences in control engineering, digital signal processing, and other applied sciences, revealing how engineers have learned to build robust systems in the face of this numerical ghost.

In our journey to understand the world, we often lean on familiar friends. Polynomials are one such friend. From the simple parabolic arc of a thrown ball to the complex characteristic equations that govern the stability of a satellite, polynomials are everywhere. We learn in school how to find their roots—those special values of $x$ that make the polynomial equal to zero. The process seems straightforward: you are given the coefficients, and you are asked to find the roots. It feels like a sturdy, reliable bridge connecting two well-defined places.

But what if the bridge is built on shaky ground? What if the numbers we use for the coefficients—numbers that come from physical measurements or are stored in the finite memory of a computer—are not perfectly precise? Our intuition tells us that a tiny nudge to a coefficient should only result in a tiny nudge to the roots. A small cause should lead to a small effect. This, however, is a beautiful and comfortable intuition that is, in certain circumstances, catastrophically wrong. The story of Wilkinson's polynomial is the story of discovering just how wrong it can be.

Let's start with a simple, tangible example. Imagine you're a control systems engineer designing a new drone. Its stability might depend on the roots of a characteristic polynomial, say, $P(x) = x^3 - 15x^2 + 74x - 120 = 0$. You've designed it perfectly, and you know the roots are exactly $4$, $5$, and $6$. These roots tell you that your drone is stable.

Now, a problem arises. A sensor in the drone has a tiny, persistent bias. This bias perturbs the coefficient of the $x^2$ term from $-15$ to $-15 - \epsilon$, where $\epsilon$ might be a minuscule number like $4.0 \times 10^{-5}$. Will the drone remain stable? Will the roots move much? We can estimate the change in the root that was originally at $x=5$. A first-order analysis reveals that the change in the root, $\delta r$, is approximately $\delta r \approx -25 \epsilon$.

Let's pause and appreciate what this means. The perturbation to the root is 25 times larger than the perturbation to the coefficient. This amplification factor is a simple example of a ​​condition number​​. It measures how sensitive a problem's output is to changes in its input. A problem with a small condition number is called ​​well-conditioned​​; our intuition about "small cause, small effect" holds true. A problem with a large condition number is called ​​ill-conditioned​​, and it's a warning sign for any scientist or engineer. In another scenario involving a satellite's control system, a polynomial of degree 3 was found to have a root sensitivity of 121. A small error in an amplifier's gain could be magnified over 100 times in the resulting system pole, a potentially dangerous situation.

Why does this amplification happen? The secret lies in a wonderfully simple formula derived from basic calculus. If we have a polynomial $P(x)$ and we introduce a small change to it, let's call it $\Delta P(x)$, then the first-order change in a root $r$ is given by:

$\delta r \approx -\frac{\Delta P(r)}{P'(r)}$

The numerator, $\Delta P(r)$, is simply the value of the perturbation evaluated at the root. The magic is in the denominator: $P'(r)$, the derivative of the polynomial at the root. Geometrically, $P'(r)$ is the slope of the polynomial's graph as it crosses the x-axis at $x=r$.

Imagine a river crossing a plain. If the riverbank is steep, moving the entire landscape up or down by a small amount will shift the water's edge horizontally by only a little. But if the river is flowing through a very flat, wide floodplain, the same small vertical shift can move the water's edge by a huge distance.

The same is true for polynomial roots. If the polynomial crosses the x-axis steeply (i.e., $|P'(r)|$ is large), the root is stable. A vertical shift $\Delta P$ won't move it far. But if the polynomial crosses the x-axis at a very shallow angle (i.e., $|P'(r)|$ is close to zero), the root is unstable. The tiniest vertical nudge can send the root flying along the x-axis. This is the mechanism of ill-conditioning. The problem becomes extremely sensitive when the derivative at the root is small.

This brings us to the main character of our story. In the 1960s, the brilliant numerical analyst James H. Wilkinson was studying the errors in computer calculations. He considered a seemingly innocuous polynomial:

$W_{20}(x) = \prod_{k=1}^{20} (x-k) = (x-1)(x-2)\cdots(x-20)$

The roots are, by construction, the integers $1, 2, 3, \ldots, 20$. They are all distinct and nicely separated by a distance of 1. It seems like the most well-behaved polynomial one could imagine. What could possibly go wrong?

Wilkinson decided to expand this polynomial into its monomial form, $W_{20}(x) = a_{20}x^{20} + a_{19}x^{19} + \dots + a_0$. The coefficients turn out to be astronomically large. For instance, the coefficient of $x^{19}$ is $a_{19} = -(1+2+\dots+20) = -210$, while the constant term $a_0$ is $20! \approx 2.43 \times 10^{18}$.

Then, he performed a thought experiment, which can be replicated today with a simple program. He took the expanded form and made a single, almost imperceptible perturbation. He changed the coefficient of $x^{19}$ from $-210$ to $-210 - 2^{-23}$. This is a tiny relative change, about one part in $10^{10}$. He then asked his computer to find the roots of this slightly modified polynomial.

The result was not a slight nudge. It was an earthquake.

• The roots at 1, 2, ..., 8 remained more or less in place.
• But the roots from 9 onwards were dramatically shifted. The roots at 10 and 11, for example, fused together and split off the real axis, becoming a pair of complex numbers. The same happened to the roots at 12 and 13, 14 and 15, and so on.
• The root that was supposed to be at 20 moved all the way to about 20.8.

A change so small it was literally invisible in the input data had completely changed the qualitative nature of the solution. This wasn't just a quantitative error; it was a numerical catastrophe.

The theory we just discussed explains this perfectly. The sensitivity of a root $j$ to a perturbation in the coefficient of $x^{n-1}$ is enormous for this polynomial. For example, the sensitivity of the root at $x=15$ is proportional to $\frac{15^{19}}{P'(15)}$. The denominator, $P'(15) = (14!) \times (-1)(-2)(-3)(-4)(-5) = -(14!)(5!)$, is a huge number. But the numerator, $15^{19}$, is incomprehensibly larger. The ratio, which acts as the condition number for this root, is gigantic. Even for a smaller degree-10 version of this polynomial, a perturbation of $\epsilon=2.0 \times 10^{-7}$ to the $x^9$ coefficient is enough to shift the root at $x=8$ by a noticeable amount. And for a degree-5 version, the relative error in a root can be amplified by a factor of nearly 80. The effect gets exponentially worse as the degree increases.

So, is the problem of finding the roots of $W_{20}(x)$ inherently cursed? Is this polynomial a mathematical monster? The final twist in the tale is perhaps the most profound. The problem is not with the polynomial itself, but with how we chose to represent it.

Think about it this way. The product form, $(x-1)(x-2)\cdots(x-20)$, is a perfectly stable and clear description. If you want to find its roots, you can just read them off! If you want to evaluate it at, say, $x=30$, you just compute $(29) \times (28) \times \cdots \times (10)$. This is a sequence of simple, stable multiplications. A computer does this with very high accuracy.

The monstrosity arises when we insist on expanding it into the monomial basis: $a_{20}x^{20} + a_{19}x^{19} + \dots + a_0$. When we evaluate this form at $x=30$, the term $a_{20}x^{20}$ is a massive positive number, and $a_{19}x^{19}$ is a massive negative number of nearly the same magnitude. We are calculating these titans and then subtracting them to get a comparatively tiny result. This is a classic numerical sin known as ​​catastrophic cancellation​​. Any tiny floating-point error in computing the large terms gets magnified enormously when their difference is taken.

The monomial basis is like a house of cards. The information about the root locations is encoded in the coefficients in an incredibly fragile, distributed, and non-intuitive way. The slightest touch to one card (one coefficient) can cause the entire structure to collapse. The product form, on the other hand, is like a solid brick wall. The information about each root is local and robust.

The lesson of Wilkinson's polynomial is therefore deep and far-reaching. The numerical stability of a problem is not an intrinsic property of the abstract mathematical question alone; it is fundamentally tied to the ​​representation​​ we choose. The coefficients of the monomial expansion are an exquisitely ill-conditioned way to represent a polynomial whose roots are clustered. This discovery was a wake-up call, teaching us that the art of scientific computing is not just about finding algorithms, but about finding the right way to represent our problems so that the answers are not lost in a sea of numerical noise. It's a powerful reminder that in the dance between mathematics and machine, we must always be mindful of our steps.

In the last chapter, we took a journey into a seemingly esoteric corner of mathematics. We met Wilkinson’s polynomial, a strange and wonderful creature whose roots, the integers from 1 to 20, are exquisitely sensitive to the tiniest whisper of a change in its coefficients. A perturbation smaller than the dust on a butterfly's wing can send the roots scattering wildly across the complex plane.

It's a fascinating story, but one might be tempted to ask, "So what?" Is this just a mathematical curiosity, a pathology confined to the pristine world of pure abstractions? Or does this ghost haunt the real-world calculations that underpin our modern technological world?

The answer, perhaps surprisingly, is that it haunts them profoundly. What Wilkinson uncovered was not an isolated monster but a fundamental truth about the nature of computation. The journey we are about to take is into the applied world—into the design of jet engines, the filtering of digital signals, and the analysis of stressed materials—to see where the echoes of Wilkinson’s ghost appear and how engineers and scientists have learned to confront them. This is the story of how an abstract mathematical insight forces us to be cleverer in how we build our world.

At its heart, the problem Wilkinson posed was about finding the roots of a polynomial. The classical approach, the one we all learn in school, involves working with the polynomial's coefficients. But as we've seen, this path is fraught with peril. So, is there a better way?

There is, and it comes from a beautiful pivot in perspective: translating the problem from the language of algebra into the language of linear algebra. For any monic polynomial, like $p(x) = x^{n} + a_{n-1} x^{n-1} + \dots + a_{0}$, one can construct a special matrix called its ​​companion matrix​​. A common form is:

The magic of this matrix is that its eigenvalues are precisely the roots of the polynomial $p(x)$. Suddenly, our root-finding problem has become an eigenvalue problem! This is a tremendous conceptual leap. But does it solve our numerical troubles?

Not by itself. The coefficients are still there, sitting in that last column. If we take Wilkinson's polynomial for $n=20$, compute its coefficients, and build the companion matrix, what happens if we introduce a tiny perturbation to one of those coefficients before computing the eigenvalues? As one might expect, the result is the same numerical catastrophe. The eigenvalues—our polynomial's roots—scatter just as before. This tells us something crucial: the problem is not in the specific root-finding algorithm, but in the representation itself. The set of coefficients is simply a treacherous, unstable way to describe a set of roots.

The real breakthrough comes from how we solve the eigenvalue problem. Powerful algorithms, like the ​​QR algorithm​​, were developed to find the eigenvalues of any matrix, and they do so in a remarkably stable way. The QR algorithm works through a series of clever matrix transformations that don't even require calculating the characteristic polynomial. It nibbles away at the matrix, iteratively revealing the eigenvalues with astonishing reliability.

So here is the first great lesson learned from Wilkinson's ghost. The modern, robust way to find the roots of a high-degree polynomial is a two-step dance: first, form the companion matrix from the coefficients, and second, apply a stable eigenvalue solver like the QR algorithm directly to that matrix. We have sidestepped the ill-conditioned problem of finding roots from coefficients by transforming it into the much more well-behaved problem of finding eigenvalues from a matrix. The ghost has been, for the moment, contained.

Let's leave the abstract world of polynomials and venture into the domain of a control engineer, who might be designing the flight control system for a supersonic jet or the process controller for a chemical plant. These systems are often described by "transfer functions," which are essentially ratios of polynomials, say $G(s) = b(s)/a(s)$. The roots of the denominator polynomial $a(s)$ are called the system's ​​poles​​, and they are of paramount importance: they determine whether the system is stable or will spiral out of control.

For analysis, engineers love to use standardized representations called "canonical forms." One of the most famous is the ​​controllable canonical form​​, which represents the system using a set of state-space matrices. And what does the main state matrix in this form look like? You guessed it—it's a companion matrix built from the coefficients of the denominator polynomial $a(s)$.

Wilkinson's ghost has just reappeared, this time in a flight controller! If a system has poles that are close together—a very common situation in, say, flexible structures like aircraft wings—then its representation in this "simple" canonical form becomes numerically poisonous. The coefficients of $a(s)$ become extremely sensitive, and any floating-point representation of the system in a computer is liable to be dangerously inaccurate.

The situation becomes even more acute in control design. A fundamental task is "pole placement," where an engineer designs a feedback law to move the system's poles to more desirable, stable locations. Some classical methods for doing this, like ​​Ackermann's formula​​, or techniques that first transform the system to the companion form, are built directly upon this treacherous polynomial-coefficient foundation. For systems that are high-order or "nearly uncontrollable"—systems that are inherently difficult to influence in certain ways—these methods can fail spectacularly. Trying to force a nearly uncontrollable system into the theoretically elegant companion form requires a numerically unstable transformation that acts like a massive amplifier for any roundoff error, leading to a completely useless result.

The lesson for generations of engineers has been a hard-won one: beware of theoretical elegance that ignores numerical reality. The companion form, so neat on the blackboard, can be a death trap in a computer. Modern, robust control design has largely abandoned these coefficient-based methods in favor of state-space techniques that rely on stable matrix operations like the QR algorithm or the Singular Value Decomposition (SVD), keeping Wilkinson's ghost safely at bay.

Our journey now takes us to the world of digital signal processing (DSP), the technology behind your phone calls, streaming music, and medical imaging. A cornerstone of DSP is the ​​digital filter​​, an algorithm designed to modify a signal, for instance, to remove noise. These filters are often described by their poles and zeros in the complex plane.

Suppose you have a high-order filter with many poles and zeros. How would you compute its frequency response, which tells you how it affects different tones? One seemingly obvious method is to take all the poles and zeros, expand them to get the coefficients of the numerator and denominator polynomials, and then evaluate these polynomials at each frequency.

An alternative is a "geometric" method, which works directly with the factored pole-zero form. It calculates the frequency response by summing up the contributions (in terms of distance and angle) from each pole and zero to the point of interest on the unit circle. A carefully designed numerical experiment pits these two methods against each other. The result is a dramatic confirmation of our theme: for high-order filters, especially those with poles and zeros clustered together, the polynomial coefficient method suffers a catastrophic loss of precision. The geometric method, by avoiding the toxic intermediate step of forming coefficients, remains perfectly accurate.

This principle extends to more advanced problems. Consider ​​spectral factorization​​, a task where one has a given power spectrum (a description of the signal's power at each frequency) and wants to find the filter that would produce such a signal. Once again, two paths present themselves. One involves representing the spectrum as a polynomial and finding its roots to construct the filter. The other involves a sophisticated state-space approach that solves a matrix equation called the ​​Discrete-time Algebraic Riccati Equation (DARE)​​. For high-order systems, the polynomial route is known to be ill-conditioned and unreliable, especially when trying to correctly identify which roots correspond to a stable filter. In contrast, the Riccati-based method, which relies on robust matrix decompositions, is the professional's tool of choice for its stability and reliability.

The message is clear and consistent. Whether we are analyzing a control system or designing a digital filter, representing the system by its polynomial coefficients is a numerically fragile choice.

You might think this is a story about electrical engineering and computers. But the principle is universal. Let's make one last stop, in the field of solid mechanics. When a steel beam in a bridge or a rock deep in the Earth's crust is under load, the internal state of force is described by a ​​Cauchy stress tensor​​—a small, $3 \times 3$ symmetric matrix. To understand if the material is close to fracturing, engineers need to find the "principal stresses," which are simply the eigenvalues of this stress matrix.

How should they be computed? One could form the characteristic cubic polynomial, $\det(\boldsymbol{\sigma} - \lambda \boldsymbol{I}) = 0$, and solve it using the cubic formula. Or, one could treat it as a matrix eigenvalue problem and use a standard numerical routine, like the QR algorithm.

Even in this seemingly simple $3 \times 3$ case, the lessons of numerical analysis hold true. The direct matrix approach via the QR algorithm is recognized as being more robust, more reliable, and more informative, as it provides the principal directions (the eigenvectors) for free. The polynomial route can suffer from loss of precision if two of the principal stresses are nearly equal (i.e., clustered eigenvalues), and it is generally seen as a less dependable approach in professional-grade software.

From finding the modes of a vibrating airplane wing to solving for energy levels in quantum mechanics, any time a problem in science or engineering boils down to finding the eigenvalues of a matrix, the same choice presents itself. And time and again, experience has shown that taking a detour through the characteristic polynomial is a path best avoided.

Wilkinson's polynomial, therefore, is far more than a mathematical curiosity. It is a profound cautionary tale, a foundational parable of the digital age that has shaped the very landscape of scientific and engineering computation. Its lesson is not that math is wrong, but that we must be exquisitely careful in how we translate our mathematical ideas into finite, physical computation.

The central theme is the critical importance of representation. The list of a polynomial's coefficients, so simple and compact on paper, is often a numerically poisonous way to store information. In its place, scientists and engineers have learned to favor more robust representations: state-space matrices, factored pole-zero forms, and direct matrix formulations. These, when paired with stable algorithms like the QR and SVD, are the tools that allow us to build reliable models of our complex world.

This journey, from a single polynomial to the breadth of modern engineering, reveals a beautiful unity in scientific computing. The same fundamental principle of numerical stability, the same ghost in the machine, echoes through discipline after discipline. It is a powerful reminder that to build reliable things in the physical world, we must first deeply understand the subtle and sometimes surprising nature of the numbers we use to describe it.