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

7 Answers 7

In this answer, I will use the Functional Paradigm to deal with triangular recursive formula in a uniform manner.

For the triangular recursive formula

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

In general, $f(x)=a x+b$, so the triangular recurisive formula can be denoted as below:

$$T_k^{(n)}=\alpha(k,n) T_{k-1}^{(n)}+\beta(k,n)T_{k-1}^{(n+1)}$$

Algorithm Schematic Diagram

enter image description here

Gerneral Solution for Triangular Recursions

1.When $\alpha(k,n)=\alpha(k)$, and $\beta(k,n)=\beta(k)$

$$ \begin{align} & T_k^{(n)}=\alpha(k)T_{k-1}^{(n)}+\beta(k)T_{k-1}^{(n+1)} \\ & \qquad =\left[\alpha(k)+\beta(k)\right]T_{k-1}^{(n)}+\beta(k)\left[T_{k-1}^{(n+1)}-T_{k-1}^{(n)}\right] \end{align}$$

Code Template

 triangular[ini_] :=
  FoldList[
    Most@#1 (α[#2] + β[#2]) + β[#2] Differences@#1 &,
     ini, Range@(Length@ini - 1)]

2.When $\alpha(k,n)=\alpha(n)$, and $\beta(k,n)=\beta(n)$

$$ T_k^{(n)}=\alpha(n)T_{k-1}^{(n)}+\beta(n)T_{k-1}^{(n+1)}\\ =\left(T_{k-1}^{(n)},T_{k-1}^{(n+1)}\right)\cdot \left(\alpha(n),\beta(n)\right) $$

Code Template

triangular[ini_] :=
 NestList[
  Dot @@@
   (Thread@
     With[{sub = Partition[#, 2, 1]},
     {sub, Table[{α[n], β[n]}, {n, 1, Length@sub}]}]) &, 
  ini, Length@ini - 1]

3.When $\alpha(k,n)=\alpha$, and $\beta(k,n)=\beta$, where $\alpha,\beta$ are constant

$$ \begin{align} & T_k^{(n)}=\alpha T_{k-1}^{(n)}+\beta T_{k-1}^{(n+1)} \\ & \qquad =\left[\alpha+\beta \right]T_{k-1}^{(n)}+\beta \left[T_{k-1}^{(n+1)}-T_{k-1}^{(n)}\right] \end{align}$$

Code Template

triangular[ini_] :=
 FoldList[
  Most@#1 (α + #2) + #2 Differences@#1 &, ini, 
  PadRight[{β}, Length@ini - 1,β]]

Examples

Case 1 for $\alpha(k)$, and $\beta(k)$ (Romberg Algorithm)

enter image description here

$$ \frac{-1}{4^1-1}T_{n}(f)+\frac{4^1}{4^1-1}T_{2n}(f)=S_n(f) \\ \frac{-1}{4^2-1}S_{n}(f)+\frac{4^2}{4^2-1}S_{2n}(f)=C_n(f) \\ \frac{-1}{4^2-1}C_{n}(f)+\frac{4^3}{4^3-1}C_{2n}(f)=R_n(f) \\ $$

where, $ T_n=\frac{h}{2}\left[ f(a)+ 2\sum_{i=1}^{n-1}f(x_i) +f(b)\right]\\ $ and $h=\frac{b-a}{n}$

Now I denote $T_{2^i}=T_0^{(i)},S_{2^i}=T_1^{(i)},C_{2^i}=T_2^{(i)},R_{2^i}=T_3^{(i)}$, then I can achieve the formula as bleow

$$T_k^{(n)}=\frac{-1}{4^k-1} T_{k-1}^{(n)}+\frac{4^k}{4^k-1} T_{k-1}^{(n+1)}$$

Implementation for Romberg Algorithm

trapezium[func_, n_, {a_, b_}] :=
 With[{h = (b - a)/n},
  Module[{f = Function[x, func@x, Listable]}, 
   h (Total@f@Range[a, b, h] - 1/2 (f@a + f@b))]]

RombergAlgorithm[ini_] :=
 FoldList[
  Most@#1 + (4^#2)/(4^#2 - 1) Differences@#1 &, ini, Range@3]

Testing for Romberg Algorithm

rombergIni = trapezium[Exp, 2^(# - 1), N[{0, 1}, 15]] & /@ Range@6;
TableForm[
 Flatten[RombergAlgorithm[rombergIni], {{2}, {1}}], 
 TableHeadings -> {2^Range[0, 5], {Subscript["T", "n"], 
  Subscript["S", "n"], Subscript["C", "n"], Subscript["R", "n"]}}, 
 TableAlignments -> Center]

enter image description here

Case 2 for $\alpha(n)$, and $\beta(n)$ (Bernoulli numbers)

$$a_{n+1}^m=(m+1)\left[a_{n}^m- a_{n}^{m+1}\right]$$

In the same way, I denote $a_{n+1}^m=T_k^{(n)}$, then I can achieve

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

where, $T_0^{(n)}=\frac{1}{n+1}$

Inplementataion for Akiyama-Tanigawa Algorithm

AkiyamaTanigawaAlgorithm[ini_] :=
 NestList[
  Dot @@@
  (Thread@
    With[{sub = Partition[#, 2, 1]},
     {sub, Table[{i, -i}, {i, 1, Length@sub}]}]) &, ini, Length@ini - 1]

Testing for Akiyama-Tanigawa Algorithm

TableForm[
 Flatten[AkiyamaTanigawaAlgorithm[1/# & /@ Range@9], {{2}, {1}}], 
 TableAlignments -> Center] // TraditionalForm

enter image description here

Case 3 B-spline basis function (C. de Boor Algorithm)

Let $\vec{U}=\{u_0,u_1,\ldots,u_m\}$ a nondecreasing sequence of real numbers,i.e, $u_i\leq u_{i+1}\quad i=0,1,2\ldots m-1$

$$N_{i,0}(u)= \begin{cases} 1 & u_i\leq u<u_{i+1}\\ 0 & otherwise \end{cases} $$ $$N_{i,p}(u)=\frac{u-u_i}{u_{i+p}-u_i}N_{i,p-1}(u)+\frac{u_{i+p+1}-u}{u_{i+p+1}-u_{i+1}}N_{i+1,p-1}(u) $$

Denoting $N_{i,p}=T_p^{(i)}$, the formula can be shown as below:

$$T_p^{(i)}=\frac{u-u_i}{u_{i+p}-u_i}T_{p-1}^{(i)}+\frac{u_{i+p+1}-u}{u_{i+p+1}-u_{i+1}}T_{p-1}^{(i+1)} $$

Seeing my answer for B-Spline basis function

Case 4 for $\alpha$, and $\beta$ are constants(deCasteljau Algorithm)

$$\vec{P}_{k,i}(u_0)=(1-u_0)\vec{P}_{k-1,i}(u_0)+u_0\vec{P}_{k-1,i+1}(u_0)$$

Simularlily, I denote $\vec{P}_{k,i}(u_0)=T_k^{(n)}$, then I can achieve

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

Inplementataion for deCasteljau Algorithm

deCasteljauAlgorithm[ini_, u0_] :=
 FoldList[
  Most@#1 + #2 Differences@#1 &, ini, PadRight[{u0}, Length@ini - 1, u0]]

Testing for deCasteljau Algorithm

Column@
 Flatten[
  deCasteljauAlogrithm[{{0, 0}, {2, 4}, {4, 5}, {6, 0}}, 2/5], {{2}, {1}}] // 
Style[#, 12, FontFamily -> Times] &

enter image description here

Specifically, when $\alpha+\beta=1$,

deCasteljau[ini_, u0_] := 
  NestList[
   MovingAverage[#, {1 - u0, u0}] &, pts, Length@ini - 1]

deCasteljau[{{0, 0}, {2, 4}, {4, 5}, {6, 0}}, 2/5]

Same result

For more detailed information, seeing my demonstration

Generating a Bezier Curve by the de Casteljau Algorithm

Summary

In this semester course, I am always encourtering many triangular recursive formula, so I deal with them one by one in functional paradigm with the help of @xzczd, @Mr.Wizard,@Michael E2 etc.


My questions about triangular recursions

share|improve this answer

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

Here's general approach, which I primarily show in order to give an answer that includes the Neville-Aitken algorithm. It peculiarly works from the bottom of the triangle up, that is $T_k^{(n)}$ or t[k, n] are generated in the order shown in the table:

Mathematica graphics

One of the distinctions to make clear is whether the function $f$ in the recursion

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

is permitted to depend explicitly on $n$ or $k$. It is convenient to do that. Extra information can be helpful. On the other hand, dependence only on the two inputs $T$ seems "purer" (forgive me if my background in CS is not sufficient to come up with accurate technical jargon). It is the purer version I present. Below f is function that takes two inputs $T$ and returns the next $T$; the function initfn takes an integer $n$ and constructs the initial object $T_0^{(n)}$. It is used to delay the construction of each object until needed. The integer k is determines which $T_k^{(0)}$ is returned.

tri[f_, initfn_, k_] :=
 With[{initial = Range[k, 0, -1]},
  Last@Fold[
    FoldList[
      f,
      initfn @ #2,
      #1
      ] &,
    {initfn @ First @ initial},
    Rest @ initial]
  ]

The objects generated by FoldList are shown in the rows of the table above. Comparing with Szabolcs', this one steps down $n$ and Szabolcs' steps up $k$.

Neville-Aitken interpolation

The difficulty seems to be coming up with a data structure for $T_k^{(n)}$ that allows the recursion rule $f$ to be written in terms of $T$ as above. What I came up with was

{xi, xj, yij}

which correspond to $x_i$, $x_j$, and $p_{i,j}(x)$ in the notation of the Wikipedia article. Translating the recursion formula in the article becomes straightforward.

Clear[na, nax, nainit]; 
na[x_] :=  (nax[{ax1_, ax2_, ay_}, {bx1_, bx2_, by_}] :=
             {ax1, bx2, ((ax1 - x) by + (x - bx2) ay)/(ax1 - bx2)}; nax);
nainit[{x_, y_}] := {x, x, y};

pts = {{0, 1}, {1/3, 1}, {1/2, 1/2}, {3/4, -1}, {1, 0}};

Last @ tri[na[1/4], nainit[pts[[# + 1]]] &, Length[pts] - 1]
(*  161/160  *)

In this application, in which the points are all given beforehand, using an initfn is inconvenient and saves no memory. So it's not a good example to demonstrate when initfn is advantageous. One can easily rewrite tri to strip out initfn and k and make tri a function of initial.

Akiyama-Tanigawa

at[n_] := Last@
  tri[{#1[[1]], #1[[1]] (#1[[2]] - #2[[2]])} &,
  {# + 1, 1/(# + 1)} &,
  n];

at[10]
(*  5/66  *)

de Casteljau

dC[BezierCurve[pts_?MatrixQ, opts___], u_?NumericQ] := BezierCurve@
  tri[
   Join[{u #2[[1]] + (1 - u) #1[[1]]}, #2] &,
   {pts[[# + 1]]} &,
   Length[pts] - 1];

dC[BezierCurve[{{0, 0, 0}, {1, 1, 1}, {2, -1, 1}, {3, 0, 2}}], 1/3]
(*  BezierCurve[{{1, 2/9, 20/27}, {5/3, 0, 10/9}, {7/3, -(2/3), 4/3}, {3, 0, 2}}]  *)
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.