Lagrange's interpolating polynomial is straightforward, but not terribly efficient; every time we change points or x we have to redo the entire calculation. Ideally we'd like to be able to perform:

1. Some computation on points that would let us efficiently evaluate the fitted polynomial for different values of x in O(n) time, or

2. A computation on a particular x that would let us efficiently add new points to the set we use to generate the interpolating polynomial.

"Neville's algorithm" lets us generate the same interpolating polynomial recursively. By flipping the recursion around and generating values from the bottom up, we can achieve goal #2 and add new points incrementally.

Neville's Algorithm

Start the recursion with a single point. Any point Loading... has a unique 0th order polynomial passing through it - the constant function Loading.... For points Loading..., Loading..., let's call this Loading..., Loading..., etc.

Loading... is the unique FIRST order polynomial (i.e., a line) going through points Loading... and Loading....

this first recursive step gives us this rule:

For higher order terms like Loading..., let's call Loading... 'P_l', and Loading... 'P_r' (the polynomial fitted through the left and right set of points).

Similarly, the left and rightmost inputs - Loading... and Loading... - will be Loading... and Loading....

Neville's algorithm states that:

This recurrence works because the two parents Loading... and Loading... already agree at all points except Loading... and Loading....

Tableau-based Methods

Neville's algorithm generates each new polynomial from Loading... and Loading..., using this recursion to incorporate the full set of points.

You can write these out these relationships in a "tableau" that grows to the right from an initial left column:

The next few functions will discuss "rows" and "columns" of the tableau. That refers to the rows and columns of this representation;

The first column here is the initial set of points. Each entry in each successive column is generated through some operation between the entry to its left, and the entry one left and one down.

Look again at Neville's algorithm:

Loading... refers to the entry in the same row, previous column, while Loading... is one row lower, previous column.

If each cell in the above tableau tracked:

• the value of P(x) for the cell
• Loading..., the x value of the leftmost point incorporated so far
• Loading..., the right point

we could build up Neville's rule incrementally. Let's attempt to build a function of this signature:

First, write a function to process each initial point into a vector that contains each of those required elements:

Next, a function that generates the next entry, given l and r:

We can use higher-order functions to turn this function into a NEW function that can transform an entire column:

neville-tableau will generate the entire tableau:

Really, we're only interested in the first row:

p0 p01 p012 p0123 p01234

So define a function to grab that:

the final piece we need is a function that will extract the estimate from our row of Loading... vectors:

Putting it all together:

How do we know this works? We can prove it by using generic arithmetic to compare the full symbolic lagrange polynomial to each entry in the successive approximation.

Generic Tableau Processing

The above pattern, of processing tableau entries, is general enough that we can abstract it out into a higher order function that takes a prepare and merge function and generates a tableau. Any method generating a tableau can use a present function to extract the first row, OR to process the tableau in any other way that they like.

This is necessarily more abstract! But we'll specialize it shortly, and rebuild neville-incremental into its final form.

I'm keeping points in the argument vector for now, vs returning a new function; if you want to do this yourself, curry the function with (partial tableau-fn prepare merge present).

Redefine neville-merge to make it slightly more efficient, with baked-in native operations:

And now, [[neville]], identical to [[neville-incremental]] except using the generic tableau generator.

The form of the tableau also makes it easy to select a particular column instead of just the first row. Columns are powerful because they allow you to successively interpolate between pairs, triplets etc of points, instead of moving onto very high order polynomials.

I'm not sure it's the best interface, but we'll add that arity here.

Modified Neville

Press's Numerical Recipes, chapter 3 (p103) ( http://phys.uri.edu/nigh/NumRec/bookfpdf/f3-1.pdf ) describes a modified version of Neville's algorithm that is slightly more efficient than the version above.

Allan Macleod, in "A comparison of algorithms for polynomial interpolation", discusses this variation under the name "Modified Neville".

By generating the delta from each previous estimate in the tableau, Modified Neville is able to swap one of the multiplications above for an addition.

To make this work, instead of tracking the previous Loading... estimate, we track two quantities:

• Loading... is the delta between Loading... and Loading..., i.e., Loading....
• Loading... is the delta between Loading... and Loading..., i.e., Loading....

We can recover the estimates generated by the original Neville's algorithm by summing C values across the first tableau row.

Equation 3.1.5 in Numerical recipes gives us the equations we need:

These equations describe a merge function for a tableau processing scheme, with state == [x_l, x_r, C, D].

Let's implement each method, and then combine them into final form. The following methods use the prefix mn for "Modified Neville".

tableau-fn allows us to assemble these pieces into a final function that has an interface identical to neville above. The implementation is more obfuscated but slightly more efficient.

Folds and Tableaus by Row

NOTE: These folds and scans seem to be higher performance than the functions above. Prefer *-sum functions when you want to consume a full sequence of points and get the full input, and *-scan functions when you want to observe all intermediate estimates. If you want a column of the tableau, stick with the versions above.

The advantage of the method described above, where we generate an entire tableau and lazily pull the first entry off of each column, is that we can pass a lazy sequence in as points and get a lazy sequence of successive estimates back. If we don't pull from the result sequence, no computation will occur.

One problem with that structure is that we have to have our full sequence of points available when we call a function like neville. As you pull an element off of the returned sequence, it just-in-time generates a new diagonal of the tableau required to realize that number.

Even if the sequence is lazy, this is limiting. What if we want to pause, save the current estimate and pick up later where we left off?

Look at the tableau again:

p0 p01 p012 p0123 p01234 p1 p12 p123 p1234 . p2 p23 p234 . . p3 p34 . . . p4 . . . . . . . . . . . . . . . . . . .

If you stare at this for a while, you might notice that it should be possible to use the merge and present functions we already have to build the tableau one row at a time, given ONLY the previous row:

This method reverses the order of the points, since rows are built from the bottom up:

The order of the points is reversed, but this is fine; polynomial interpolation doesn't care about the order of points. (NOTE that this WILL be something we have to consider in the fold version of Richardson extrapolation, in emmy.polynomial.richardson!)

Notice that the /diagonal/ of this tableau is identical to the top row of the tableau before the points were reversed.

Here's something close, using our previous merge and prepare definitions:

There's a problem here. reduce only returns the FINAL value of the aggregation:

(let [f (generate-new-row* prepare present)] (f [p1 p12 p123 p1234] [x0 fx0])) ;; => p01234

If we want to continue building rows, we need the entire new row! Lucky for us, Clojure has a version of reduce, called reductions, that returns each intermediate aggregation result:

(let [f (generate-new-row prepare present)] (f [p1 p12 p123 p1234] [x0 fx0])) ;; => [p0 p01 p012 p0123 p01234]

Quick aside here, as we've stumbled across a familiar pattern. The discussion above suggests the idea of a "fold" from functional programming: https://en.wikipedia.org/wiki/Fold_(higher-order_function)

Folds are explored (and implemented!) in detail in the [[emmy.algebra.fold]] namespace. See that namespace for an introduction, but I'll add a reminder here.

A fold consists of:

• init, an function that returns an initial piece of state called an "accumulator"

• a binary merge function that combines ("folds") a new element x into the accumulator, and returns a value of the same shape / type as the accumulator returned by init.

• a present function that transforms the accumulator into a final value.

In Clojure, you perform a fold on a sequence with the reduce function:

For example:

Our generate-new-row function from above is exactly the merge function of a fold. The accumulator is the latest tableau row:

init == a function that returns [], the initial empty row.

present == a function similar to neville-present or mn-present that simply sums up the estimate deltas for the entire row, instead of returning the running tally.

present only needs to return the final value in each row because that is the best current estimate given all points supplied so far.

If you want to recover the previous behavior of a lazy sequence of all estimates, the lazy "scan" pattern of emmy.util.aggregate/scan allows you to observe each element of the diagonal of the tableau as it's generated. This is identical to the "first row" of the non-fold tableau.

Now that we've identified this new pattern, we can rewrite generate-new-row to return a new function matching the fold interface described in [[emmy.algebra.fold]]. The new function is called [[tableau-fold-fn]]:

Next, we can use this to generate specialized fold functions for our two incremental algorithms above - neville and modified-neville.

Instead of using neville-present as above, neville-fold returns only the final estimate. This is because folds are meant to consume their entire input sequence. If you want to observe intermediate values as they're generated, you can use the "scan" pattern, implemented in neville-scan.

Instead of using mn-present as above, modified-neville-fold uses the dynamically bound emmy.util.aggregate/*fold* to sum all deltas and return the current best estimate taking all points into account.

If you want to observe intermediate values as they're generated, you can use the "scan" pattern, implemented in neville-scan.

Fold Utilities

af/fold->scan will return a function that acts identically to the non-fold, column-wise version of the interpolators. It does this by folding in one point at a time, but processing EVERY intermediate value through the presentation function.

Using this function, we specialize to our two incremental methods.

Next, check out:

• rational.cljc to learn how to interpolate rational functions
• richardson.cljc for a specialized implementation of polynomial interpolation, when you know something about the ratios between successive x elements in the point sequence.