Tetration

Gottfried Helms - Univ Kassel 07' 2007 - 2008

Eigendecomposition of triangular matrix-operators (here: U_t)

An iterative procedure for Pari/GP (pseudocode)

Version 2. Some changes in formulas were made for convenience, which resulted in different signs-convention for the A-matrices and denominators (no overall effect of changes)

Gottfried Helms,11.1.2010
first version: 9.8.2008

Intro

The reason for the computation of the eigendecomposition of the U_t-matrices is, that fractional iterates of U_t-tetration (x->t^x -1) can be computed using fractional powers of the matrix(-operator) U_t.

This concept is not new; it is known in the context of Bell- and Carleman-matrices for iteration of functions. However, the diagonalization method is only one nethod to make use of these matrices to arrive at fractional iteration as meaningful interpolation of the integer iteration. Another method involves the use of the binomial expansion instead.

The general question of a method-independent, external appropriateness-criterion for the selection of an interpolation-method for fractional iteration of the exponential-function (exp_b(x)) or of the "decremented exponential-function" (dxp_t(x) or U_t(x)as I call it here) is not yet answered and the decision for the "unique" appropriate interpolation-method is thus still open.

An "external" argument could be, that the interpolations for fractional iterates should be invariant across different fixpoint-shifts. I should explain this a bit.

Assume the usual tetrational function

            T_b(x) = b^x
            T_b°¹(x) = T_b(x)
            T_b°^h(x) = T_b°^h-1(T_b(x))

then T_b(x)has inifinitely many fixpoints t₀,t₁,t₂,... (almost all are complex), such that

T_b°^h(t_k) = t_k or: b^t^k = t_k or: b = t_k^1/t^k

Then we may translate the problem

T_b^oh(x) = (U_tk^oh(x/t_k – 1 ) +1)*t_k .

The advantage of this fixpoint-shift is then, that the powerseries for U_tk(x) has no constant term (then matrix-operator for U_tk() is triangular) such that it can be –for instance – inverted. Also fractional iterates can then formally be computed using matrix-logarithm or diagonalization.

The problem with this is that it is possible (there was made an argument in a discussion), that these fixpoint-shifts give different results using different fixpoints. If this is indeed the case, then a unique definition of tetration must include the convention which fixpoint has to be taken (besides of the question, whether this method is appropriate at all for the complete mathematical context).

However; the matrix-based methods and also diagonalization has been, and is, in use in the research, so I leave this problem open and I'm going to discuss the method for the cases where no fixpoint-shift is needed further to exhibit some deeper details.

The diagonalization gives consistent results in the sense that, for instance,

U_t^o0.5(U_t^o0.5(x)) = U_t^o0.5+0.5(x) = U_t^o1(x)

or more generally

U_t^oa(U_t^ob(x)) = U_t^oa+b(x) .

as far as always the same method is applied.

Here I give technical details, how to compute the diagonalization matrices (eigen-matrices) for such triangular matrix-operators and some sample values for comparision, which I haven't seen elsewhere.

Further links: a more detailed description can be found at

https://go.helms-net.de/math/tetdocs/ContinuousfunctionalIteration.pdf

A different approach, based on polynomial interpolation which however gives the same result can be found at

https://go.helms-net.de/math/tetdocs/ExponentialPolynomialInterpolation.pdf

For an intro in tetration see also [1]; a list of valuable references can be found in [2] and [3].

Eigendecomposition and Schröder-function

If we denote

V(x)~ * U_t^h = V(U_t^oh(x))~

(see all matrix-definitions in appendix 1) then fractional powers of U_t can be computed by fractional powers of the diagonal matrix D of the eigen-decomposition of U_t:

V(x)~ * U_t^h = V(x)~ * (W * D^h*W^-1) = V(U_t^oh(x))~

where D is the diagonal-matrix of eigenvalues and W is the triangular matrix of eigenvectors. (Note, that W and W^-1 are only unique besides column-scaling of W if they are assumed to be triangular)

In my computations I use a column-normed version of W having unit-diagonal; then W seems to be also a matrix-operator. It provides then the coefficients of the Schröder-function σ which decomposes the iteration of a function into

U_t^oh(x) = σ^-1(u^h*σ(x))

(see a short definition of the Schröder-function at [1](superlogarithm); Schröders article at [4])

In the light of Schröder's result, that

lim_h->inf f°^h(x)/f'(0)^h = σ(x)

which is also

lim_h->inf f°^h(x)/u^h = σ(x)

and where we find the relevant 2'nd column of the matrix W just being that limit, the assumed property, that W is in fact a matrix-operator, is then evident.

The matrix D is diagonal, containing the consecutive powers of u = log(t).

Note: the given eigen-decomposition can only be computed, if t≠exp(1), u≠1, so we discuss only iterations of U_t(x) = t^x –1 ≠ exp(x) – 1here.

The Pari/GP-pseudocode

Besides numerical matrices Pari/GP can handle symbolic coefficients in a basic way, so we can feed the following procedure with the U_t-matrix containing u as unevaluated symbol as well to get the symbolic solution. However, the symbolic description of terms is much consumptive with increasing index; matrix-size of 16x16 gives the result in a blink; 32x32 needs a bit of time, for 64x64 matrices computation time and memory requirements are already excessive (see a zoomed image coefficients-matrix at a₁₆ in appendix 2)

The input for the following procedure is some triangular matrix (here U_t, where u=log(t)≠1 ). The matrices W and W^-1 are initialized with the unit-diagonal and then new entries are computed using the previously computed ones. Here in W we can compute the columns and in W^-1 the rows separately, always beginning at the diagonalelement (which is 1) proceeding to the element at position r,c vertically resp. horizontally right-to-left.

{ TRI_EW(U) =

local(u=U[2,2], W, D, Wi, t=exp(u));

dim=rows(U);

D=diag(U); \\ for operator-matrices U_t this is 1,u,u²,u³,...

W=matid(dim);

for(c=1,dim-1,

for(r=c+1,dim,

W_r,c=sum(k=c,r-1,U_r,k*W_k,c) / (D_c-D_r) ))

Wi=matid(dim); \\ means matrix W^-1

for(r=2,dim,

forstep(c=r-1,1,-1,

Wi_r,c=sum(k=0,r-1-c,Wi_r,r-k*U_r-k,c) / (D_r-D_c) ))

return([[u,t,exp(u/t)],W,D,Wi]); }

Coefficients of the powerseries for fractional iterates

Then the powerseries for U_t(x) = t^x-1 and its h'th iterate U_t°^h(x) due to the eigendecomposition

U_t^h = W * ^dV(u^h) * W^-1

The denominators d_k are products

d_k =

The a_h;k are polynomials in u and u^h; in the following the a_h;k are expressed as matrices of coefficients A_k of these polynomials:

A₁ =		=> a_h,1= 1u^h* d₁ = 1
A₂ =		=> a_h;2 = (1 u^h – 1u^2h) d₂= (u-1)
A₃ =		=> a_h;3 = (1* u^h – 3u^2h+ 2u^3h) + u(2* u^h – 3u^2h+ 1u^3h) = (1 + 2u) u^h – (3 + 3u) u^2h + (2 + 1u) u^3h d₃= (u-1)(u²-1)
A₄ =		d₄= (u-1)(u²-1)(u³-1)
A₅ =
A₆ =
A₇ =
A₈ =

The number of rows at a_h;k are binomial(k-1,2)+1, the number of columns is k+1.

The highest exponent to which a power of u appears is thus binomial(k-1,2) due to the rows plus k*h due to the columns.

For instance in a_h;8 for height h=1 is thus (binomial(7,2))+8*h=21+8 = 29 and for h=2 is 37.

However, for the complete coefficient at x^k this is increased by the coefficient u^k-1 and also decreased by exponents of u in the denominator d₈; here the highest exponent is binomial(k,2)=28 so the remaining most significant power of u at x⁸ is u^{21+8*h+(8-1)-28} which is u⁸ for h=1 and u¹⁶ for h=2 . In general, the ratio of most significant powers of u in numerator and denominator for each k is simply u^k*h .

If u is in the near of 1, and h is fractional, then additional hyperbolic effects occur since the denominators approach zero if integer powers uⁿ approach 1, and these are not evenly cancelled by the numerators if h is noninteger.

For integer heights the coefficients a_h;k are evenly divisible by the denominators and only few different powers of u are involved in the final coefficient a_h;k/d_k; for fractional h this is not the case. In that case we have a much more complicated structure of powers of u (however the most significant power is still u^k*h.)

The explicite computation of the coefficients a_k is

a_h;k = V(u)~*A_k*V(u^h)=

where the index of matrices are beginning at 1 as used in Pari/GP

The coefficients-matrices for heights 1 and 2

Heights h introduce higher powers of u into the coefficients. For the coefficients-matrices A_k this means an increasing rowshift of each column according to height h. Let now D denote the coefficients at powers of u of the denominator, which are independent of the height-parameter. Then if the A_k matrices are rewritten according to integer heights, they occur to be linear combinations of (shifted) D-vectors. This shifting expresses multiplications with powers of u.

The rows of the following matrices agree to the resp. powers of u, so I added the accordingly vectors of u-powers in the most-left column to make this explicit. Since the row-shifting "renormalized" the powers of u depending on the h-parameter, we can simply use the row-sums to arrive at the final values for the scalar coefficient a_h,k.

	A₄ dep. on h	rowsums	decomposition into multiples of D₄	a_h,k =
~ *				= V(u)~ * D₄ 0*
~*				= u * V(u)~ * D₄ [1]*
~*				= u² V(u)~ D₄ [1, 7u, 6u², 1u³]*

The coefficients for the Schröder-function

The coefficients of the powerseries for the Schröder-function for U_t(x) , σ_Ut , occur in the second column of W, so that

V(x)~ * W = V(σ_Ut(x))~

and for the inverse

V(x)~ * W^-1 = V(σ_Ut^-1(x))~

so that finally

            V(x)~               * W                 = V(σ_Ut(x))~
            V(σ_Ut(x))~        * ^dV(u^h)           = V(u^h * σ_Ut(x))~
            V(u^h * σ_Ut(x))~ * W^-1              = V(σ_Ut^-1(u^h σ_Ut(x)))~
                                                            = V(U_t^°h(x))~

The powerseries is

The denominators d_k are the same as the in the previous

d_k =

The s_k are polynomials in u only; in the following they are expressed as rowvectors of their coefficients (S_k ), where the leftmost entry is associated with u⁰:

S₁ =		=> s₁= 1 d₁ = 1
S₂ =		=> s₂ = 1 d₂= (u-1)/u
S₃ =		=> s₃ = (1 + 2u) d₃= (u-1)(u²-1)/u²
S₄ =
S₅ =
S₆ =
S₇ =

The number of columns at s_k are binomial(k-1,2)+1, so the highest exponent to which a power of u appears in s_k is thus binomial(k-1,2).

Note, that these coefficients occur also in the corresponding A-matrices in the second column, which agrees to the fact, that the schröder-function is assumed to be the limit of the case h->inf: then all other columns of that A-matrices get shifted to infinity and –if |u|<1 – vanish by approximation of u^h to zero.

We have from this

s_k = S_k *V(u) =

To see the ratio of the highest exponent of u in numerator and denominator we add binomial(k-1,2)+1 and (k-1) according to the main formula and subtract binomial(k,2) according to the denominator which comes out to be zero for each k.

The coefficients for the inverse of the Schröder-function are in the same way taken from W^-1 instead.

Appendix 1: Definitions for matrix-notation

For U-tetration the matrices U and U_tare the triangular Stirling-matrix scaled as the following:

U = ^dF^-1*St2 * ^dF // performs x - >exp(x)-1
U_t = ^dV(u)* U // where u=log(t); performs x -> t^x– 1

where the involved (infinite size assumed) matrices are (only the top-left segment is displayed):

V(u) = ^dV(u) = diag(V(u))	columnvector([1, u, u²,u³,u⁴,...]) diagonal arrangement
^dF =	diag(0!,1!,2!,3!,...)
St2 =	Stirlingnumbers 2^nd kind
^dF^-1St2 ^dF = U =
U_t =	the numerical value for u=log(t) is actually used here