The Newton interpolation formula, with more variables Alain Lascoux ^{⋆}

Newton addressed the question of transforming discrete sets of data, say the positions of planets at different times, into algebraic functions, before submitting them to the differential calculus that he had just constructed.

Of course, anybody faced with the sequence 1, 2, 4, 8, 16, … will exclaim 32, except for the unwise Khalif who pledged to fill the squares of a chessboard with wheat, doubling the number of grains from one square to the other. To recognize the function f(n) = an^{2} + bn + c , n = 1, 2, …, which is not much more complicated, one already needs to have recourse to finite differences, that is, to evaluate the functions g(n) := f(n) − f(n − 1) and iterate. Indeed, polynomials can be characterized by the fact that iterating finite differences ultimately produces the constant sequence 0, 0, 0, . . ..

However, comets are not likely to appear at regularly spaced lapses of time, and to handle their seemingly erratic apparitions, Newton found the solution of normalizing differences of positions by the interval of time to which they correspond.

In other words, he writes, starting from the table f(t_{1}), f(t_{2}), f(t_{3}), … the normalized differences

[1, 2] = \frac{f(t_{1}) − f(t_{2})}{t_{1} − t_{2}} , [2, 3] = \frac{f(t_{2}) − f(t_{3})}{t_{2} − t_{3}} , [3, 4] = \frac{f(t_{3}) − f(t_{4})}{t_{3} − t_{4}} …

then the differences

[1, 2, 3] = \frac{[1, 2] − [2, 3]}{t_{1} − t_{3}} , [2, 3, 4] = \frac{[2, 3] − [3, 4]}{t_{2} − t_{4}} , …

Since already at the first step one obtains functions of two variables, the proper way to interpret Newton’s operations is to define divided differences ∂_{i} as operators on functions of, say b_{1}, b_{2}, b_{3}, …:

∂_{i} : f ⟶ \frac{f − f^{s_{i}}}{b_{i} − b_{i+1}} , i = 1, 2, …,

denoting by s_{i} the transposition of b_{i} and b_{i+1}.

Divided differences are “local” (they act only on two variables at a time), rational (the coefficients 1/(b_{i} − b_{i+1}) are rational), and decrease degree by 1. More general operators of the same kind can be found in [LS2].

⋆ text written during the Conference Applications of the Macdonald Polynomials, at the Newton Institute in April 2001.

They provide a discrete version of differential calculus, and one would have avoided two centuries of controversy, initiated by Bishop Berkeley, by sticking to them instead of introducing such puzzling entities as epsilons vanishing at different orders.

Given a function of one variable, let

f^{∂} , f^{∂∂} , f^{∂∂∂} , …

denote

∂_{1}(f(b_{1})) , ∂_{2}(∂_{1}(f(b_{1}))) , ∂_{3}(∂_{2}(∂_{1}(f(b_{1})))) , …

respectively.

Then the formula of Newton to interpolate f(t) from its values at time b_{1}, b_{2}, … is f(t) =

f(b_{1})+f^{∂} (t−b_{1})+f^{∂∂} (t−b_{1})(t−b_{2})+f^{∂∂∂} (t−b_{1})(t−b_{2})(t−b_{3})+⋯ (⋆)

the expansion being exact if n + 1 points b_{i} are used and f is a polynomial of degree ≤ n.

Modernists will object that the set of polynomials

{t^{0}, (t − b_{1}), (t − b_{1})(t − b_{2}), …}

is triangular in the usual basis {t^{0}, t^{1}, …}, and therefore, there is no mystery in the fact that it is a linear basis of polynomials. However, what is remarkable in Newton’s formula is the process by which he obtained the coefficients of his polynomials (t − b_{1}) ⋯ (t − b_{n}).

When b_{1}, b_{2}, … all collapse to 0, then (⋆) becomes Taylor’s formula

f(t) = f(0) + f′(0) t + f′′(0)/2! t^{2} + ⋯,

the factorials in the denominator being clear by putting b_{i} = iϵ and letting ϵ tend to 0.

Since divided differences are the discrete analogues of derivatives, the question is :

What is the discrete analogue of Taylor’s formula in several variables ?

In other words, what are the coefficients of the images of f(a_{1}, a_{2}, …) under the possible different divided differences (evaluated in b_{1}, b_{2}, b_{3}, …) in the expansion of f(a_{1}, a_{2}, …) ?

The answer, that I obtained with Marcel-Paul Schützenberger some years ago (cf. for example [LS1]), is as simple as in the case of one variable. The universal coefficients, called Schubert polynomials due to their relevance to geometry, can be defined recursively as follows:

1) They are polynomials Y_{v}(𝔸, 𝔹) in two infinite totally ordered sets of variables 𝔸 = {a_{1}, a_{2}, …}, 𝔹 = {b_{1}, b_{2}, …}, indexed by vectors v ∈ ℕ^{∞} (from now on, we shall suppose that they have only a finite number of non-zero components).

2) They are globally stable under divided differences in the a_{i}’s, that is, the image of a Schubert polynomial is either 0 or a Schubert polynomial.

3) When v is weakly decreasing (one says v is dominant, or is a partition, deleting the terminal zeros), then

Y_{v} = ∏_{(i,j)∈Diagr(v)}(a_{i} − b_{j}) ,  (⋄)

product on all boxes in the diagram of v: one stacks v_{1}, v_{2}, v_{3}, … boxes in successive rows, packing them to the left. A box in position (i, j) gives a factor (a_{i} − b_{j}).

4) Given v ∈ ℕ^{∞} and i ∈ ℕ such that v_{i} > v_{i+1}, writing vs_{i} for [v_{1}, …, v_{i−1}, v_{i+1}, v_{i} − 1, v_{i+2}, …], one has

∂_{i}(Y_{v}) = Y_{vs_{i}} .  (⋄⋄)

One does not write what happens when v_{i} ≤ v_{i+1}, because in that case Y_{v} is the image of some Y_{v′} under ∂_{i}, (precisely v′ = [v_{1}, …, v_{i−1}, v_{i+1} + 1, v_{i}, v_{i+2}, …]). The vanishing ∂_{i}^{2} = 0 implies that ∂_{i}(Y_{v}) = 0 in that case.

On the indices of Schubert polynomials, the action of divided differences amounts to sorting (and decreasing), in other words, amounts to use the symmetric group. Thus to check whether the equations (⋄⋄) are consistent essentially reduces to check that

∂_{i}∂_{i+1}∂_{i}(Y_{v}) = ∂_{i+1}∂_{i}∂_{i+1}(Y_{v})  (♣)

for any i, any v, because s_{i}s_{i+1}s_{i} = s_{i+1}s_{i}s_{i+1}, and s_{i}s_{j} = s_{j}s_{i}, |j − i| ≠ 1 are the defining equations of the symmetric group, as generated by simple transpositions s_{i}. There is no need to check ∂_{i}∂_{j}(Y_{v}) = ∂_{j}∂_{i}(Y_{v}) because operations on two disjoint pairs of variables clearly commute.

Relation (♣) is in fact a property of divided differences, substituting in it Schubert polynomials is irrelevant . One indeed has

∂_{i}∂_{i+1}∂_{i} = ∂_{i+1}∂_{i}∂_{i+1} , i ≥ 1  (♣♣)

To check this, one can without loss of generality take i to be 1. Functions of a_{1}, a_{2}, a_{3}, … can be expressed as linear combinations of the six monomials

1, a_{1}, a_{2}, a_{1}^{2}, a_{1}a_{2}, a_{1}^{2}a_{2} ,

with coefficients which are functions symmetrical in a_{1}, a_{2}, a_{3} and arbitrary in a_{4}, a_{5}, . . ..

Now, the divided differences ∂_{1}, ∂_{2} commute with multiplication by symmetric functions in a_{1}, a_{2}, a_{3}, and thus the identity ∂_{1}∂_{2}∂_{1} = ∂_{2}∂_{1}∂_{2} needs only to be tested on the above set. However, both ∂_{1}∂_{2}∂_{1} and ∂_{2}∂_{1}∂_{2} decrease the degree by 3. Thus it is only necessary to check, and we leave this pleasure to the reader, that ∂_{1}∂_{2}∂_{1}(a_{1}^{2}a_{2}) and ∂_{2}∂_{1}∂_{2}(a_{1}^{2}a_{2}) are the same constant (which is 1), the other images being null for reason of degree.

To get a general v, one can start from different possible dominant vectors. Because the set of partitions is a lattice (with respect to intersection and union of diagrams), compatibility of the different choices boils down to the following lemma. Lemma. Given a dominant v ∈ ℕ^{∞}, and n ∈ ℕ such that v ≤ ρ := [n − 1, n − 2, …, 1, 0, 0, …] (componentwise), then there exists at least one chain of divided differences such at each step a ∂_{i} is applied to a dominant Schubert polynomial of the type (a_{i} − b_{j}) g, with g symmetrical in a_{i}, a_{i+1} (and therefore the image is a dominant Schubert polynomial corresponding to a diagram with the box of coordinate (i, j) erased, the factor (a_{i} − b_{j}) becoming 1 under ∂_{i}).

The proof of the lemma consists in peeling off boxes of a staircase diagram, in such a way as to have, at a step where ∂_{i} will be used, columns i and i + 1 of lengths differing by 1. This can be realized by erasing boxes from top to bottom, in successive diagonals, as indicated by the following example for v = [6, 6, 6, 2, 2, 0, 0, …], n = 9 (boxes are numbered in the order they are peeled off).

1

8 2
13 9 3
10 4
14 11 5

6

12 7

Now, if Y_{ρ} was the product ∏_{i,j,i+j≤n}(a_{i} − b_{j}), then Y_{v}, as obtained by the above process, will be the product of factors (a_{i} − b_{j}) corresponding to the boxes which are left. Therefore, the definition of Schubert polynomials is consistent.

Newton’s polynomials could be characterized by their vanishing properties, because they are written in terms of their roots. Similarly, when v is dominant, and different from [0, 0, …], one can also easily write the vanishing properties of Y_{v}.

Indeed we shall only need the property that the specialization a_{1} = b_{1}, a_{2} = b_{2}, a_{3} = b_{3}, … in Y_{v}(𝔸, 𝔹), which we shall write as Y_{v}(𝔹, 𝔹), vanishes (all the factors in the main diagonal of the diagram vanish).

More generally, one can prove (this will be the only property in this text that we admit. Otherwise, we would need to appeal to properties of the Ehresmann-Bruhat order on the symmetric group, for which we refer to Macdonald [M]), that

v ∈ ℕ^{∞}, v ≠ [0, 0, …] ⇒ Y_{v}(𝔹, 𝔹) = 0 .  (⋆⋆)

The relation ∂_{1}∂_{2}∂_{1} = ∂_{2}∂_{1}∂_{2} shows that products of divided differences are not independent. We shall only use the following products (which are, in fact, canonical representatives of all possible products of divided differences). Given K = [k_{1}, k_{2}, …] ∈ ℕ^{n}, let

∂^{k} := (∂_{1} ⋯ ∂_{k_{1}}) (∂_{2} ⋯ ∂_{k_{2}+1}) (∂_{3} ⋯ ∂_{k_{3}+2}) …,

that is, ∂^{k} consists in products of blocks of consecutive divided differences, of respective lengths k_{1}, k_{2}, …, starting with ∂_{1}, ∂_{2}, … respectively, divided differences operating on their left.

We can now state the multivariate Newton interpolation formula. For any polynomial in a_{1}, a_{2}, …, one has :

f(𝔸) = ∑_{K∈ℕ^{∞}} f^{∂^{K}} (𝔹) Y_{K}(𝔸, 𝔹) .  (⋆ ⋆ ⋆)

Proof. One has to test (⋆ ⋆ ⋆) on a linear basis of polynomials in 𝔸 (with coefficients in anything, for example functions of 𝔹). One takes the Schubert polynomials {Y_{J}} as such a basis. The question is : Is is true that

Y_{J}(𝔸, 𝔹) = ∑_{K∈ℕ^{∞}}(Y_{J})^{∂^{K}} (𝔹) Y_{K}(𝔸, 𝔹) .  (??)

But the Y_{J}^{∂^{K}} are either zero, or Schubert polynomials, and they all vanish under the specialization 𝔸 = 𝔹, except for Y_{[0,0,...]} = 1. Now, Y_{J}^{∂^{K}} can be Y_{[0,0,...]} only when k_{1} + k_{2} + ⋯ = j_{1} + j_{2} + ⋯, for degree reason. By recursion on the righmost non-zero component of J, say j_{ℓ}, one sees that k_{i} = 0, i > ℓ and k_{ℓ} = j_{ℓ}. Finally K = J and equation (??) becomes the irrefutable identity

Y_{J}(𝔸, 𝔹) = Y_{[0,0,...]}(𝔹, 𝔹) Y_{J}(𝔸, 𝔹) .

This proves (⋆ ⋆ ⋆) in full generality.

The original Newton formula has not been lost. When f depends only on a_{1}, then f^{∂^{K}} is 0 if K is not of the type K = [n, 0, 0, …] for some integer n. In that case, Y_{K} is dominant, corresponding to a row diagram with n boxes, and indeed the Newton polynomial is

Y_{[n,0,0,...]} = (a_{1} − b_{1}) (a_{1} − b_{2}) ⋯ (a_{1} − b_{n}) .

References

[LS1] A. Lascoux, M.P. Schützenberger, Interpolation de Newton à plusieurs variables, Séminaire d’Algèbre M.P.Malliavin 1983-84, Springer L.N. 1146 (1985) 161-175.

[LS2] A. Lascoux, M.P. Schützenberger, Symmetrization operators on polynomial rings, Funct.Anal. 21(1987) 77-78.

[M] I.G. Macdonald, Notes on Schubert polynomials, Publications du LACIM, Montréal, (1991).

[N] I. Newton, Philosophiæ Naturalis Principia Mathematica, Liber III, Lemma V, London (1687).

C.N.R.S., Institut Gaspard Monge Université de Marne-la-Vallée,

5 Bd Descartes, Champs sur Marne, 77454 Marne La Vallée Cedex 2 FRANCE

                                    Alain.Lascoux@univ-mlv.fr
                           http://phalanstere.univ-mlv.fr/∼al