Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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]

Bézier curve truncated with the de Casteljau algorithm

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?

share|improve this question
    
Ideally I'd be expecting an answer demonstrating all three examples... –  J. M. 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. –  J. M. Feb 14 '12 at 23:07

5 Answers 5

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] :=
 BezierCurve@Transpose@MapThread[triangular2[fb[u], {##}] &, pts]

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.

share|improve this answer
    
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.

share|improve this answer
    
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. –  J. M. 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.

share|improve this 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...

share|improve this answer

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.

share|improve this answer

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.