# Efficient implementation of the divided difference

Question. I am looking for more efficients ways to compute Newton's divided difference, i.e. the operator $$\partial : f(x,y) \mapsto \frac{f(x,y)-f(y,x)}{x-y} \ ,$$

for the case where $f$ is a polynomial in $x,y$.

Of course the above can be implemented as

 Par = (# - ReplaceAll[#,{x->y,y->x}]) /(x - y) & ;


but for $f$ a polynomial I suspect that there are faster ways to get the same result. Indeed, since the numerator is antisymmetrized in $x,y$ it is divisible by $x-y$, so if $f$ is a polynomial then so is $\partial f$.

It is really this quotient-free form that I am after, and one might expect the division to be relatively time-costly.

NB. The reason I would like to speed this up is that I really want to work with the case of more variables and various $\partial$s, where things start to become slower quite rapidly. See below for an example. For my question, though, the case with two variables suffices.

To get things started let me give one alternative implementation (although it actually appears to be slower than the above Par in examples, see below).

One can check that on monomials we have $$\partial(x^m y^n) = \text{sgn}(m-n) \, (x \, y)^{\min(m,n)} \, \sum_{i=1}^{|m-n|} x^{|m-n|-i} \, y^{i-1} \ ,$$ where $\text{sgn}(0)=0$:

 Table[
Par[x^m y^n] == Sign[m-n] (x y)^Min[m,n] Sum[x^(Abs[m-n] - i) y^(i-1),{i,Abs[m-n]}],
{m, 0, 3}, {n, 0, 3}
] //Simplify //Flatten //Apply[And]
(* Out: True *)


par[expr_Plus] := par /@ expr


Next we reduce to monomials in $x,y$:

par[term_] /; FreeQ[term, x] && FreeQ[term,y] := 0 (* discard terms not involving x or y *)
par[a___ b_ c___] /; FreeQ[b, x] && FreeQ[b, y] := b par[a c] (* extract factors independent of x,y *)
par[term_Times] := par @ Expand @ term (* expand remaing products involving x,y *)


To extract factors of $x\,y$ we can use the comment of AccidentalFourierTransform to improve my previous code:

par[x^m_. y^n_.] := Sign[m-n] (x y)^Min[m,n] par[x^Abs[m-n]]


where I use $\partial(y^n) = -\partial(x^n)$. Finally,

par[y^n_.] := -par[x^n]
par[x^m_.] := Sum[x^(m - i) y^(i - 1), {i,m}]


completes the definition.

Bonus question. Is there a nicer way to implement the above, which in particular gets rid of the need to distinguish between $x$ and higher powers of $x$, and likewise for $y$?

The example that I actually have in mind relates to the Hecke algebra which acts on the space of polynomials in $N$ variables $z_i$ as follows. Writing $\partial_i$ for the divided difference with respect to $z_i$ and $z_{i+1}$, the generators act as $T_i = -(q \,z_i-q^{-1}\,z_{i+1})\,\partial_i + q\$ for $1\leq i\leq N-1$: these obey the braid relations along with the 'Hecke condition' $(T_i - q)(T_i+q^{-1})=0$. In any case, we define

Par[i_] := (# - ReplaceAll[#,{z[i] -> z[i+1],z[i+1]->z[i]}]) /(z[i] - z[i+1]) & ;
T[i_] := -(q z[i] - q^-1 z[i + 1]) Cancel @ Par[i] @ # + q # &


where I note that Factor is slower; Simplify and FullSimplify are much slower still and don't seem to remove all denominators. The above alternative now becomes

par[i_][expr_Plus] := par[i] /@ expr (* reduce to terms; the Expand avoids getting stuck with e.g. (z_i + q)(2 z_i +q) *)
par[i_][term_] /; FreeQ[term, z[i]] && FreeQ[term, z[i+1]] := 0 (* reduce to terms involving z_i or z_{i+1} *)
par[i_][a___ b_ c___] /; FreeQ[b, z[i]] && FreeQ[b, z[i+1]] := b par[i][a c](* reduce to monomials in z_i,z_{i+1} *)
par[i_][term_Times]:= par[i] @ Expand @ term (* expand products still involving z_i,z_{i+1} *)
par[i_][z[j_]^m_. z[k_]^n_.] /; j==i && k==i+1 := Sign[m-n] (z[i] z[i+1])^Min[m, n] par[i][z[i]^Abs[m - n]](* remove factors of (z_i z_{i+1}) *)
par[i_][z[j_]^m_.] := (-1)^(i-j) Sum[z[i]^(m-n) z[i+1]^(n-1), {n, m}]


and

t[i_] := -(q z[i] - q^-1 z[i + 1]) par[i] @ # + q # &


Now I want to compute things like $T_{j-1} \cdots T_2 \, T_1 \, z_1^r$. On my laptop (Win 10 x64, Mma v11.3) I get

j = 10; r=8;
z^r // Composition @@ Table[ T[k], {k,j-1,1,-1}] // AbsoluteTiming ;
z^r // Composition @@ Table[ t[k], {k,j-1,1,-1}] // AbsoluteTiming ;
{%[], %%[], %[] == %%[] // Factor}
Clear[j,r]
(* Out:
{8.59344, 13.2408, True}
*)

• I'm not sure why you find your code inefficient. Your very first piece of code, Simplify[Par[...]] returns the simplified polynomial in under 0.1 seconds, for random polynomials of degree $x^{10}y^{10}$, with around 90 terms. For larger polynomials, it takes more time, but it is still pretty fast. How complex are your polynomials? What degree do you have in mind? – AccidentalFourierTransform Jun 20 '18 at 15:47
• @AccidentalFourierTransform I guess the extra computation time in my case comes from having multiple variables, which for the purpose of this question we can treat as having coefficients involving lots of additional arbitary parameters – Jules Lamers Jun 20 '18 at 18:40
• And thanks for the suggestion of the default patterns, I didn't know about those! – Jules Lamers Jun 20 '18 at 18:42
• How many terms do your polynomials have, on average? – AccidentalFourierTransform Jun 20 '18 at 18:52
• Yes, I assume so :-) the problem is that for all the polynomials I've tried, your code was rather efficient. It would be useful if you could include a representative example of the kind of polynomials you're working with. – AccidentalFourierTransform Jun 20 '18 at 20:03

Okay, here goes a wild attempt at using linear algebra to deal with the problem. The idea is to work with the CoefficientList of the polynomial instead of the symbolic expression. In addition, I wish to turn the CoefficientList of the polynomial in $x,y$, which will be a 2D List, into a vector representing the polynomial:

demoPoly = Sum[c[i, j] x^i y^j, {i, 0, 8}, {j, 0, 8}];
polyVector = Flatten[CoefficientList[demoPoly, {x, y}], 1];


The position of $c_{i, j}$ in the vector is (8 + 1)*i + j + 1.

Now, we know that if $m > n$, then $\partial(x^m y^n) = \sum_{i=1}^{m-n} x^{m-i} y^{n+i-1}$. We also know that if $m < n$, we get the same, except for a minus sign. Thus, we can find the matrix of the operator $\partial$ in the space spanned by monomials $x^m y^n$ with maximum exponent $\alpha-1$:

dMatrix[maxExpo_] := Block[
{alpha, dim, mat, j, j2, sign, min, diff},
alpha = maxExpo + 1;
dim = alpha^2;
mat = ConstantArray[0, {dim, dim}];
Do[
j = alpha m + n + 1;
Do[
mat[[alpha (m - i) + n + i, j]] = 1
, {i, 1, m - n}
];
j2 = alpha n + m + 1;
mat[[All, j2]] = -mat[[All, j]];
, {n, 0, maxExpo - 1}, {m, n + 1, maxExpo}
];
mat
]


The nice thing about this matrix is that it can be applied to lots of different polynomial vectors as long as the largest exponent doesn't exceed $\alpha - 1$.

We can make a helper function that applies this to a symbolic polynomial and returns a symbolic polynomial:

fromCoefficientList[cl_] := Expand[Fold[FromDigits[Reverse[#1], #2] &, cl, {x, y}]

divDiff[maxExpo_][poly_] := fromCoefficientList[
Partition[dMatrix[maxExpo].Flatten[
CoefficientList[poly, {x, y}], {maxExpo + 1, maxExpo + 1}]
, 1], maxExpo + 1]]


The AbsoluteTiming of Expand@Cancel[Par@demoPoly] is

0.088310

and the AbsoluteTiming of divDiff[demoPoly] is only

0.006377

Let me know if this helps you at all! I thought it was a fun problem anyways :)

• I should also add that it's probably not that hard to generalize this to you actual use case, by simply making the variables $x,y$ in the helper function arguments to the helper function instead :) Then you should probably take dMatrix outside the helper function, since it only needs to be computed once. – Marius Ladegård Meyer Jun 20 '18 at 21:36
• Thanks for your answer, that is a nice idea! I will need some time to play around with this a little, but hope to have time for that soon. I'll get back to you when I have, but for now thanks again! – Jules Lamers Jun 22 '18 at 21:14
• No problem! I may add that I tested my approach vs. your Par for Sum[z^5 // Composition @@ Flatten@Table[ Cancel@*T[k], {m, em, 1, -1}, {k, j[[m]] - 1, m, -1}], {j, Subsets[Range[en], {em}]}] with em = 4, en = 8, and clocked in at 6 seconds compared to your 22 on my machine, so there is at least some improvement for your actual use-case :) I haven't tested this, but my gut is that most of those 6 seconds are spent going back and forth between symbolic poly and CoefficientList, but don't know a way around yet... – Marius Ladegård Meyer Jun 22 '18 at 22:01