Open main menu
Barycentric Interpolation

Barycentric Interpolation

In this research article, I review the barycentric interpolation method as a more efficient form for working with Lagrange bases especially when paired with nth roots of unity.

Written by

Seun LanlegeSeun Lanlege

Published on

Do not index
Do not index
Special thanks to Alistair Stewart for the helpful discussions
 
In this research article, I review the barycentric interpolation method
[1]^{[1]}
 as a more efficient form for working with Lagrange bases especially when paired with nth roots of unity.
 
First we define
ℓ(x)\ell(x)
 as the product of terms:
 
ℓ(x)=∏i=0n(x−xi)\ell(x) = \prod^n_{i=0} (x - x_i)
 
Next we define the barycentric weight
wiw_i
 as:
 
wi=∏j=0,j≠in(xi−xj)−1w_i = \prod^{n}_{j=0, j \ne i} (x_i -x_j) ^{-1}
 
Such that the Lagrange base
Li(x)\mathcal{L}_i(x)
 can be written as:
 
Li(x)=ℓ(x)wix−xi\begin{equation} \mathcal{L}_i(x) = \ell(x) \frac{w_i}{x-x_i} \tag{1} \end{equation}
This works because:
 
Li(x)=δij=∏j=0,i≠jn(x−xjxi−xj)=∏i=0n(x−xi)∏j=0,j≠in(xi−xj)(x−xi) \mathcal{L}_i(x) = \delta_{ij} = \prod\limits_{j = 0, { i \ne j}}^{n} (\frac{x-x_j}{x_i - x_j}) = \frac{\prod^n_{i=0} (x - x_i)}{\prod^{n}_{j=0, j \ne i} (x_i -x_j)(x-x_i)}
 
(The
(x−xi)(x - x_i)
 term in the numerator and denominator both cancel out)
 
Where
δij\delta_{ij}
 is the Kronecker delta function. And so our new Lagrange polynomial becomes:
 
ϕ(x)=∑i=0nℓ(x)⋅wix−xi⋅vi\phi(x) = \sum^{n}_{i = 0} \frac{\ell(x) \cdot w_i}{x-x_i} \cdot v_i
 
Since
ℓ(x)\ell(x)
 is a common term we can factor it out such that we have:
 
ϕ(x)=ℓ(x)∑i=0nwix−xi⋅vi\begin{equation} \phi(x) = \ell(x) \sum^{n}_{i = 0} \frac{ w_i}{x-x_i} \cdot v_i \tag{2} \end{equation}
 
This form allows us to precompute
ℓ(x)\ell(x)
 and only have to compute the individual barycentric weights for each
ii
-th term. This gives us a much better complexity than the naive Lagrange interpolation polynomial.
 
Equation
(2)(2)
 is also known as the “First form of the barycentric interpolation formula”.
 
Polynomial interpolation breathing
Polynomial interpolation breathing
 
But we’re not done, we need to re-write the barycentric weight
wiw_i
 as a function, rather than an expression. Notice that we can rewrite
wiw_i
 as:
 
wi=1∏j=0,j≠in(xi−xj)=1ℓ′(xi)\begin{align*} w_i &= \frac{1}{\prod^{n}_{j=0, j \ne i} (x_i -x_j)} \\ &= \frac{1}{\ell^{\prime}(x_i)} \end{align*}
 
Where
ℓ′(x)\ell^{\prime}(x)
 is the derivative of
ℓ(x)\ell(x)
. In order to understand how we arrived at this derivative expression, Let’s examine an example for
n=4n = 4
 
ℓ(x)=∏i=04(x−xi)\ell(x) = \prod^{4}_{i = 0}(x - x_i)
 
In order to compute the derivative of this function we leverage the product rule:
 
(u⋅v)′=u′⋅v+u⋅v′(u \cdot v)^{\prime} = u^{\prime} \cdot v + u \cdot v^{\prime} \\
 
Thus the derivative becomes:
 
ℓ′(x)=(x−x2)(x−x3)(x−x4)+(x−x1)(x−x3)(x−x4)+(x−x1)(x−x2)(x−x4)+(x−x1)(x−x2)(x−x3) \begin{split} \ell^{\prime}(x) &= (x - x_2)(x - x_3)(x - x_4)\\&+(x - x_1)(x - x_3)(x - x_4)\\&+(x - x_1)(x - x_2)(x - x_4) \\&+(x - x_1)(x - x_2)(x - x_3)\\ \end{split}
 
Evaluating the derivative
ℓ′(x)\ell^{\prime}(x)
 at
x1x_1
 gives us:
 
ℓ′(x1)=(x1−x2)(x1−x3)(x1−x4)\ell^{\prime}(x_1) = (x_1 - x_2)(x_1 - x_3)(x_1 - x_4)
 
From observing how the derivative behaves at other coordinates
x1,x2…xnx_1, x_2 \dots x_n
, we can generalise the derivative as:
 
ℓ′(x)=∑i=0n(∏j=0,i≠jn(x−xj)) \ell^{\prime}(x) = \sum^{n}_{i = 0} (\prod^{n}_{j = 0, i \ne j} (x -x_j)) \\
 
So we see now how our barycentric weight
wiw_i
 can be re-written as a function of
xix_i
:
 
wi=∏j=0,i≠jn(xi−xj)−1=ℓ′(xi)−1\begin{split} w_i &= \prod^{n}_{j = 0, i \ne j} (x_i -x_j)^{-1} \\ &= \ell^{\prime}(x_i)^{-1} \end{split}
 

nth Roots of Unity

Another optimisation we can make here is to reduce the computation time for the terms
ℓ(x)\ell(x)
 and
ℓ′(xi)\ell^{\prime}(x_i)
 to be logarithmic, by leveraging the
nn
th roots of unity. The
nn
th Roots of unity are simply the solutions to the equation:
 
Xn=1 X^n = 1
 
For example, the solutions to the equation at
n=2n = 2
, are
,x=1,−1, x = 1, -1
. This is because,
12=(−1)21^2 = (-1)^2
. At
n=4n=4
, we run into a problem, the first two solutions are
11
 &
−1-1
 same as before, but other what two numbers raised to
44
 will give us
11
? Well logically that would be
1\sqrt{1}
 and
−1\sqrt{-1}
 but there are no solutions for this expression in the real numbers. Instead, they are referred to as
−i-i
 and
ii
, imaginary units which lie in the complex plane. Usually we’d be working in finite fields or cyclic groups where there are indeed solutions to these expressions. The fundamental theorem of algebra states that any
nn
th degree polynomial has
nn
 roots. What this means is that we can factor the equation above into the form:
 
ℓ(x)=Xn−1=∏i=0n(x−ωi)\begin{equation} \ell(x) = X^n - 1 = \prod^{n}_{i = 0}(x - \omega^i) \tag{3} \\ \end{equation}
 
Where
ω\omega
 is some primitive root of unity such that
ωn=1\omega^n = 1
. This primitive root of unity has an order
nn
, that is to say that
ω\omega
 generates a multiplicative group of order
nn
. It's worth noting that this group is an Abelian group. Equation
(3)(3)
 is also referred to as the Zero or Vanishing Polynomial. We can show this is true by recalling basic group theory, for any positive integer
k<nk < n
:
 
ωkn=(ωk)n=(ωn)k=1k=1\omega^{k^n} = (\omega^k)^n = (\omega^n)^k = 1^k = 1
 
So this means equation
(3)(3)
 will “vanish” on all powers of
ωi ∀i∈[n]\omega^i \space \forall i\in [n]
. Finally, the barycentric weight, is just the derivative of equation
(3)(3)
 :
 
ℓ′(x)=nx(n−1)\ell^{\prime}(x) = nx^{(n-1)}
 
So our new Lagrange base becomes:
 
Li(x)=xn−1nωi(n−1)(x−ωi)=ωi(xn−1)n(x−ωi)\begin{split} \mathcal{L}_i(x) &= \frac{x^n - 1}{n\omega^{i(n-1)}(x - \omega^i)} \\ &= \frac{\omega^i(x^n - 1)}{n(x - \omega^i)} \end{split}
 
 

References

[1]^{[1]}
 Barycentric Lagrange Interpolation

Written by

Seun Lanlege
Seun Lanlege

Mad scientist.