Barycentric Lagrange Interpolation

In my studies of the Remez algorithm, I learned about the barycentric Lagrange interpolation formula.

The context is finding a polynomial of degree at most $n$ that passes through $n+1$ points $(x_0, y_0), \dots, (x_n, y_n)$. The classical Lagrange interpolation formula is what you’d write down if you “just did it.”

$$f(x) = \sum_{i=0}^n y_i \cdot \prod_{j \neq i}\frac{x - x_j}{x_i - x_j}$$

I wrote a 2014 article deriving this more gently, and implementing it in Haskell for secret sharing. An even gentler, Python-based version was the “pracitcal application” of chapter 2 of PIM.

I put practical application in quotes because the classical formula is numerically unstable. I acutally didn’t know this at the time I wrote PIM, but for the application of secret sharing it doesn’t really matter because you do the math over a finite field and numerical accuracy is irrelevant.

But now that I’m in my polynomial approximation era, I actually do care about numerical precision and stability. And for that purpose, the barycentric formula is much better.

The main observation here is that the terms of the classical formula can be thought of as divisors of a single polynomial $\ell(x) = \prod_{i=0}^n (x - x_i)$.

Then each term is $\ell_j(x) = \ell(x) / (x-x_j)$, multiplied by the term-specific weight $w_j = \prod_{i \neq j} \frac{1}{(x_i - x_j)}$. Then you can factor out $\ell(x)$, precompute the $w_j$, and get

$$ f(x) = \ell(x) \sum_{i=0}^n \frac{y_i w_i}{x - x_i} $$

This is more efficient, but you can also avoid computing $\ell(x)$ entirely by using the mathematician’s favorite trick of dividing by a cleverly-disguised 1. In this case, the 1 is represented by its interpolation formula (including $\ell(x)$), which is achieved by just setting all $y_i$ to 1.

When you divide, the $\ell(x)$ cancels out and you get

$$ f(x) = \left [ \sum_{i=0}^n \frac{y_i w_i}{x - x_i} \right ] / \left [ \sum_{i=0}^n \frac{w_i}{x - x_i} \right ] $$

This is the barycentric formula. It is numerically stable for well-chosen point sets $x_i$ (see section 7 of the paper of Trefethen linked below). I have to admit, I have no idea what this has to do with barycentric coordinates as I originally learned about them. I guess the lagrange terms $\ell_j(x)$ are supposed to be the vertices of a simplex, and the $w_j$ the barycentric weights, but the calling this a simplex seems weird to me.

The other nice thing about this formula is that the weights are trivial to compute for certain point sets. For uniform spaced points on $[-1, 1]$, $w_j = (-1)^j\binom{n}{j}$, and for Chebyshev points on the same domain it’s $w_j = (-1)^j$ with values halved at $j=0, j=n$.

I wrote some Python code to implement this here, along with one extra tweak to improve numerical stability, which boils down to multiplying the numerator and denominator of each term in the sum by 2 to avoid underflow. The test includes an interpolation of degree 100,000 that runs in about a second in vanilla Python.

Popularizing barycentric Lagrange interpolation seems to be something of a passion project for Lloyd N. Trefethen. Trefethen, who until 2023 was a member of the numerical analysis group at Oxford, is an author on basically every paper I can find on this topic, in addition to the applications of this to Remez-like algorithms. In this paper, section 10 gives a history, of the curious obscurity of this formula. Among other interesting facts, Trefethen states that Richard Hamming included it in the first edition of Numerical Methods for Scientists and Engineers, and then dropped it in the second edition with no explanation.