Alternative ways to implement a triangular recursion

Question

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?

Answer 1

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
  ]

triangular 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.

Answer 2

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

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)

$$ \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]

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

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] &

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

• How to deal with the condition $u_i=u_{i+1}$ in B-Spline basis function?
• Calculating an integral by the Romberg Algorithm
• A problem about recursion formula of de Casteljau algorithm
• How to deal with recursion formula in Mathematica?

Answer 3

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.

Answer 4

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:

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}}]  *)

Answer 5

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...

Answer 6

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.

Answer 7

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.