Lagrange Interpolation Primer

You may remember, when in High School, being asked to find the equation of a line that goes through two points. The process was straightforward: you would calculate the slope $a$ and the intercept $b$ of the line $f(x) = a x + b$.

This is a simple example of polynomial interpolation, where we are trying to find a polynomial that goes through a set of points. While the process is simple for a line, it becomes slightly more complex for higher degree polynomials.

Polynomials

One of the most, if not the most important property of polynomials, is that each one of them can uniquely be described by $n+1$ points for a polynomial of degree $n$, noted $f(x) in PP_n$.

The most straightforward way of seeing this is by setting an equation for each point $A_i (x_i, y_i)$ we have:

We can see that we have $n+1$ equations for $n+1$ factors of $x$, noted $u_i | 0 <= i <= n$, which could also be represented with the following Vandermonde matrix equation:

Written as $V u = y$, we can solve this by inverting $V$:

This is doable by applying the Jordan-Gauss algorithm to the augmented matrix $V | I$, where $I$ is the identity matrix (diagonal filled with $1$s), until we only have pivots equal to $1$:

If you prefer to see this in a graphical way, let's consider the following example:

We have two points $A_0 (1, 2)$ and $A_1 (-2, -3)$. We can represent them as follows:

#import "helpers.typ": *

#school(
  ..pin(1, 2, [$A_0$], align: top + right),
  ..pin(-2, -3, [$A_1$], align: top + right),
)

It is visually obvious that we can only draw a single line that goes through both points.

#import "helpers.typ": *

#let m = (-3 - 2) / (-2 - 1)
#let f(x) = m * x + 2 - m

#school(
  curve(f, -5, 5),
  ..pin(1, 2, [$A_0$], align: bottom + right),
  ..pin(-2, -3, [$A_1$], align: bottom + right),
)

A line is essentially just $f(x) = a x + b$, which is of degree $1$. We can guess and deduce we need at least $n+1$ points to describe a polynomial function $f(x) in PP_n$.

The same goes if we add a third point $A_2 (0, 4)$:

#import "helpers.typ": *

#let f(x) = -(11 / 6) * calc.pow(x, 2) - x / 6 + 4

#school(
  curve(f, -5, 5),
  ..pin(1, 2, [$A_0$]),
  ..pin(-2, -3, [$A_1$], align: top + left),
  ..pin(0, 4, [$A_2$], align: top + right),
)

If we did not add this third point and still tried to find a polynomial $P_2 (x)$, we would have an infinite number of solutions:

#import "helpers.typ": *

#let f(x, b) = ((3 * b - 5) / 3) * calc.pow(x, 2) + b * x + ((11 - 6 * b) / 3)

#school(
  curve(x => f(x, 1), -5, 5, color: theme.series.at(0)),
  curve(x => f(x, 2), -5, 5, color: theme.series.at(1)),
  curve(x => f(x, 3), -5, 5, color: theme.series.at(2)),
  ..pin(1, 2, [$A_0$], align: top + left),
  ..pin(-2, -3, [$A_1$], align: bottom + left, gap: 7pt),
)

But now the question is: how do we actually find the polynomial that goes through all the points? This is where Lagrange interpolation comes into play.

Lagrange Interpolation

Let's take $A_0 (0, 1)$, $A_1 (1, 3)$, $A_2 (2, 2)$ and $A_3 (3, 4)$ to demonstrate this method.

#import "helpers.typ": *

#school(
  xlim: (-0.5, 5), ypos: left, ylim: (0, 5),
  ..pin(0, 1, [$A_0$], align: left),
  ..pin(1, 3, [$A_1$]),
  ..pin(2, 2, [$A_2$], align: bottom),
  ..pin(3, 4, [$A_3$], align: bottom + right),
)

The main principle behind this is to split the wanted function into multiple sub-functions $l_i (x)$, that each contribute to one given point, also called "node".

We want to construct $l_i (x)$ so that:

Making a function equal to $0$ at a certain point $x_j$ is as simple as multiplying it by $(x - x_j)$. Based on this property, we can already "guess" $l_1^* (x)$:

#import "helpers.typ": *

#let l1(x) = (x - 0) * (x - 2) * (x - 3)

#school(
  xlim: (-0.5, 4), ypos: left, ylim: (-1, 5), legend: top + left,
  curve(l1, 0, 5, label: $l_1^* (x)$),
  ..pin(0, 1, [$A_0$], align: left),
  ..pin(1, 3, [$A_1$]),
  ..pin(2, 2, [$A_2$], align: bottom),
  ..pin(3, 4, [$A_3$], align: bottom + right),
)

$(x-x_j)$ add a root ?

For a function $f(x) = (x - x_j) dot g(x) | g(x) <- RR[x]$, when $x = x_j$:

That's a good start! But we have a problem: $l_1^* (x)$ is not equal to $1$ at $x = 1$. We can fix this by dividing $l_1^* (x)$ by $l_1^* (1)$:

And we effectively have $l_1 (1) = 1$. We can now repeat this process for $l_0 (x)$, $l_2 (x)$ and $l_3 (x)$:

#import "helpers.typ": *

#let l0(x) = (x - 1) * (x - 2) * (x - 3) / -6
#let l1(x) = (x - 0) * (x - 2) * (x - 3) / 2
#let l2(x) = (x - 0) * (x - 1) * (x - 3) / -2
#let l3(x) = (x - 0) * (x - 1) * (x - 2) / 6

#school(
  xlim: (-0.5, 4), ypos: left, ylim: (-1, 5), legend: top + left,
  curve(l0, 0, 5, color: theme.series.at(0), label: $l_0 (x)$),
  curve(l1, 0, 5, color: theme.series.at(1), label: $l_1 (x)$),
  curve(l2, 0, 5, color: theme.series.at(2), label: $l_2 (x)$),
  curve(l3, 0, 5, color: theme.series.at(3), label: $l_3 (x)$),
  ..pin(0, 1, [$A_0$], color: theme.series.at(0), align: left),
  ..pin(1, 3, [$A_1$], color: theme.series.at(1), align: top + right),
  ..pin(2, 2, [$A_2$], color: theme.series.at(2), align: top),
  ..pin(3, 4, [$A_3$], color: theme.series.at(3), align: bottom + right),
)

The polynomial $f(x) in PP_3$ is computed by summing all the $l_i (x)$ together and multiplying them by the corresponding $y_i$:

#import "helpers.typ": *

#let l0(x) = (x - 1) * (x - 2) * (x - 3) / -6 * 1
#let l1(x) = (x - 0) * (x - 2) * (x - 3) / 2 * 3
#let l2(x) = (x - 0) * (x - 1) * (x - 3) / -2 * 2
#let l3(x) = (x - 0) * (x - 1) * (x - 2) / 6 * 4

#school(
  xlim: (-0.5, 4), ypos: left, ylim: (-1, 5), legend: top + left,
  curve(l0, 0, 5, color: theme.series.at(0), label: $l_0 (x) dot y_0$),
  curve(l1, 0, 5, color: theme.series.at(1), label: $l_1 (x) dot y_1$),
  curve(l2, 0, 5, color: theme.series.at(2), label: $l_2 (x) dot y_2$),
  curve(l3, 0, 5, color: theme.series.at(3), label: $l_3 (x) dot y_3$),
  ..pin(0, 1, [$A_0$], color: theme.series.at(0), align: left),
  ..pin(1, 3, [$A_1$], color: theme.series.at(1), align: top + right),
  ..pin(2, 2, [$A_2$], color: theme.series.at(2), align: top),
  ..pin(3, 4, [$A_3$], color: theme.series.at(3), align: bottom + right),
)

And we have our final polynomial which effectively goes through all the points:

#import "helpers.typ": *

#let l0(x) = (x - 1) * (x - 2) * (x - 3) / -6 * 1
#let l1(x) = (x - 0) * (x - 2) * (x - 3) / 2 * 3
#let l2(x) = (x - 0) * (x - 1) * (x - 3) / -2 * 2
#let l3(x) = (x - 0) * (x - 1) * (x - 2) / 6 * 4
#let P3(x) = l0(x) + l1(x) + l2(x) + l3(x)

#school(
  xlim: (-0.5, 4), ypos: left, ylim: (-1, 5), legend: top + left,
  curve(P3, 0, 5, label: $P_3 (x)$),
  ..pin(0, 1, [$A_0$], align: left),
  ..pin(1, 3, [$A_1$]),
  ..pin(2, 2, [$A_2$], align: bottom),
  ..pin(3, 4, [$A_3$], align: bottom + right),
)

Generalizing this process, we can write the function $f(x) in PP_n (x)$ that goes through $n+1$ points $A_i (x_i, y_i)$: