Alternative ways to implement a triangular recursion

Triangular recursions are a class of algorithms that frequently turn up in computational mathematics. These recursions are expressible in the general form

$$T_k^{(n)}=f(T_{k-1}^{(n)},T_{k-1}^{(n+1)})$$

for some function $f$ and initial values $T_0^{(n)}$. The "triangular" adjective is easily visualized if the $T_k^{(n)}$ are arranged in an array:

$$\begin{array}{}T_0^{(0)}&T_1^{(0)}&T_2^{(0)}&T_3^{(0)}\\T_0^{(1)}&T_1^{(1)}&T_2^{(1)}&\\T_0^{(2)}&T_1^{(2)}&&\\T_0^{(3)}&&&\end{array}$$

Due to the nice theory behind triangular recursions, it is known that such algorithms can be implemented using only a one-dimensional scratch array instead of a two-dimensional array.

So much for theory. I shall now give various examples of triangular recursions that occur in computational practice, as implemented in Mathematica, with the underlying triangular recursion flanked by (* ------ *) comment lines.

Here is the Akiyama-Tanigawa algorithm for the Bernoulli numbers, with the redefinition $B_1=\frac12$:

myBernoulliB[n_Integer] := Module[{atArray = 1/Range[n + 1]},

(* ------ *)
Do[
atArray[[j]] = j (atArray[[j]] - atArray[[j + 1]]),
{k, n + 1}, {j, k - 1, 1, -1}];
(* ------ *)

First[atArray]]

myBernoulliB /@ Range[0, 20]
{1, 1/2, 1/6, 0, -(1/30), 0, 1/42, 0, -(1/30), 0, 5/66}

BernoulliB[Range[0, 20]]
{1, -(1/2), 1/6, 0, -(1/30), 0, 1/42, 0, -(1/30), 0, 5/66}


Here's the Neville-Aitken algorithm for polynomial interpolation:

iPolyVal[pts_?MatrixQ, x_] := Module[{n = Length[pts] - 1, xa, ya, temp},
{xa, ya} = Transpose[pts];

(* ------ *)
Do[
temp = (x - xa[[j]])/(x - xa[[k]]) - 1;
ya[[j]] = ya[[j + 1]] + (ya[[j + 1]] - ya[[j]])/temp,
{k, n + 1}, {j, k - 1, 1, -1}];
(* ------ *)

First[ya]]

iPolyVal[{{0, 1}, {1/3, 1}, {1/2, 1/2}, {3/4, -1}, {1, 0}}, 1/4]
161/160

InterpolatingPolynomial[{{0, 1}, {1/3, 1}, {1/2, 1/2}, {3/4, -1}, {1, 0}}, 1/4]
161/160


Here's (a simplified version of) de Casteljau's algorithm for splitting Bézier curves:

bezierChop[BezierCurve[pts_?MatrixQ, opts___], u_?NumericQ] :=
Module[{n = Length[pts]},
BezierCurve[Transpose[
Function[{vec}, Block[{ta = vec},

(* ------ *)
Do[
ta[[j]] = u ta[[j + 1]] + (1 - u) ta[[j]],
{k, n - 1}, {j, k, 1, -1}];
(* ------ *)

ta]] /@ Transpose[pts]]]] /; 0 < u < 1

Graphics3D[{AbsoluteThickness[5],
BezierCurve[{{0, 0, 0}, {1, 1, 1}, {2, -1, 1}, {3, 0, 2}}],
Directive[Blue, AbsoluteThickness[3]],
Translate[
bezierChop[
BezierCurve[{{0, 0, 0}, {1, 1, 1}, {2, -1, 1}, {3, 0, 2}}],
1/3], {0, 0, 1/10}]}, Boxed -> False]


There are many more examples of situations that make use of triangular recursions, such as Romberg quadrature, divided differences for interpolating polynomials, the Levin transformation for summing series, the Cox-de Boor algorithm for B-splines... etc.

Having shown that they are important, note that the common core of the algorithms that use triangular recursions is a double-index Do[] loop working on a preset scratch array, modifying elements as needed.

That should be sufficient preamble. My question now is: are there alternative methods (e.g. non-procedural methods) for implementing triangular recursions?

-
Ideally I'd be expecting an answer demonstrating all three examples... –  Ｊ. Ｍ. Feb 13 '12 at 14:41
Apparently people are having a hard time with Neville-Aitken. It's probably a good thing I had chosen Akiyama-Tanigawa and de Casteljau as the other examples; Cox-de Boor will apparently make people weep. –  Ｊ. Ｍ. Feb 14 '12 at 23:07

The following is a naive but general implementation of the recursion formula

$$T_k^{(n)}=f(T_{k-1}^{(n)},T_{k-1}^{(n+1)})$$

triangular[f_, initial_] :=
First@Nest[
f @@@ MapIndexed[Join, Partition[#, 2, 1]] &,
initial,
Length[initial] - 1
]


traingular takes a function f[tn0, tn1, n] where tn0 corresponds to $T_{k-1}^{(n)}$, tn1 corresponds to $T_{k-1}^{(n+1)}$ and n corresponds to $n$.

Then the Bernoulli number function can be implemented as

myBernoulliB[n_] := triangular[#3 (#1 - #2) &, 1/Range[n+1]]


We can write it a bit more readably as

triangular[Function[{Tn0, Tn1, n}, n (Tn0 - Tn1)], 1/Range[5]]


The other two algorithms can also be written in terms of triangular[], but I don't think this is really better than your Do loop. (triangular[] would need to be extended to pass $k$ to the function as well.)

This can be generalized, based on the same principle, to return multiple results (${T_0^{(3)}, T_1^{(2)}, T_2^{(1)}, T_3^{(0)} }$ above) like this:

triangular2[f_, initial_] := Module[{tag, a, b},
{{a}, {b}} = Reap[
Nest[
Function[arg,
Sow[Last[arg], tag];
f[##, Last[arg]] & @@@ MapIndexed[Join, Partition[arg, 2, 1]]
],
initial,
Length[initial] - 1
],
tag
];
Append[b, a]
]


Then bezierChop can be implemented using:

fb[u_] := Function[{tn0, tn1}, u tn1 + (1 - u) tn0];

bezierChop[BezierCurve[pts_?MatrixQ, opts___], u_?NumericQ] :=


To make this easier to understand using the original triangular function, MapThread[triangular[fb[u], {##}]&, pts] would return only the first point out of all points included in the resulting BezierCurve.

-
This does not use a two-dimensional array, but it does keep re-allocating new and new arrays behind the scenes which is pretty much the same. –  Szabolcs Feb 13 '12 at 14:47
It's not the same, because you don't need to have a 2D array in memory at the same time. So it now depends on how smart the heap manager is, but at least you give it a chance to have less simultaneous memory consumption. So, with a smart heap manager, you trade memory for run-time, which seems to be the objective of the question. +1. –  Leonid Shifrin Feb 13 '12 at 15:33

I don't know if this can be done for all cases, but for the first example (Bernoulli numbers), here is a possible functional implementation:

Clear[myBernoulliBF];
myBernoulliBF[n_Integer] :=
With[{atArray = 1/Range[n + 1], len = n + 1},
First@Fold[
(# + Range[len]*Append[Differences[#], 0]) *
UnitStep[Range[len] - #2] - Range[len]*Append[Differences[#], 0] &,
atArray,
Reverse[Range[len]]]];


It is functional because everything is immutable, and I don't have to ever change the original atArray variable.

Perhaps, similar methods can be used for other cases.

-
Looks good, but I'm still holding out for a general technique, as the Akiyama-Tanigawa algorithm isn't the only triangular method I use. –  Ｊ. Ｍ. Feb 13 '12 at 13:54

Yet another implementation of the general expression for term $T_k^{(n)}$ using ReplaceRepeated is as follows:

term[k_, n_, func_, initial_] := Module[{t},
t[i_, 1] := initial[[i]];
t[k, n] //. t[i_, j_] /; j != 1 :> func[t[i, j - 1], t[i + 1, j - 1]]
]


Then, the final step is merely term[1, n, func, initial]. The modified BernoulliB function can be implemented as:

ClearAll[myBernoulliB];
SetAttributes[myBernoulliB, Listable]
myBernoulliB[n_Integer] := Module[{initial = 1/Range[n + 1]},
t[i_, 1] := initial[[i]];
t[1, n + 1] //. t[i_, j_] /; j != 1 :> i (t[i, j - 1] - t[i + 1, j - 1])
]

myBernoulliB[Range[10]]
(*Out[1]= {1/2, 1/6, 0, -(1/30), 0, 1/42, 0, -(1/30), 0, 5/66}*)


However, do note that ReplaceRepeated is the devil incarnate when it comes to recursive replacements, and since the above approach walks backwards from $T_k^{(n)}$ all the way to $T_0^{(i)}$, the initial values, it is computationally inefficient for anything other than moderately small values of n. In general, it is better to walk forward like in Szabolcs's answer.

-

One more way of implementing recurrences is RecurrenceTable. For the first and simplest of your examples its application is trivial:

RecurrenceTable[{a[0, m] == 1/(m + 1),
a[n + 1, m] == (m + 1) (a[n, m] - a[n, m + 1])},
a, {n, 0, 10}, {m, 0, 10}
][[All, 1]]

{1, 1/2, 1/6, 0, -(1/30), 0, 1/42, 0, -(1/30), 0, 5/66}


I spent some time trying to implement Neville's algorithm but did not succeed. Although I feel that it can be done in this way too. This is a good challenge...

-

I couldn't leave until I posted something so here it is. Just pieces and ideas for you to work with, and relating to the first example only at this point.

Nest[
#2[[1]]*(#[[1]] - #[[2]]) & ~MapIndexed~ Partition[#, 2, 1] &,
1/Range[n + 1],
n
]


Or:

f[{x_, y_, n_}] := (x - y) n

Nest[f /@ Transpose[{Most@#, Rest@#, Range[Length@# - 1]}] &, 1/Range[n + 1], n]

Fold[f /@ Transpose[{Most@#, Rest@#, Range@#2}] &, 1/Range[n + 1], Range[n, 1, -1]]


These last two could potentially be quite fast with a Compiled f.

-