Carleman Linearization and Systems of Arbitrary Depth Polynomial Recursions

Mikołaj Myszkowski

Physics Department, University of Oxford, Oxford, UK.

DOI: 10.4236/alamt.2022.121001                

Abstract

New approach to systems of polynomial recursions is developed based on the Carleman linearization procedure. The article is divided into two main sections: firstly, we focus on the case of uni-variable depth-one polynomial recurrences. Subsequently, the systems of depth-one polynomial recurrence relations are discussed. The corresponding transition matrix is constructed and upper triangularized. Furthermore, the powers of the transition matrix are calculated using the back substitution procedure. The explicit expression for a solution to a broad family of recurrence relations is obtained. We investigate to which recurrences the framework can be applied and construct sufficient conditions for the method to work. It is shown how introduction of auxiliary variables can be used to reduce arbitrary depth systems to the depth-one system of recurrences dealt with earlier. Finally, the limitations of the method are discussed, outlining possible directions for future research.

Keywords

Polynomial Recursion, Carleman Linearization, Transfer Matrices

1. Introduction

Recurrence relations arise in various fields of mathematics and prove to be a powerful tool for studying probability [1] and combinatorics [2]. It is often the case that construction of a solution to a mathematical problem simplifies to solving a recurrence relation. The theory of linear recurrences is well studied, leading to various results such as the celebrated Binet’s formula or the characteristic polynomial method [3] [4]. The dynamics of non-linear recurrence sequences are of special interest, as there are various examples of non-linear (in particular, polynomial) maps exhibiting chaotic behaviour. The examples include the logistic, quadratic and exponential maps, as well as the Ricatti recurrence [5] - [10].

Definition 1.1. [11] Let $u_i^k\in ℂ\mathrm{,}i\ge 0$ be k sequences of complex numbers such that for every $i\ge n$ the following relation holds:

$(\begin{array}{l}{u}_i^1=F_1\left(u_{i-1}^1,\cdots ,u_{i-1}^k,u_{i-2}^1,\cdots ,u_{i-2}^k,\cdots ,u_{i-n}^1,\cdots ,u_{i-n}^k\right)\\ u_i^2=F_2\left(u_{i-1}^1,\cdots ,u_{i-1}^k,u_{i-2}^1,\cdots ,u_{i-2}^k,\cdots ,u_{i-n}^1,\cdots ,u_{i-n}^k\right)\\ ⋮\\ u_i^k=F_k\left(u_{i-1}^1,\cdots ,u_{i-1}^k,u_{i-2}^1,\cdots ,u_{i-2}^k,\cdots ,u_{i-n}^1,\cdots ,u_{i-n}^k\right)\end{array}$ (1.1)

Then, we say that sequences $u_i^k\in ℂ$ are a solution to the depth-n system of recurrences defined by functions $F_i$.

In case when functions $F_i$ are linear, the system of Equations (1.1) can be rewritten in the matrix form [12]. If the auxiliary vectors $y^i\mathrm{,}i\ge n$ and $e_j$ are defined as:

$y^i={\left[1,u_i^1,\cdots ,u_i^k,u_{i-1}^1,\cdots ,u_{i-1}^k,u_{i-2}^1,\cdots ,u_{i-2}^k,\cdots ,u_{i-n+1}^1,\cdots ,u_{i-n+1}^k\right]}^{\text{T}}$ (1.2)

$e_j=\left[{\delta }_{0,j},{\delta }_{1,j},\cdots ,{\delta }_{nk,j}\right]$ (1.3)

then the system of recurrences (1.1) can be rewritten as recurrence relation for the sequence of vectors $y^i$ :

$y^i=M{y}^{i-1}$ (1.4)

where $\left(1+nk\right)\times \left(1+nk\right)$ matrix M is called the transition matrix associated with (1.1). By induction, the solution to (1.4) is given by powers of matrix M:

$y^i=M^{i+1-n}{y}^{n-1}$ (1.5)

where $y^{n-1}$ is defined by initial conditions. Equation (1.5) can be then multiplied by vectors $e_j$ from the left, hence extracting the original sequences from $y^i$ :

$u_i^k=e_k{M}^{i+1-n}{y}^{n-1}$ (1.6)

The above method of solving linear recurrences is fundamentally equivalent to the method of characteristic polynomials [13].

Even if functions $F_i$ are non-linear (but analytic), the solution to the system still can be formulated in terms of powers of the appropriately chosen transition matrix, in analogy with the linear counterpart [14]. The method of describing dynamical systems (recursions in particular) in terms of infinite-dimensional linear algebra is referred to as Carleman embedding or Carleman linearization [15] [16].

Theorem 1.2. [9] [16] Let $k=1$, and $F_1=F_1\left(u_{i-1}^1\right)$ be analytic in $u_{i-1}^1$. Let the (infinite) auxiliary vectors be defined as:

$y^i={\left[1,u_i^1,{\left(u_i^1\right)}^2,{\left(u_i^1\right)}^3,\cdots \right]}^{\text{T}},e_j=\left[{\delta }_{0,j},{\delta }_{1,j},{\delta }_{2,j},\cdots \right]$ (1.7)

Then, the solution can be written in the form:

$u_i^1=e_1{M}^i{y}^0$ (1.8)

where $y^0$ is specified by initial conditions, and the infinite-dimensional matrix M is given by coefficients of Maclaurin expansion of $F_1$ :

$M_{ab}=\frac{1}{b!}{\frac{{\text{d}}^b}{\text{d}{x}^b}{\left(F_1\left(x\right)\right)}^a|}_{x=0}$ (1.9)

Theorem 1.2 allows for an analytic function to be iterated by computing powers of the infinite transition matrix associated with the function. In practice, however, it is troublesome to give an explicit expression for the powers of matrix M [17]. Nevertheless, in some cases, the transition matrix can be diagonalized and the corresponding eigenvectors can be found [16]. We focus on systems of recurrence relations (1.1) with $F_i$ restricted to be finite degree polynomials. The article is divided into two main sections.

Firstly, the k = 1 (uni-variable) depth-one case is discussed, leading to the exact expression for the solution in terms of the initial condition and coefficients of polynomials $F_i$. As it turns out, in most cases the transition matrix can be diagonalized by a suitable choice of variables.

The method of the Carleman embedding is then applied to the multi-variable depth one recurrence. With minor modifications of the auxiliary vectors, the transition matrix can be compactly written as a block matrix. Again, M is diagonalized and the eigenvectors are found. The main theorem of the article is given, solving the recursion.

Finally, we discuss the general case (1.1). Provided that $F_i$ is polynomial in $u_i^k$, the system of Equations (1.1) is embedded in an infinite linear space. It is shown that by introducing appropriate auxiliary variables, any multi-variable arbitrary depth recurrence can be reduced to multi-variable depth-one recurrence. Therefore, the recurrence (1.1) can be solved using the main theorem presented in Section 3.

Throughout the article, we assume that the functions $F_i$ are polynomial in all variables. For convenience, the first elements of vectors are indexed from 0 instead of 1. A similar rule applies to matrices.

2. Uni-Variable Recurrences of Depth One

In the uni-variable case, there is only one sequence $u_i^1\equiv u_i$ and one function $F_i\equiv F$ (we drop the unnecessary indices for simplicity). Furthermore, F is a function of a single variable. The recurrence relations of this kind can be written as:

$u_i=F\left(u_{i-1}\right)={\displaystyle \underset{j=0}{\overset{m}{\sum }}}{c}_j{\left(u_{i-1}\right)}^j$ (2.1)

where $c_j$ are constant coefficients of the polynomial F. The recurrence (2.1) is of special interest, as it includes many well-known nonlinear chaotic maps, i.e. the logistic map, which is used as an example at the end of this section [18].

Using Theorem 1.2 and the multinomial expansion theorem [19], the infinite dimensional transition matrix can be obtained:

$\begin{array}{c}{T}_{ab}=\frac{1}{b!}{\frac{{\text{d}}^b}{\text{d}{x}^b}{\left(F_1\left(x\right)\right)}^a|}_{x=0}=\frac{1}{b!}{\frac{{\text{d}}^b}{\text{d}{x}^b}{\left({\displaystyle \underset{j=0}{\overset{m}{\sum }}}{c}_j{x}^j\right)}^a|}_{x=0}\\ =\frac{1}{b!}{\frac{{\text{d}}^b}{\text{d}{x}^b}{\displaystyle \underset{{\displaystyle {\sum }_{l=0}^m{k}_l}=a}{\sum }}\left(\begin{array}{c}a\\ k_0,k_1,\cdots ,k_m\end{array}\right){\displaystyle \underset{t=0}{\overset{m}{\prod }}}{\left(c_t{x}^t\right)}^{k_t}|}_{x=0}\end{array}$ (2.2)

where $k_l$ are assumed to be non-negative integers. The condition ${\displaystyle {\sum }_{l=0}^m{k}_l}=a$ denotes the fact that the bigger sum runs only over such non-negative integers $k_l$ that sum up to a. This further simplifies to:

$\begin{array}{c}{T}_{ab}=\frac{1}{b!}{\frac{{\text{d}}^b}{\text{d}{x}^b}{\displaystyle \underset{{\displaystyle {\sum }_{l=0}^m{k}_l}=a}{\sum }}\left(\begin{array}{c}a\\ k_0,k_1,\cdots ,k_m\end{array}\right){\displaystyle \underset{t=0}{\overset{m}{\prod }}}{\left(c_t{x}^t\right)}^{k_t}|}_{x=0}\\ ={\displaystyle \underset{{\displaystyle {\sum }_{l=0}^m{k}_l}=a,{\displaystyle {\sum }_{l=0}^{m}l{k}_l}=b}{\sum }}\left(\begin{array}{c}a\\ k_0,k_1,\cdots ,k_m\end{array}\right){\displaystyle \underset{t=0}{\overset{m}{\prod }}}{c}_t^{k_t}\\ ={\displaystyle \underset{{\displaystyle {\sum }_{l=0}^m{k}_l}=a,{\displaystyle {\sum }_{l=0}^{m}l{k}_l}=b}{\sum }}\frac{a!}{{\displaystyle {\prod }_{s=0}^m{k}_s!}}{\displaystyle \underset{t=0}{\overset{m}{\prod }}}{c}_t^{k_t}\end{array}$ (2.3)

Again, the conditions under the bigger sum narrow down the possible combinations of $k_l$ to those that obey ${\displaystyle {\sum }_{l=0}^m{k}_l}=a$, ${\displaystyle {\sum }_{l=0}^{m}l{k}_l}=b$. If no such integers exist, the sum vanishes and $T_{ab}=0$. The transition matrix (2.3), together with the auxiliary vectors from Equation (1.7) defines the Carleman embedding of the recurrence (2.1) in the infinite dimensional linear space. We now prove several facts that will become useful later on in calculating powers of $T_{ab}$.

Theorem 2.1. Let $T_{ab}$ be the infinite dimensional transition matrix corresponding to the recurrence (2.1). Then, the following holds true:

1) $T_{0b}={\delta }_{0,b}$

2) If $c_0=0$, $T_{ab}$ is upper triangular

Proof. (1) We obtain:

$T_{0b}={\displaystyle \underset{{\displaystyle {\sum }_{l=0}^m{k}_l}=a,{\displaystyle {\sum }_{l=0}^{m}l{k}_l}=b}{\sum }}\left(\begin{array}{c}0\\ k_0,k_1,\cdots ,k_m\end{array}\right){\displaystyle \underset{t=0}{\overset{m}{\prod }}}{c}_t^{k_t}$ (2.4)

However, since $k_l\ge 0$ and ${\displaystyle {\sum }_{l=0}^m{k}_l}=0$, all of the integers $k_l$ have to vanish. Therefore ${\displaystyle {\sum }_{l=0}^{m}l{k}_l}=b=0$. If $b\ne 0$, there are no $k_l$ satisfying the conditions under the sum, and $T_{0b}=0$. If $b=0$ :

$T_{00}=\left(\begin{array}{c}0\\ 0,\cdots ,0\end{array}\right){\displaystyle \underset{t=0}{\overset{m}{\prod }}}{c}_t^0=1$ (2.5)

(2) If $c_0=0$, then the lowest non-zero power of x appearing in polynomial $F\left(x\right)$ is $x^1$. By extension, the lowest power of x appearing in powers $F^a$ is of order a. As a result:

$T_{ab}=\frac{1}{b!}{\frac{{\text{d}}^b}{\text{d}{x}^b}{\left({\displaystyle \underset{j=1}{\overset{m}{\sum }}}{c}_j{x}^j\right)}^a|}_{x=0}=\frac{1}{b!}{\frac{{\text{d}}^b}{\text{d}{x}^b}\left(c_1^a{x}^a+\cdots \right)|}_{x=0}=0$ (2.6)

for $b<a$.

It is worth noting that the above mentioned polynomials have also been studied in the setting of graphs [20]. The fact that the transition matrix can take the upper triangular form makes diagonalization feasible. One is always free to work with new shifted variables ${u^{\prime }}_i=u_i-d$, which satisfy a new recurrence:

${u^{\prime }}_i=F^{\prime }\left({u^{\prime }}_{i-1}\right)={\displaystyle \underset{j=0}{\overset{m}{\sum }}}{c^{\prime }}_j{\left({u^{\prime }}_{i-1}\right)}^j$ (2.7)

The coefficients ${c^{\prime }}_j$ can be obtained by substituting ${u^{\prime }}_i=u_i-d$ into (2.1):

$\begin{array}{c}{u^{\prime }}_i+d={\displaystyle \underset{j=0}{\overset{m}{\sum }}}{c}_j{\left({u^{\prime }}_{i-1}+d\right)}^j={\displaystyle \underset{j=0}{\overset{m}{\sum }}}{c}_j{\displaystyle \underset{l=0}{\overset{j}{\sum }}}\left(\begin{array}{c}j\\ l\end{array}\right)d^l{u^{\prime }}_{i-1}^{j-l}\\ ={\displaystyle \underset{l=0}{\overset{m}{\sum }}}{\displaystyle \underset{j=l}{\overset{m}{\sum }}}{c}_j\left(\begin{array}{c}j\\ l\end{array}\right)d^{j-l}{u^{\prime }}_{i-1}^l={\displaystyle \underset{j=0}{\overset{m}{\sum }}}{c^{\prime }}_j{u^{\prime }}_{i-1}^j+d\end{array}$ (2.8)

Which implies:

${c^{\prime }}_0={\displaystyle \underset{j=0}{\overset{m}{\sum }}}{c}_j{d}^j-d,{c^{\prime }}_l={\displaystyle \underset{j=l}{\overset{m}{\sum }}}{c}_j\left(\begin{array}{c}j\\ l\end{array}\right)d^{j-l}\text{for}l\ne 0$ (2.9)

Therefore, new coefficients ${c^{\prime }}_j$ are functions of old coefficients $c_j$ and the shift parameter d. It turns out that there always exists d such that the transition matrix of the transformed recurrence becomes upper triangular.

Lemma 2.2. Let the recurrence on $u_i$ be defined as in (2.1). For every $F\left(u_{i-1}\right)$ there exists $d\in ℂ$ such that the transition matrix ${T^{\prime }}_{ab}$ corresponding to the new sequence ${u^{\prime }}_i=u_i-d$ is upper triangular.

Proof. The new sequence ${u^{\prime }}_i$ is obtained by performing a shift transformation with a parameter d. As already derived in (2.7), the recurrence of the shifted sequence is:

${u^{\prime }}_i={\displaystyle \underset{j=0}{\overset{m}{\sum }}}{c^{\prime }}_j{u^{\prime }}_{i-1}^j$ (2.10)

where the constant term (term near 0-th power of ${u^{\prime }}_{i-1}$ ) reads:

${c^{\prime }}_0={\displaystyle \underset{j=0}{\overset{m}{\sum }}}{c}_j{d}^j-d$ (2.11)

Theorem 2.1 states that the transition matrix ${T^{\prime }}_{ab}$ corresponding to the shifted recurrence is upper triangular if the constant term ${c^{\prime }}_0=0$ vanishes. This leads to the following polynomial equation in d:

$0={\displaystyle \underset{j=0}{\overset{m}{\sum }}}{c}_j{d}^j-d$ (2.12)

By the fundamental theorem of algebra, there always exists $d\in C$ such that ${c^{\prime }}_0=0$ [21]. Therefore, the transition matrix $T_{ab}$ can be always brought to an upper triangular form by a suitable shift of variables.

As a result, we can consider only the recurrences satisfying $c_0=0$ without loss of generality.

Because the transition matrix is upper triangular, the diagonal elements $T_{aa}$ are its eigenvalues [22]. Using (2.2):

${\lambda }_a=T_{aa}=\frac{1}{a!}{\frac{{\text{d}}^a}{\text{d}{x}^a}{\left({\displaystyle \underset{j=1}{\overset{m}{\sum }}}{c}_j{x}^j\right)}^a|}_{x=0}=\frac{1}{a!}{\frac{{\text{d}}^a}{\text{d}{x}^a}\left(c_1^a{x}^a+\cdots \right)|}_{x=0}=c_1^a$ (2.13)

The transition matrix is diagonalizable if all its eigenvalues are different, i.e. when $c_1^a=c_1^b\Rightarrow a=b$ [23]. Since in general $c_1\in ℂ$, this condition can be rewritten as $c_1\ne 0\wedge \left(\left|c_1\right|\ne 1\vee arg\left(c_1\right)\notin \mathbb{Q}\right)$. It should be noted that the right hand side of Equation (2.11) may have more than one distinct root, which means that there exist recurrences that can be shift transformed from $c_1=0\vee \left(\left|c_1\right|=1\wedge arg\left(c_1\right)\in \mathbb{Q}\right)$ to $c_1\ne 0\wedge \left(\left|c_1\right|\ne 1\vee arg\left(c_1\right)\notin \mathbb{Q}\right)$.

Example 2.3. Let $F\left(x\right)=x^3+2x^2+x$. The recurrence reads:

$u_i=F\left(u_{i-1}\right)={\left(u_{i-1}\right)}^3+2{\left(u_{i-1}\right)}^2+u_{i-1}$ (2.14)

The linear coefficient $c_1=1$, hence there is no guarantee that the corresponding transition matrix is diagonalizable. Nevertheless, by performing the shift transformation with $d=-2$, we obtain (using formulas (2.9)):

${u^{\prime }}_i={\left({u^{\prime }}_{i-1}\right)}^3-4{\left({u^{\prime }}_{i-1}\right)}^2+5{u^{\prime }}_{i-1}$ (2.15)

with ${c^{\prime }}_0=0,\left|{c^{\prime }}_1\right|\ne 1\ne 0$.

The above example shows that the condition $c_1\ne 0\wedge \left(\left|c_1\right|\ne 1\vee arg\left(c_1\right)\notin \mathbb{Q}\right)$ is sufficient, but not necessary for the diagonalization of $T_{ab}$. Some transition matrices that are hard to deal with can be diagonalized after shift of variables, thus the recurrence can be solved.

The next two Lemmas regard the procedure of diagonalization of infinite upper triangular matrices. The set of linearly independent eigenvectors is obtained. Even though the matrix is infinite, the eigenvector corresponding to a-th eigenvalue ${\lambda }_a$ has at most $a+1$ non-zero components. An explicit expression for the inverse of an upper triangular matrix is given (provided that it exists), which allows for calculation of the inverse of the corresponding modal matrix. The transition matrix is then diagonalized, and the powers of the matrix are subsequently calculated.

Lemma 2.4. Let $M_{ab}$ be an infinite upper triangular matrix with $\forall a$ ${\lambda }_a=M_{aa}\ne 0$ and different eigenvalues, i.e. ${\lambda }_a={\lambda }_b\Rightarrow a=b$. Let $v_a^b$ denote the b-th component of the (not necessary normalized) eigenvector corresponding to the a-th eigenvalue ${\lambda }_a$. Then:

1) $\forall b>a$, $v_a^b=0$ ;

2) $\forall b=a$, $v_a^b=1$ ;

3) $\forall l_0<l_{p+1}$ :

$v_{l_{p+1}}^{l_0}={\displaystyle \underset{p=0}{\overset{l_{p+1}-l_0-1}{\sum }}}{\displaystyle \underset{l_0<l_1<\cdots <l_p<l_{p+1}}{\sum }}{\left(-1\right)}^{p+1}{\displaystyle \underset{j=0}{\overset{p}{\prod }}}\frac{M_{l_j,l_{j+1}}}{M_{l_j,l_j}-{\lambda }_{l_{p+1}}}$ (2.16)

Proof. The eigenvector equation $M{v}_a={\lambda }_a{v}_a$ in the matrix form reads:

$\left[\begin{array}{ccccccc}{M}_{00}-{\lambda }_a& M_{01}& M_{02}& M_{03}& M_{04}& M_{05}& \cdots \\ 0& \ddots & ⋮& ⋮& ⋮& ⋮& \\ 0& \cdots & M_{a-1,a-1}-{\lambda }_a& M_{a-1,a}& M_{a-1,a+1}& M_{a-1,a+2}& \cdots \\ 0& \cdots & 0& 0& M_{a,a+1}& M_{a,a+2}& \cdots \\ 0& \cdots & 0& 0& M_{a+1,a+1}& M_{a+1,a+2}& \cdots \\ ⋮& & ⋮& ⋮& ⋮& ⋮& \ddots \end{array}\right]\left[\begin{array}{c}{v}_a^0\\ ⋮\\ v_a^{a-1}\\ v_a^a\\ v_a^{a+1}\\ ⋮\end{array}\right]=0$ (2.17)

Lets consider a finite upper triangular system of equations that arises from truncating the above matrix at $\left(a+1\right)\times \left(a+1\right)$ :

$\left[\begin{array}{cccc}{M}_{00}-{\lambda }_a& M_{01}& \cdots & M_{0a}\\ 0& \ddots & ⋮& ⋮\\ 0& \cdots & M_{a-1,a-1}-{\lambda }_a& M_{a-1,a}\\ 0& \cdots & 0& 0\end{array}\right]\left[\begin{array}{c}{v}_a^0\\ ⋮\\ v_a^{a-1}\\ v_a^a\end{array}\right]=0$ (2.18)

Since all of the eigenvalues are different, the diagonal elements $\forall b\ne a$ ${\lambda }_b-{\lambda }_a\ne 0$ and the matrix M are diagonalizable. The determinant of (2.16) vanishes, and thus the solution is not unique (up to a scaling factor), as expected. The explicit expression for $v_a^b$ as a function of $v_a^a$ can be obtained by the back substitution method [24]:

$v_{l_{p+1}}^{l_0}=v_{l_{p+1}}^{l_{p+1}}{\displaystyle \underset{p=0}{\overset{l_{p+1}-l_0-1}{\sum }}}{\displaystyle \underset{l_0<l_1<\cdots <l_p<l_{p+1}}{\sum }}{\left(-1\right)}^{p+1}{\displaystyle \underset{j=0}{\overset{p}{\prod }}}\frac{M_{l_j,l_{j+1}}}{M_{l_j,l_j}-{\lambda }_{l_{p+1}}}$ (2.19)

For simplicity, we set the scaling factor $v_{l_{p+1}}^{l_{p+1}}=1$ for all $l_{p+1}$. It can be seen that $l_0>l_{p+1}$ would contradict with $l_0<l_1<\cdots <l_p<l_{p+1}$, hence (2.19) vanishes and the identity (1) holds (for the eigenvectors of the truncated matrix). Since we know that $\forall b>a$ $v_a^b=0$, the eigenvectors of the truncated matrix (2.18) with zeros substituted for the rest of the components are also the solution to the infinite Equation (2.17). Therefore, expressions (1), (2), (3) and (2.16) also hold for the matrix M.

The following formula is an explicit expression for the inverse of an upper triangular matrix:

Lemma 2.5. Let U be an infinite upper triangular matrix with complex entries:

$U_{km}=\left[\begin{array}{ccccc}{U}_{00}& U_{01}& U_{02}& U_{03}& \cdots \\ 0& U_{11}& U_{12}& U_{13}& \cdots \\ 0& 0& U_{22}& U_{23}& \cdots \\ 0& 0& 0& U_{33}& \cdots \\ ⋮& ⋮& ⋮& ⋮& \ddots \end{array}\right]$ (2.20)

and let $\forall n$ $U_{nn}\ne 0$. The explicit expression for the inverse of U is given by:

$U_{km}^{-1}=(\begin{array}{l}0\text{for}k>m\\ \frac{1}{U_{km}}\text{for}k=m\end{array}$ (2.21)

and:

$U_{l_0,l_{p+1}}^{-1}=\frac{1}{U_{l_{p+1},l_{p+1}}}{\displaystyle \underset{p=0}{\overset{l_{p+1}-l_0-1}{\sum }}}{\displaystyle \underset{l_0<l_1<\cdots <l_p<l_{p+1}}{\sum }}{\left(-1\right)}^{p+1}{\displaystyle \underset{j=0}{\overset{p}{\prod }}}\frac{U_{l_j,l_{j+1}}}{U_{l_j,l_j}}$ (2.22)

for $l_0<l_{p+1}$.

Proof. The inverse of the upper triangular matrix is also upper triangular [25], therefore the equation $U{U}^{-1}=I$ reads:

$\left[\begin{array}{ccccc}{U}_{00}& U_{01}& U_{02}& U_{03}& \cdots \\ 0& U_{11}& U_{12}& U_{13}& \cdots \\ 0& 0& U_{22}& U_{23}& \cdots \\ 0& 0& 0& U_{33}& \cdots \\ ⋮& ⋮& ⋮& ⋮& \ddots \end{array}\right]\left[\begin{array}{ccccc}{U}_{00}^{-1}& U_{01}^{-1}& U_{02}^{-1}& U_{03}^{-1}& \cdots \\ 0& U_{11}^{-1}& U_{12}^{-1}& U_{13}^{-1}& \cdots \\ 0& 0& U_{22}^{-1}& U_{23}^{-1}& \cdots \\ 0& 0& 0& U_{33}^{-1}& \cdots \\ ⋮& ⋮& ⋮& ⋮& \ddots \end{array}\right]=I$ (2.23)

For the diagonal elements, we have:

$U_{nn}^{-1}=\frac{1}{U_{nn}}$ (2.24)

Elements further off the diagonal can be obtained by the back substitution algorithm. Using (2.23):

$U_{nn}{U}_{n,n+1}^{-1}+U_{n,n+1}{U}_{n+1,n+1}^{-1}=0$ (2.25)

Which together with (2.24) gives:

$U_{n,n+1}^{-1}=-\frac{U_{n,n+1}{U}_{n+1,n+1}^{-1}}{U_{nn}}=-\frac{U_{n,n+1}}{U_{nn}{U}_{n+1,n+1}}$ (2.26)

By repeating the same process for the rest of the elements of $U_{km}^{-1}$, the relation (2.22) for $k<m$ is obtained.

The powers of the transition matrix can be calculated by diagonalizing T:

$T^n={\left(PD{P}^{-1}\right)}^n=P{D}^n{P}^{-1}$ (2.27)

where D denotes the infinite matrix with the eigenvalues of T as diagonal entries, while P is the infinite modal matrix composed of eigenvectors of T:

$D=diag\left[{\lambda }_0,{\lambda }_1,{\lambda }_2,\cdots \right],P=\left[v_0,v_1,v_2,\cdots \right]$ (2.28)

We now prove the main theorem of this section:

Theorem 2.6. Let $u_i\in ℂ\mathrm{,}i\ge 0$ be a sequence of numbers satisfying the recurrence relation (2.1). Let $F\left(u_{i-1}\right)={\displaystyle {\sum }_{j=0}^m}{c}_j{\left(u_{i-1}\right)}^j$ be polynomial in $u_{i-1}$ such that the corresponding polynomial equation:

${\displaystyle {\sum }_{j=0}^m{c}_j{d}^j}=d$ (2.29)

has a root $d\in ℂ$ satisfying:

$\left|{\displaystyle \underset{j=1}{\overset{m}{\sum }}}j{c}_j{d}^{j-1}\right|\ne 0\wedge \left(\left|{\displaystyle \underset{j=1}{\overset{m}{\sum }}}j{c}_j{d}^{j-1}\right|\ne 1\vee arg\left({\displaystyle \underset{j=1}{\overset{m}{\sum }}}j{c}_j{d}^{j-1}\right)\notin \mathbb{Q}\right)$ (2.30)

Then, the solution to the recurrence (2.1) is given by:

$u_i=d+{\displaystyle {\sum }_{l=1}^{\infty }{\displaystyle {\sum }_{j=1}^l{v^{\prime }}_j^1{c^{\prime }}_1^{ij}{P^{\prime }}_{jl}^{-1}{\left(u_0-d\right)}^l}}$ (2.31)

where:

${P^{\prime }}_{l_0,l_{p+1}}^{-1}=(\begin{array}{l}0\text{for}{l}_0>l_{p+1}\\ 1\text{for}{l}_0=l_{p+1}\\ {\displaystyle {\sum }_{p=0}^{l_{p+1}-l_0-1}{\displaystyle {\sum }_{l_0<l_1<\cdots <l_p<l_{p+1}}{\left(-1\right)}^{p+1}}}{\displaystyle {\prod }_{j=0}^p\frac{{v^{\prime }}_{l_{j+1}}^{l_j}}{{v^{\prime }}_{l_j}^{l_j}}}\text{for}{l}_0<l_{p+1}\end{array}$

$\begin{array}{l}{v^{\prime }}_{l_{p+1}}^{l_0}=(\begin{array}{l}0\text{for}{l}_0>l_{p+1}\\ 1\text{for}{l}_0=l_{p+1}\\ {\displaystyle {\sum }_{p=0}^{l_{p+1}-l_0-1}{\displaystyle {\sum }_{l_0<l_1<\cdots <l_p<l_{p+1}}{\left(-1\right)}^{p+1}}}{\displaystyle {\prod }_{j=0}^p\frac{{T^{\prime }}_{l_j,l_{j+1}}}{{T^{\prime }}_{l_j,l_j}-{\lambda }_{l_{p+1}}}}\text{for}{l}_0<l_{p+1}\end{array}\\ {T^{\prime }}_{ab}={\displaystyle \underset{{\displaystyle {\sum }_{l=0}^m{k}_l}=a,{\displaystyle {\sum }_{l=0}^{m}l{k}_l}=b}{\sum }}\frac{a!}{{\displaystyle {\prod }_{s=0}^m{k}_s!}}{\displaystyle \underset{t=0}{\overset{m}{\prod }}}{c^{\prime }}_t^{k_t}\\ {c^{\prime }}_0=0,{c^{\prime }}_l={\displaystyle \underset{j=l}{\overset{m}{\sum }}}{c}_j\left(\begin{array}{c}j\\ l\end{array}\right)d^{j-l}\text{for}l\ne 0\end{array}$ (2.32)

Proof. Because the Equation (2.29) is equivalent to requiring ${c^{\prime }}_0=0$ ( ${c^{\prime }}_0$ is the free term of the polynomial $F^{\prime }$ obtained after shift of variables by d), the variables ${u^{\prime }}_i=u_i-d$ are the ones in which the free term vanishes. In the same way, from (2.9) we have:

${c^{\prime }}_1={\displaystyle \underset{j=1}{\overset{m}{\sum }}}{c}_j\left(\begin{array}{c}j\\ 1\end{array}\right)d^{j-1}={\displaystyle \underset{j=1}{\overset{m}{\sum }}}j{c}_j{d}^{j-1}$ (2.33)

Therefore condition (2.30) expresses the fact that all of the powers of ${c^{\prime }}_1$ have different values.

The transition matrix T is (2.3):

${T^{\prime }}_{ab}={\displaystyle \underset{{\displaystyle {\sum }_{l=0}^m{k}_l}=a,{\displaystyle {\sum }_{l=0}^{m}l{k}_l}=b}{\sum }}\frac{a!}{{\displaystyle {\prod }_{s=0}^m{k}_s!}}{\displaystyle \underset{t=0}{\overset{m}{\prod }}}{c^{\prime }}_t^{k_t}$ (2.34)

By Lemma 2.2, ${c^{\prime }}_0=0$ and the infinite transition matrix is upper triangular, while (2.30) implies that T is diagonalizable. Furthermore, by Lemma 2.4 the a-th eigenvector has at most $a+1$ non-zero elements, and the (infinite) matrix P can be chosen to be upper triangular (the inverse of P is then upper triangular as well, Lemma 2.5). The powers of transition matrix T can be obtained from (2.27). Hence, using (1.8), (2.34) and Lemma 2.4:

$\begin{array}{l}{u^{\prime }}_i=e_1{T^{\prime }}^i{y^{\prime }}^0=e_{1}P{D^{\prime }}^i{P}^{-1}{y^{\prime }}^0\\ ={\left[\begin{array}{c}0\\ 1\\ 0\\ ⋮\end{array}\right]}^{\text{T}}\left[\begin{array}{cccc}{v^{\prime }}_0^0& {v^{\prime }}_1^0& {v^{\prime }}_2^0& \cdots \\ 0& {v^{\prime }}_1^1& {v^{\prime }}_2^1& \cdots \\ 0& 0& {v^{\prime }}_2^2& \cdots \\ ⋮& ⋮& ⋮& \ddots \end{array}\right]{\left[\begin{array}{cccc}1& 0& 0& \cdots \\ 0& {c^{\prime }}_1& 0& \cdots \\ 0& 0& {c^{\prime }}_1^2& \cdots \\ ⋮& ⋮& ⋮& \ddots \end{array}\right]}^i\left[\begin{array}{cccc}{P^{\prime }}_{00}^{-1}& {P^{\prime }}_{01}^{-1}& {P^{\prime }}_{02}^{-1}& \cdots \\ 0& {P^{\prime }}_{11}^{-1}& {P^{\prime }}_{12}^{-1}& \cdots \\ 0& 0& {P^{\prime }}_{22}^{-1}& \cdots \\ ⋮& ⋮& ⋮& \ddots \end{array}\right]\left[\begin{array}{c}1\\ {u^{\prime }}_0\\ {u^{\prime }}_0^2\\ ⋮\end{array}\right]\end{array}$ (2.35)

where $P_{ab}^{-1}$ denote the elements of the inverse of P that can be calculated using Lemma 2.5, while the formula for ${v^{\prime }}_l^j$ is given in (2.16). The i-th power of the diagonal matrix is just a matrix with i-th power of its elements on the diagonal. Carrying out the matrix multiplication, we obtain:

${u^{\prime }}_i={\displaystyle \underset{l=1}{\overset{\infty }{\sum }}}{u^{\prime }}_0^l{\displaystyle \underset{j=1}{\overset{l}{\sum }}}{v^{\prime }}_j^1{\left({c^{\prime }}_1^j\right)}^i{P^{\prime }}_{jl}^{-1}={\displaystyle \underset{l=1}{\overset{\infty }{\sum }}}{\displaystyle \underset{j=1}{\overset{l}{\sum }}}{v^{\prime }}_j^1{c^{\prime }}_1^{ij}{P^{\prime }}_{jl}^{-1}{u^{\prime }}_0^l$ (2.36)

By retrieving the original sequence $u_i={u^{\prime }}_i+d$ and substituting back the original coefficients $c_j$ for the shifted parameters ${c^{\prime }}_j$, we arrive at the result.

Theorem 2.6 allows for solving a large family of nonlinear recurrences. It is apparent from (2.31) that the general solution has a qualitative form:

$u_i={\displaystyle \underset{j=0}{\overset{\infty }{\sum }}}{f}_j\left(i\right)u_0^j$ (2.37)

where $f_j$ are some (in general complex) functions of number of iterations i. Therefore, Theorem 2.6 can be used to determine the form of functions $f_j$. As an example, the results are applied to the well-known logistic mapping:

Example 2.7. The logistic map is a recurrence relation of the form [26]:

$u_i=r{u}_{i-1}-r{u}_{i-1}^2$ (2.38)

with the coefficient $r\in ℝ$. Conditions (2.29) and (2.30) are trivially fulfilled for $r\ne 1$ and $d=0$. We obtain:

$\begin{array}{l}P=\left[\begin{array}{ccccc}1& 0& 0& 0& \cdots \\ 0& 1& \frac{1}{1-r}& \frac{2}{r^3-r^2-r+1}& \cdots \\ 0& 0& 1& \frac{2}{1-r}& \cdots \\ 0& 0& 0& 1& \cdots \\ ⋮& ⋮& ⋮& ⋮& \ddots \end{array}\right],\\ P^{-1}=\left[\begin{array}{ccccc}1& 0& 0& 0& \cdots \\ 0& 1& \frac{1}{r-1}& \frac{2r}{r^3-r^2-r+1}& \cdots \\ 0& 0& 1& \frac{2}{r-1}& \cdots \\ 0& 0& 0& 1& \cdots \\ ⋮& ⋮& ⋮& ⋮& \ddots \end{array}\right]\end{array}$ (2.39)

Thus, by Theorem 2.6 the solution takes the form:

$u_i=r^i{u}_0+\frac{r^i-r^{2i}}{r-1}{u}_0^2+\frac{2r{r}^i-2\left(r+1\right)r^{2i}+2r^{3i}}{r^3-r^2-r+1}{u}_0^3+\cdots$ (2.40)

It can be seen that the “coefficient functions” $f_j$ in front of the powers of initial value $u_0$ become progressively more complicated as the order of $u_0^j$ increases.

3. Multi-Variable Recurrences of Depth One

The treatment of polynomial multi-variable recurrences of depth one follows the same structure as Section 2. With appropriately chosen auxiliary vectors, the recurrence can be Carleman-linearized and represented by an appropriate transition matrix. Let:

$z^i={\left[u_i^1,u_i^2,\cdots ,u_i^k\right]}^{\text{T}}$ (3.1)

be a vector composed of k variables $u_i^k$. Then, a general polynomial recurrence of depth one in k variables can be compactly written in the form:

$\begin{array}{c}{u}_i^{p+1}=z_p^i=F_{p+1}\left(u_{i-1}^1,u_{i-1}^2,\cdots ,u_{i-1}^k\right)={}_p{}^{0}C+{\displaystyle \underset{j=1}{\overset{m}{\sum }}}{\displaystyle \underset{l_1,\cdots ,l_j=0}{\overset{k-1}{\sum }}}{}_p{}^{j}C{}_{l_1,\cdots ,l_j}{\displaystyle \underset{s=1}{\overset{j}{\prod }}}{z}_{l_s}^{i-1}\\ ={}_p{}^{0}C+{\displaystyle \underset{l_1=0}{\overset{k-1}{\sum }}}{}_p{}^{1}C{}_{l_1}{z}_{l_1}^{i-1}+{\displaystyle \underset{l_1=0}{\overset{k-1}{\sum }}}{\displaystyle \underset{l_2=0}{\overset{k-1}{\sum }}}{}_p{}^{2}C{}_{l_1,l_2}{z}_{l_1}^{i-1}{z}_{l_2}^{i-1}+\cdots \end{array}$ (3.2)

where symbols ${}_p{}^{j}C{}_{l_1,\cdots ,l_j}$ represent j + 1 dimensional arrays. It is apparent that the description in terms of ${}_p{}^{j}C{}_{l_1,\cdots ,l_j}$ is somewhat redundant: every permutation of lower right indices $l_1,\cdots ,l_j$ is the coefficient of the same combination of powers of variables $u_i^k$. For convenience, we choose ${}_p{}^{j}C{}_{l_1,\cdots ,l_j}$ to vanish unless the lower right indices are in the non-decreasing order, i.e. ${}_p{}^{j}C{}_{l_1,\cdots ,l_j}=0$ if $\neg \left(l_1\le \cdots \le l_j\right)$. As a result, the symbol ${}_p{}^{j}C{}_{l_1,\cdots ,l_j}$ has only $k\left(\begin{array}{c}k+j-1\\ j\end{array}\right)$ independent components, when all of the lower indices are taken into account [27].

Theorem 3.1. Let the polynomial multi-variable depth one recurrence be defined as in (3.2), and let the Kronecker powers of $z^i$ be defined as in [13]:

${\left(z^i\right)}^{\alpha }={\left[u_i^1{\left(z^i\right)}^{\alpha -1}\mathrm{,}{u}_i^2{\left(z^i\right)}^{\alpha -1}\mathrm{,}\cdots \mathrm{,}{u}_i^k{\left(z^i\right)}^{\alpha -1}\right]}^{\text{T}}\mathrm{,}{\left(z^i\right)}^0\mathrm{<mo>\; :\; </mo>}=1$ (3.3)

Then, the Carleman linearization of the system can be defined on infinite auxiliary vectors:

$y^i={\left[{\left(z^i\right)}^0\mathrm{,}{\left(z^i\right)}^1\mathrm{,}{\left(z^i\right)}^2\mathrm{,}\cdots \right]}^{\text{T}}$ (3.4)

and the solution to recurrence (3.2) can be written as:

$u_i^k=e_k{T}^i{y}^0$ (3.5)

where $y^0$ is defined by initial conditions and the infinite dimensional transition matrix is:

$T_{ab}=\left[\begin{array}{cccc}1& 0& 0& \cdots \\ {}_p{}^{0}C& {}_p{}^{1}C{}_{l_1}& {}_p{}^{2}C{}_{\mathrm{1,}{l}_1}& \cdots \\ {}_p{}^{0}C\otimes {}_p{}^{0}C& {}_p{}^{0}C\otimes {}_p{}^{1}C{}_{l_1}+{}_p{}^{1}C{}_{l_1}\otimes {}_p{}^{0}C& \cdots & \cdots \\ {}_p{}^{0}C\otimes {}_p{}^{0}C\otimes {}_p{}^{0}C& ⋮& \ddots & \ddots \\ ⋮& ⋮& \ddots & \ddots \end{array}\right]$ (3.6)

The Kronecker product in (3.6) acts on the subscripts only.

Proof. Since $y^i$ is composed of consecutive Kronecker powers of $z^i$, any polynomial $F\left(u_i^1\mathrm{,}{u}_i^2\mathrm{,}\cdots \mathrm{,}{u}_i^k\right)$ in k variables can be rewritten as linear combination of elements of $y^i$. The recurrence is defined solely by the arrays ${}_p{}^{j}C{}_{l_1,\cdots ,l_j}$, therefore the transition matrix consists of the numbers appearing in ${}_p{}^{j}C{}_{l_1,\cdots ,l_j}$. The elements of vector $y^i$ can be identified with the product of elements of $z^i$. From the definitions (3.3), (3.4) and the properties of the Kronecker product we have:

$\begin{array}{l}{y}_l^i=(\begin{array}{l}1\text{for}l=0\\ {\displaystyle {\prod }_{s=1}^{\lfloor {\log }_k\left(1+l\left(k-1\right)\right)\rfloor }{u}_i^{\left(\lfloor \left(l-\frac{1-k^s}{1-k}\right)k^{1-s}\rfloor \%k\right)+1}}\end{array}\\ {\displaystyle \underset{j=1}{\overset{s}{\prod }}}{u}_i^{l_j}=y_{\frac{1-k^s}{1-k}+{\displaystyle {\sum }_{j=1}^s{k}^{s-j}\left(l_j-1\right)}}^i=y_{{\displaystyle {\sum }_{j=1}^s{k}^{s-j}{l}_j}}^i\end{array}$ (3.7)

The product of elements of $z^i$ can be expressed in terms of elements of $z^{i-1}$ :

${\displaystyle \underset{j=0}{\overset{k-1}{\prod }}}{\left(z_j^i\right)}^{t_j}={\displaystyle \underset{j=0}{\overset{k-1}{\prod }}}{\left(F_j\left(z^{i-1}\right)\right)}^{t_j}$ (3.8)

Analogously to (1.9), the infinite transition matrix can be defined as:

$\begin{array}{l}{T}_{ab}=R_b\left({\displaystyle \underset{s=1}{\overset{\lfloor {\log }_k\left(1+b\left(k-1\right)\right)\rfloor }{\prod }}}\frac{\partial }{\partial u_{i-1}^{\left(\lfloor \left(b-\frac{1-k^s}{1-k}\right)k^{1-s}\rfloor \%k\right)+1}}\right)\\ \cdot {\left({\displaystyle \underset{s=1}{\overset{\lfloor {\log }_k\left(1+a\left(k-1\right)\right)\rfloor }{\prod }}}{F}_{\left(\lfloor \left(a-\frac{1-k^s}{1-k}\right)k^{1-s}\rfloor \%k\right)+1}\left(u_{i-1}^1,\cdots ,u_{i-1}^k\right)\right)|}_{u_{i-1}^k=0}\end{array}$ (3.9)

The constants $R_b$ do not depend on $F_i$, and are defined as a product:

$R_b=\frac{1}{c_b}{\displaystyle \underset{s=1}{\overset{k}{\prod }}}\frac{1}{a_s!}$ (3.10)

where $a_s$ is the number of times $u_{i-1}^s$ appears in $y_b^{i-1}$, while $c_b$ is the number of times the element $y_b^i$ appears in $y^i$. Equivalently, $T_{ab}$ can be written in the matrix form:

$\left[\begin{array}{ccccc}1& 0& \cdots & 0& \cdots \\ {}_0{}^{0}C& {}_0{}^{1}C{}_0& \cdots & {}_0{}^{1}C{}_{k-1}& \cdots \\ ⋮& ⋮& \ddots & ⋮& \cdots \\ {}_{k-1}{}^{0}C& {}_{k-1}{}^{1}C{}_0& \cdots & {}_{k-1}{}^{1}C{}_{k-1}& \cdots \\ {}_0{}^{0}C{}_0{}^{0}C& {}_0{}^{0}C{}_0{}^{1}C{}_0+{}_0{}^{0}C{}_0{}^{1}C{}_0& \cdots & {}_0{}^{0}C{}_0{}^{1}C{}_{k-1}+{}_0{}^{0}C{}_0{}^{1}C{}_{k-1}& \cdots \\ ⋮& ⋮& \ddots & ⋮& \cdots \\ {}_0{}^{0}C{}_{k-1}{}^{0}C& {}_0{}^{0}C{}_{k-1}{}^{1}C{}_0+{}_{k-1}{}^{0}C{}_0{}^{1}C{}_0& \cdots & {}_0{}^{0}C{}_{k-1}{}^{1}C{}_{k-1}+{}_{k-1}{}^{0}C{}_0{}^{1}C{}_{k-1}& \cdots \\ {}_1{}^{0}C{}_0{}^{0}C& {}_1{}^{0}C{}_0{}^{1}C{}_0+{}_0{}^{0}C{}_1{}^{1}C{}_0& \cdots & {}_1{}^{0}C{}_0{}^{1}C{}_{k-1}+{}_0{}^{0}C{}_1{}^{1}C{}_{k-1}& \cdots \\ ⋮& ⋮& \ddots & ⋮& \cdots \\ {}_1{}^{0}C{}_{k-1}{}^{0}C& {}_1{}^{0}C{}_{k-1}{}^{1}C{}_0+{}_{k-1}{}^{0}C{}_1{}^{1}C{}_0& \cdots & {}_1{}^{0}C{}_{k-1}{}^{1}C{}_{k-1}+{}_{k-1}{}^{0}C{}_1{}^{1}C{}_{k-1}& \cdots \\ ⋮& ⋮& \ddots & ⋮& \cdots \\ {}_{k-1}{}^{0}C{}_{k-1}{}^{0}C& {}_{k-1}{}^{0}C{}_{k-1}{}^{1}C{}_0+{}_{k-1}{}^{0}C{}_{k-1}{}^{1}C{}_0& \cdots & {}_{k-1}{}^{0}C{}_{k-1}{}^{1}C{}_{k-1}+{}_{k-1}{}^{0}C{}_{k-1}{}^{1}C{}_{k-1}& \cdots \\ ⋮& ⋮& ⋮& ⋮& \ddots \end{array}\right]$ (3.11)

By replacing the groups of elements from (3.11) with the block-matrix notation of Kronecker powers of arrays ${}_p{}^{j}C{}_{l_1,\cdots ,l_j}$, we arrive at (3.6).

For $k=1$, expression (3.6) reduces to (2.2). Furthermore, the analogue of the Theorem 2.1 can be proven. If ${}_p{}^{0}C=0$ for all p, then all of the Kronecker products involving ${}_p{}^{0}C$ vanish. Since all of the block matrices below the diagonal in (3.6) include the ${}_p{}^{0}C$ in the Kronecker product, the lower part of the block transition matrix vanishes and the matrix becomes block upper triangular. The block matrices on the diagonal of $T_{ab}$ are the Kronecker powers of the array standing near the linear term ${}_p{}^{1}C{}_l$. The Kronecker powers of ${}_p{}^{1}C{}_l$ are defined as in [13]:

${\left({}^{1}C\right)}^{\alpha }=\left[\begin{array}{cccc}{}_0{}^{1}C{}_0{\left({}^{1}C\right)}^{\alpha -1}& {}_0{}^{1}C{}_1{\left({}^{1}C\right)}^{\alpha -1}& \cdots & {}_0{}^{1}C{}_{k-1}{\left({}^{1}C\right)}^{\alpha -1}\\ {}_1{}^{1}C{}_0{\left({}^{1}C\right)}^{\alpha -1}& {}_1{}^{1}C{}_1{\left({}^{1}C\right)}^{\alpha -1}& \cdots & {}_1{}^{1}C{}_{k-1}{\left({}^{1}C\right)}^{\alpha -1}\\ ⋮& ⋮& \ddots & ⋮\\ {}_{k-1}{}^{1}C{}_0{\left({}^{1}C\right)}^{\alpha -1}& {}_{k-1}{}^{1}C{}_1{\left({}^{1}C\right)}^{\alpha -1}& \cdots & {}_{k-1}{}^{1}C{}_{k-1}{\left({}^{1}C\right)}^{\alpha -1}\end{array}\right],{\left({}^{1}C\right)}^0=1$ (3.12)

Therefore, the transition matrix can be brought to the upper triangular form if $\forall p{}_p{}^{0}C=0\wedge \forall p>l{}_p{}^{1}C{}_l=0$. Even though in general recurrences do not fulfill those conditions, variables can be (in some cases) transformed in such a way that the new transition matrix is upper triangular.

Lemma 3.2. Let the recurrence on $u_i^k$ be defined as in (3.2). There exists a linear transformation of variables for $F_k$ such that the transition matrix $T_{ab}$ corresponding to the new sequence is upper triangular if and only if there exists k dimensional vector B satisfying:

$B_p={\displaystyle \underset{s=0}{\overset{m}{\sum }}}\frac{1}{s!}{\displaystyle \underset{l_1,\cdots ,l_s=0}{\overset{k-1}{\sum }}}{\displaystyle \underset{\sigma }{\sum }}{}_p{}^{s}C{}_{l_{{\sigma }_1},\cdots ,l_{{\sigma }_s}}{\displaystyle \underset{t=1}{\overset{s}{\prod }}}{B}_{l_t}$ (3.13)

Proof. The general (invertible) linear transformation of k variables $u_i^k$ can be written in the matrix form:

$z^i=A{z^{\prime }}^i+B\iff {z^{\prime }}^i=A^{-1}\left(z^i-B\right)$ (3.14)

with $det\left(A\right)\ne 0$. The new recurrence defined on primed variables ${z^{\prime }}^i$ satisfies:

${z^{\prime }}_{p-1}^i={}_p{}^0{C}^{\prime }+{\displaystyle \underset{j=1}{\overset{m}{\sum }}}{\displaystyle \underset{l_1,\cdots ,l_j=0}{\overset{k-1}{\sum }}}{}_p{}^j{C}^{\prime }{}_{l_1,\cdots ,l_j}{\displaystyle \underset{s=1}{\overset{j}{\prod }}}{z^{\prime }}_{l_s}^{i-1}$ (3.15)

where:

$\begin{array}{l}{}_k{}^0{C}^{\prime }=-A_{kj}^{-1}{B}_j+{\displaystyle \underset{s=0}{\overset{m}{\sum }}}{\displaystyle \underset{j,l_1,\cdots ,l_s=0}{\overset{k-1}{\sum }}}{A}_{kj}^{-1}{}_j{}^s\hat{C}{}_{l_1,\cdots ,l_s}{\displaystyle \underset{t=1}{\overset{s}{\prod }}}{B}_{l_t}\\ {}_k{}^p{C}^{\prime }{}_{x_1,\cdots ,x_p}={\displaystyle \underset{s=p}{\overset{m}{\sum }}}\left(\begin{array}{c}s\\ p\end{array}\right){\displaystyle \underset{j,l_1,\cdots ,l_s=0}{\overset{k-1}{\sum }}}{A}_{kj}^{-1}{}_j{}^s\hat{C}{}_{l_1,\cdots ,l_s}\left({\displaystyle \underset{t_1=1}{\overset{p}{\prod }}}{A}_{l_{t_1}{x}_{t_1}}\right)\left({\displaystyle \underset{t_2=p+1}{\overset{s}{\prod }}}{B}_{l_{t_2}}\right)\\ {}_j{}^s\hat{C}{}_{l_1,\cdots ,l_s}=\frac{1}{s!}{\displaystyle \underset{\sigma }{\sum }}{}_j{}^{s}C{}_{l_{{\sigma }_1},\cdots ,l_{{\sigma }_s}}\end{array}$ (3.16)

The symbol ${}_j{}^s\hat{C}{}_{l_1,\cdots ,l_s}$ denotes the symmetrized array that is equal to ${}_j{}^{s}C{}_{l_1,\cdots ,l_s}$ averaged over all of the permutations $\sigma$ of the right subscripts.

It is trivial that convolution of linear transformations is a linear transformation. We firstly perform a shift of the variables by the vector B. From (3.16) we have:

${}_p{}^0{C}^{\prime }=-B_k+{\displaystyle \underset{s=0}{\overset{m}{\sum }}}\frac{1}{s!}{\displaystyle \underset{l_1,\cdots ,l_s=0}{\overset{k-1}{\sum }}}{\displaystyle \underset{\sigma }{\sum }}{}_p{}^{s}C{}_{l_{{\sigma }_1},\cdots ,l_{{\sigma }_s}}{\displaystyle \underset{t=1}{\overset{s}{\prod }}}{B}_{l_t}$ (3.17)

Setting ${}_p{}^0{C}^{\prime }=0$, above equation reduces to a set of coupled polynomial equations. Unfortunately, the resulting system of equations is not always solvable. Assuming that (3.17) has a solution, the other arrays transform according to (3.16).

Subsequently a linear transformation is performed with $B=0$. The constant term ${}_p{}^0{C}^{″}={}_p{}^0{C}^{\prime }=0$ (i.e. still vanishes, regardless of the choice of A), while the linear term is:

${}_p{}^1{C}^{″}{}_t={\displaystyle \underset{j,l_1=0}{\overset{k-1}{\sum }}}{A}_{pj}^{-1}{}_j{}^1{C}^{\prime }{}_l{A}_{lt}$ (3.18)

Similarly, there always exists a suitable choice of $k\times k$ dimensional matrix A such that ${}_p{}^{1}C{}_t$ becomes upper triangular in lower indices [13] and thus the transition matrix can be brought to an upper triangular form provided (3.13) has a solution.

If the shift vector B satisfying (3.13) exists, it’s not necessary unique. The situation is analogous to the one found for a uni-variable recurrence. Different choices of B satisfying (3.13) with ${}_p{}^0{C}^{\prime }=0$ will lead to ${}_p{}^{1}C{}_l$ with different eigenvalues. Therefore, it is suitable to choose B such that it satisfies the diagonalizability conditions stated in Lemma 3.3.

It is worth noting that even though the arrays ${}_p{}^{j}C{}_{l_1,\cdots ,l_j}$ vanish for a non-decreasing order of lower right indices, the same may not be true for the transformed arrays ${}_p{}^j{C}^{\prime }{}_{l_1,\cdots ,l_j}$. Therefore, after transforming the coefficient arrays with the lineary transformation (3.13), it is necessary to upper triangularize the transformed coefficients so that they also satisfy ${}_p{}^j{C}^{\prime }{}_{l_1,\cdots ,l_j}=0$ if $\neg \left(l_1\le \cdots \le l_j\right)$.

Since the eigenvalues of the Kronecker product of two matrices are the products of eigenvalues of the two matrices, we have:

$Spec\left({\left({}^{1}C\right)}^{\alpha }\right)=Spec{\left({}^{1}C\right)}^{\alpha }$ (3.19)

where the right-hand side represents the Cartesian powers of the set of eigenvalues of ${}^{1}C$ [24]. As a result, the general eigenvalue is of the form:

$\lambda ={{\lambda }^{\prime }}_0^{a_0}{{\lambda }^{\prime }}_1^{a_1}\cdots {{\lambda }^{\prime }}_{k-1}^{a_{k-1}}$ (3.20)

for some numbers $a_p\in \mathbb{N}$, ${{\lambda }^{\prime }}_p\in Spec\left({}^{1}C\right)$. In the same way as in (3.7), the l-th eigenvalue of the transition can be written as:

${\lambda }_l=(\begin{array}{l}1\text{for}l=0\\ {\displaystyle {\prod }_{s=1}^{\lfloor {\log }_k\left(1+l\left(k-1\right)\right)\rfloor }{{\lambda }^{\prime }}_{\lfloor \left(l-\frac{1-k^s}{1-k}\right)k^{1-s}\rfloor \%k}}\end{array}$ (3.21)

The matrix is diagonalizable if all of its eigenvalues are different, thus a sufficient condition for the existence of diagonal form can be constructed.

Although this is true, a problem arises due to redundancy in description in terms of the auxiliary vectors (3.4). Since $y^i$ are composed of Kronecker powers of $z^i$, there are components that are redundant, i.e. that represent the same product of variables $u_i^k$. For example, if $k=2$, the infinite auxiliary vectors are:

$y^i={\left[1,u_i^1,u_i^2,{\left(u_i^1\right)}^2,u_i^1{u}_i^2,u_i^2{u}_i^1,{\left(u_i^2\right)}^2,\cdots \right]}^{\text{T}}$ (3.22)

and it is apparent that the 4-th and 5-th components are the same.

Due to the redundancy, the repeating copies of the same eigenvalues appear in the upper triangular form of the transition matrix (3.6). Lemmas 2.4 and 2.5 are applicable only if the eigenvalues of $T_{ab}$ are non-degenerate, hence the Lemmas cannot be used to diagonalize the transition matrix. Nevertheless, the special choice of ${}_p{}^{j}C{}_{l_1,\cdots ,l_j}$ allows one to “reduce” the transition matrix and get rid of the undesired copies.

The redundant terms arise from the permutation of lower right indices of ${}_p{}^j{C}^{\prime }{}_{l_1,\cdots ,l_j}$. By definition (the begining of Section 3) the arrays vanish for a decreasing sequence of lower indices. All of the columns in $T_{ab}$ associated with the copies of the same power of variables $u_i^k$ vanish except for the diagonal element. For example, since the 4-th and 5-th element of the auxiliary vectors (3.22) are similar, all of the elements $T_{a5}=0$, $a\ne 5$ vanish. This observation allows one to get rid of the redundant columns (and rows) of $T_{ab}$, hence making the calculation of the powers of the transition matrix feasible. In other words, the redundancy has to be taken into account when assessing the diagonalizability of $T_{ab}$.

In order to differentiate between quantities related to the redundant (original) and non-redundant system, tilde notation is introduced. The quantities related to the system without redundancy (vectors without redundant components and matrices without redundant rows and columns) are denoted with tilde above the symbol. For example, the equivalent of the auxiliary vector (3.22) for the non-redundant Carleman embedding is:

${\tilde{y}}^i={\left[\mathrm{1,}{u}_i^1\mathrm{,}{u}_i^2\mathrm{,}{\left(u_i^1\right)}^2\mathrm{,}{u}_i^1{u}_i^2\mathrm{,}{\left(u_i^2\right)}^2\mathrm{,}\cdots \right]}^{\text{T}}$ (3.23)

Similarly, the transition matrix for ${\tilde{y}}^i$ is denoted by ${\tilde{T}}_{ab}$ respectively. Although it is the non-redundant Carleman embedding defined by ${\tilde{T}}_{ab}$ and ${\tilde{y}}^i$ that is eventually diagonalized, it is hard to express it using a concise formula. Equivalents of expression (3.4) for ${\tilde{y}}^i$ and (3.6) for ${\tilde{T}}_{ab}$ have not been found. As opposed to the redundant system, the transition matrix ${\tilde{T}}_{ab}$ cannot be simply expressed in terms of Kronecker products. On the other hand, it is much easier to calculate powers of ${\tilde{T}}_{ab}$ rather than $T_{ab}$. The following can be said about diagonalizability of the “non-redundant” transition matrix:

Lemma 3.3. A sufficient condition for the infinite upper triangular transition matrix ${\tilde{T}}_{ab}$ to be diagonalizable is for the corresponding $k\times k$ dimensional matrix ${}^{1}C$ to have eigenvalues ${{\lambda }^{\prime }}_p$ such that the equation:

$(\begin{array}{l}{x}_0+x_1{\log }_{\left|{\lambda }_0\right|}\left|{\lambda }_1\right|+\cdots +x_{k-1}{\log }_{\left|{\lambda }_0\right|}\left|{\lambda }_{k-1}\right|=0\\ Arg\left({\lambda }_0\right)x_0+Arg\left({\lambda }_1\right)x_1+\cdots +Arg\left({\lambda }_{k-1}\right)x_{k-1}=2\pi t\end{array}$ (3.24)

has no (non-trivial) solution $x_0\mathrm{,}\cdots \mathrm{,}{x}_{k-1}$, $t\in \mathbb{Z}$. In case of two variables $k=2$, this can be further simplified to:

${\log }_{\left|{\lambda }_0\right|}\left|{\lambda }_1\right|\notin \mathbb{Q}\vee \frac{\pi }{Arg\left({\lambda }_1\right)-Arg\left({\lambda }_0\right){\log }_{\left|{\lambda }_0\right|}\left|{\lambda }_1\right|}\notin \mathbb{Q}$ (3.25)

Proof. As stated before, the redundancy in the auxiliary vectors $y^i$ implies that there always will be an eigenvector with algebraic multiplicity of at least two. However, since all of the elements of the corresponding column except the diagonal vanish, a sufficient condition for diagonalization can be still formulated based on the tilded transition matrix ${\tilde{T}}_{ab}$.

The condition for ${\tilde{T}}_{ab}$ to have different eigenvalues can be written as:

${{\lambda }^{\prime }}_0^{a_0}{{\lambda }^{\prime }}_1^{a_1}\cdots {{\lambda }^{\prime }}_{k-1}^{a_{k-1}}={{\lambda }^{\prime }}_0^{{a^{\prime }}_0}{{\lambda }^{\prime }}_1^{{a^{\prime }}_1}\cdots {{\lambda }^{\prime }}_{k-1}^{{a^{\prime }}_{k-1}}\Rightarrow \forall p{a}_p={a^{\prime }}_p$ (3.26)

Dividing by the primed side:

${{\lambda }^{\prime }}_0^{a_0-{a^{\prime }}_0}{{\lambda }^{\prime }}_1^{a_1-{a^{\prime }}_1}\cdots {{\lambda }^{\prime }}_{k-1}^{a_{k-1}-{a^{\prime }}_{k-1}}=1$ (3.27)

Therefore, the condition is satisfied if and only if the only solution $a_p-{a^{\prime }}_p\in \mathbb{Z}$ to (3.27) is $a_p-{a^{\prime }}_p=0$ for all p. By rewriting the eigenvalues in the polar form ${{\lambda }^{\prime }}_p=r_p{\text{e}}^{i{\theta }_p}$ :

$(\begin{array}{l}{r}_0^{a_0-{a^{\prime }}_0}{r}_1^{a_1-{a^{\prime }}_1}\cdots r_{k-1}^{a_{k-1}-{a^{\prime }}_{k-1}}=1\\ {\text{e}}^{i{\theta }_0\left(a_0-{a^{\prime }}_0\right)}{\text{e}}^{i{\theta }_1\left(a_1-{a^{\prime }}_1\right)}\cdots {\text{e}}^{i{\theta }_{k-1}\left(a_{k-1}-{a^{\prime }}_{k-1}\right)}=1\end{array}$ (3.28)

and taking the logarithm of both sides of (3.28), we obtain a set of linear equations:

$(\begin{array}{l}\left(a_0-{a^{\prime }}_0\right)+\left(a_1-{a^{\prime }}_1\right){\log }_{r_0}{r}_1+\cdots +\left(a_{k-1}-{a^{\prime }}_{k-1}\right){\log }_{r_0}{r}_{k-1}=0\\ {\theta }_0\left(a_0-{a^{\prime }}_0\right)+{\theta }_1\left(a_1-{a^{\prime }}_1\right)+\cdots +{\theta }_{k-1}\left(a_{k-1}-{a^{\prime }}_{k-1}\right)=2\pi t,t\in \mathbb{Z}\end{array}$ (3.29)

There always exists a trivial solution $\forall p{a}_p-{a^{\prime }}_p=0$.

In case $k=2$, the equations reduce to:

$(\begin{array}{l}\left(a_0-{a^{\prime }}_0\right)+\left(a_1-{a^{\prime }}_1\right){\log }_{r_0}{r}_1=0\\ {\theta }_0\left(a_0-{a^{\prime }}_0\right)+{\theta }_1\left(a_1-{a^{\prime }}_1\right)=2\pi t\mathrm{,}t\in \mathbb{Z}\end{array}$ (3.30)

The upper equation has integer solutions if and only if ${\log }_{r_0}{r}_1\in \mathbb{Q}$. Substituting the result to the lower equation, we conclude that the above system of equations does not have an integer solution if:

${\log }_{r_0}{r}_1\notin \mathbb{Q}\vee \frac{\pi }{{\theta }_1-{\theta }_0{\log }_{r_0}{r}_1}\notin \mathbb{Q}$ (3.31)

If the transition matrix ${\tilde{T}}_{ab}$ can be brought into an upper triangular form, Lemmas 2.4 and 2.5 can be used to diagonalize it. This allows one to give an explicit expression for powers of ${\tilde{T}}_{ab}$, thus solving the recurrence. We now prove the main theorem of the article:

Theorem 3.4. Let $u_i^n\in ℂ\mathrm{,}i\ge \mathrm{0,1}\le n\le k$ be k sequences of numbers satisfying the recurrence relation (3.2) and let A and B denote the parameters of a linear transformation that brings the transition matrix to an upper triangular form, i.e. $k\times k$ matrix A and k dimensional vector B such that:

$B_p={\displaystyle \underset{s=0}{\overset{m}{\sum }}}\frac{1}{s!}{\displaystyle \underset{l_1,\cdots ,l_s=0}{\overset{k-1}{\sum }}}{\displaystyle \underset{\sigma }{\sum }}{}_p{}^{s}C{}_{l_{{\sigma }_1},\cdots ,l_{{\sigma }_s}}{\displaystyle \underset{t=1}{\overset{s}{\prod }}}{B}_{l_t}$ (3.32)

Furthermore, let the $k\times k$ array ${}_p{}^1{C}^{″}{}_x$ be upper triangular, where:

$\begin{array}{l}{}_p{}^1{C}^{″}{}_x={\displaystyle \underset{j,l_1=0}{\overset{k-1}{\sum }}}{A}_{pj}^{-1}{}_j{}^1{C}^{\prime }{}_l{A}_{lx}\\ {}_p{}^1{C}^{\prime }{}_{l_1}={\displaystyle \underset{s=1}{\overset{m}{\sum }}}{\displaystyle \underset{l_2,\cdots ,l_s=0}{\overset{k-1}{\sum }}}\frac{1}{\left(s-1\right)!}{\displaystyle \underset{\sigma }{\sum }}{}_p{}^{s}C{}_{l_{{\sigma }_1},\cdots ,l_{{\sigma }_s}}\left({\displaystyle \underset{t=2}{\overset{s}{\prod }}}{B}_{l_t}\right)\end{array}$ (3.33)

Let the eigenvalues ${}_p{}^1{C}^{″}{}_x$ (by eigenvalues of the array we understand the eigenvalues of the corresponding matrix created by considering only two lower indices) ${{\lambda }^{\prime }}_p$ be such that the equation:

$(\begin{array}{l}{x}_0+x_1{\log }_{\left|{\lambda }_0\right|}\left|{\lambda }_1\right|+\cdots +x_{k-1}{\log }_{\left|{\lambda }_0\right|}\left|{\lambda }_{k-1}\right|=0\\ Arg\left({\lambda }_0\right)x_0+Arg\left({\lambda }_1\right)x_1+\cdots +Arg\left({\lambda }_{k-1}\right)x_{k-1}=2\pi t\end{array}$ (3.34)

has no non-trivial solutions $x_0\mathrm{,}\cdots \mathrm{,}{x}_{k-1}\mathrm{,}t\in \mathbb{Z}$. Then the solution to the recurrence is given in terms of the initial condition ${\tilde{y}}^0$ by:

$u_i^k=B_{k-1}+{\displaystyle \underset{s=1}{\overset{k}{\sum }}}{\displaystyle \underset{l=1}{\overset{\infty }{\sum }}}{\displaystyle \underset{j=1}{\overset{l}{\sum }}}{A}_{k-1,s-1}{{\tilde{v}}^{\prime }}_j^s{{\tilde{T}}^{\prime }}_{jj}^i{{\tilde{P}}^{\prime }}_{jl}^{-1}{{\tilde{y}}^{\prime }}_l^0$ (3.35)

where:

$\begin{array}{l}{P^{\prime }}_{l_0,l_{p+1}}^{-1}=(\begin{array}{l}0\text{for}{l}_0>l_{p+1}\\ 1\text{for}{l}_0=l_{p+1}\\ {\displaystyle {\sum }_{p=0}^{l_{p+1}-l_0-1}{\displaystyle {\sum }_{l_0<l_1<\cdots <l_p<l_{p+1}}{\left(-1\right)}^{p+1}}}{\displaystyle {\prod }_{j=0}^p\frac{\tilde{v}{\text{'}}_{l_{j+1}}^{l_j}}{\tilde{v}{\text{'}}_{l_j}^{l_j}}}\text{for}{l}_0<l_{p+1}\end{array}\\ {v^{\prime }}_{l_{p+1}}^{l_0}=(\begin{array}{l}0\text{for}{l}_0>l_{p+1}\\ 1\text{for}{l}_0=l_{p+1}\\ {\displaystyle {\sum }_{p=0}^{l_{p+1}-l_0-1}{\displaystyle {\sum }_{l_0<l_1<\cdots <l_p<l_{p+1}}{\left(-1\right)}^{p+1}}}{\displaystyle {\prod }_{j=0}^p\frac{{{\tilde{T}}^{\prime }}_{l_j,l_{j+1}}}{{{\tilde{T}}^{\prime }}_{l_j,l_j}-{{\lambda }^{\prime }}_{l_{p+1}}}}\text{for}{l}_0<l_{p+1}\end{array}\end{array}$ (3.36)

Proof. Conditions (3.32) and (3.33) ensure that after the linear transformation ${z^{\prime }}^i=A^{-1}\left(z^i-B\right)$ the corresponding transition matrix ${T^{\prime }}_{ab}$ is upper triangular. Similarly, condition (3.34) is simply a sufficient condition for the diagonalization of ${{\tilde{T}}^{\prime }}_{ab}$ stated in Lemma 3.3.

Starting from Theorem 3.1, we have ${u^{\prime }}_i^k=e_k{T^{\prime }}^i{y}^0$. However, as mentioned before, due to redundancy of description in terms of $T_{ab}$, the corresponding ${\tilde{T}}_{ab}$ is used instead. Furthermore, since there is no redundancy in the first $k+1$ elements of $y^i$, the solution can be equivalently written as:

${u^{\prime }}_i^k=e_k{T^{\prime }}^i{y}^0=e_k{{\tilde{T}}^{\prime }}^i{\tilde{y}}^0={\displaystyle \underset{l=1}{\overset{\infty }{\sum }}}{\displaystyle \underset{j=1}{\overset{l}{\sum }}}{{\tilde{v}}^{\prime }}_j^k{{\tilde{T}}^{\prime }}_{jj}^i{{\tilde{P}}^{\prime }}_{jl}^{-1}{{\tilde{y}}^{\prime }}_l^0$ (3.37)

where the symbols with tilde are symbols related to ${\tilde{T}}_{ab}$ free from the redundancy mentioned before. The prime refers to the quantities related to the recurrence obtained after the linear transformation of variables. Transforming (3.37) back to the original sequences:

$u_i^k={\displaystyle \underset{s=1}{\overset{k}{\sum }}}{A}_{k-1,s-1}{u^{\prime }}_i^s+B_{k-1}=B_{k-1}+{\displaystyle \underset{s=1}{\overset{k}{\sum }}}{\displaystyle \underset{l=1}{\overset{\infty }{\sum }}}{\displaystyle \underset{j=1}{\overset{l}{\sum }}}{A}_{k-1,s-1}{{\tilde{v}}^{\prime }}_j^s{{\tilde{T}}^{\prime }}_{jj}^i{{\tilde{P}}^{\prime }}_{jl}^{-1}{{\tilde{y}}^{\prime }}_l^0$ (3.38)

Thus, we arrive at expression (3.35). The eigenvectors and the inverse of the infinite modal matrix $P^{-1}$ are given in Lemmas 2.4 and 2.5.

The solutions to the recurrence (3.2) have a general form:

$u_i={\displaystyle \underset{j_1,\cdots ,j_k=0}{\overset{\infty }{\sum }}}{f}_{j_1,\cdots ,j_k}\left(i\right){\displaystyle \underset{s=1}{\overset{k}{\prod }}}{\left(u_0^s\right)}^{j_s}$ (3.39)

which closely resembles Formula (2.37). A noticeable difference is that for the multi-variable recurrence, the solution takes the form of a multi-variable (instead of uni-variable) power series in the initial condition.

Example 3.5. As an example, we consider the simple recurrence defined by:

$(\begin{array}{l}{u}_i^1=8u_{i-1}^1+10u_{i-1}^2+{\left(u_{i-1}^1\right)}^2+3u_{i-1}^1{u}_{i-1}^2+{\left(u_{i-1}^2\right)}^2\\ u_i^2=-3u_{i-1}^1-3u_{i-1}^2+{\left(u_{i-1}^1\right)}^2-u_{i-1}^1{u}_{i-1}^2+{\left(u_{i-1}^2\right)}^2\end{array}$ (3.40)

The arrays ${}_p{}^{j}C{}_{l_1,\cdots ,l_j}$ can be obtained by comparing (3.40) with (3.2):

${}_p{}^{0}C=\left[\begin{array}{c}0\\ 0\end{array}\right],{}_p{}^{1}C{}_l=\left[\begin{array}{cc}8& 10\\ -3& -3\end{array}\right],{}_0{}^{2}C{}_{l_1,l_2}=\left[\begin{array}{cc}1& 3\\ 0& 1\end{array}\right],{}_1{}^{2}C{}_{l_1,l_2}=\left[\begin{array}{cc}1& -1\\ 0& 1\end{array}\right]$ (3.41)

In order to bring the transition matrix into the upper triangular form, the variables are transformed by a $2\times 2$ matrix A:

$A=\left[\begin{array}{cc}1& 2\\ -3& -5\end{array}\right]$ (3.42)

to obtain the new recurrence:

$(\begin{array}{l}{u^{\prime }}_i^1=2{u^{\prime }}_{i-1}^1+87{\left({u^{\prime }}_{i-1}^1\right)}^2+67{u^{\prime }}_{i-1}^1{u^{\prime }}_{i-1}^2+13{\left({u^{\prime }}_{i-1}^2\right)}^2\\ {u^{\prime }}_i^2=3{u^{\prime }}_{i-1}^2-212{\left({u^{\prime }}_{i-1}^1\right)}^2-164{u^{\prime }}_{i-1}^1{u^{\prime }}_{i-1}^2-32{\left({u^{\prime }}_{i-1}^2\right)}^2\end{array}$ (3.43)

with ${}_p{}^0{C}^{\prime }=0$ and:

${}_p{}^1{C}^{\prime }{}_l=\left[\begin{array}{cc}2& 0\\ 0& 3\end{array}\right],{}_0{}^2{C}^{\prime }{}_{l_1,l_2}=\left[\begin{array}{cc}87& 67\\ 0& 13\end{array}\right],{}_1{}^2{C}^{\prime }{}_{l_1,l_2}=\left[\begin{array}{cc}-212& -164\\ 0& -32\end{array}\right]$ (3.44)

The eigenvalues satisfy conditions (3.34), therefore ${{\tilde{T}}^{\prime }}_{ab}$ can be diagonalized and the powers of the transition matrix can be calculated. The infinite transition matrix is:

$T_{ab}=\left[\begin{array}{ccccccccccccc}1& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& \cdots \\ 0& 2& 0& 87& 67& 0& 13& 0& 0& 0& 0& 0& \cdots \\ 0& 0& 3& -212& -164& 0& -32& 0& 0& 0& 0& 0& \cdots \\ 0& 0& 0& 4& 0& 0& 0& 348& 268& 0& 52& 0& \cdots \\ 0& 0& 0& 0& 6& 0& 0& -424& 67& 0& 137& 0& \cdots \\ 0& 0& 0& 0& 0& 6& 0& -424& 67& 0& 137& 0& \cdots \\ 0& 0& 0& 0& 0& 0& 9& 0& -1272& 0& -987& 0& \cdots \\ 0& 0& 0& 0& 0& 0& 0& 8& 0& 0& 0& 0& \cdots \\ 0& 0& 0& 0& 0& 0& 0& 0& 12& 0& 0& 0& \cdots \\ 0& 0& 0& 0& 0& 0& 0& 0& 0& 12& 0& 0& \cdots \\ 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 18& 0& \cdots \\ 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 12& \cdots \\ ⋮& ⋮& ⋮& ⋮& ⋮& ⋮& ⋮& ⋮& ⋮& ⋮& ⋮& ⋮& \ddots \end{array}\right]$ (3.45)

The 5-th, 9-th and 11-th rows and columns are examples of the arising redundancy. Lemmas 2.4 and 2.5 can now be applied to the reduced transition matrix ${\tilde{T}}_{ab}$, which is obtained by eliminating all of the redundant rows and columns:

${\tilde{T}}_{ab}=\left[\begin{array}{ccccccc}1& 0& 0& 0& 0& 0& \cdots \\ 0& 2& 0& 87& 67& 13& \cdots \\ 0& 0& 3& -212& -164& -32& \cdots \\ 0& 0& 0& 4& 0& 0& \cdots \\ 0& 0& 0& 0& 6& 0& \cdots \\ 0& 0& 0& 0& 0& 9& \cdots \\ ⋮& ⋮& ⋮& ⋮& ⋮& ⋮& \ddots \end{array}\right]$ (3.46)

This can be further rewritten in the transformed form:

$\begin{array}{c}{\tilde{T}}_{ab}=\left[\begin{array}{ccccccc}1& 0& 0& 0& 0& 0& \cdots \\ 0& 1& 0& 87& 201& 39& \cdots \\ 0& 0& 1& -424& -656& -112& \cdots \\ 0& 0& 0& 2& 0& 0& \cdots \\ 0& 0& 0& 0& 12& 0& \cdots \\ 0& 0& 0& 0& 0& 21& \cdots \\ ⋮& ⋮& ⋮& ⋮& ⋮& ⋮& \ddots \end{array}\right]\left[\begin{array}{ccccccc}1& 0& 0& 0& 0& 0& \cdots \\ 0& 2& 0& 0& 0& 0& \cdots \\ 0& 0& 3& 0& 0& 0& \cdots \\ 0& 0& 0& 4& 0& 0& \cdots \\ 0& 0& 0& 0& 6& 0& \cdots \\ 0& 0& 0& 0& 0& 9& \cdots \\ ⋮& ⋮& ⋮& ⋮& ⋮& ⋮& \ddots \end{array}\right]\\ \times \left[\begin{array}{ccccccc}1& 0& 0& 0& 0& 0& \cdots \\ 0& 1& 0& -\frac{87}{2}& -\frac{67}{4}& -\frac{13}{7}& \cdots \\ 0& 0& 1& 212& \frac{164}{3}& \frac{16}{3}& \cdots \\ 0& 0& 0& \frac{1}{2}& 0& 0& \cdots \\ 0& 0& 0& 0& \frac{1}{12}& 0& \cdots \\ 0& 0& 0& 0& 0& \frac{1}{21}& \cdots \\ ⋮& ⋮& ⋮& ⋮& ⋮& ⋮& \ddots \end{array}\right]\end{array}$ (3.47)

The matrices to the left and right of the diagonal matrix are $P$ and $P^{-1}$ respectively. The powers of ${\tilde{T}}_{ab}$ can now be easily calculated. Using (3.5), we finally arrive at:

$\begin{array}{l}{u^{\prime }}_i^1=2^i{u^{\prime }}_0^1+\frac{87}{2}\left(4^i-2^i\right){\left({u^{\prime }}_0^1\right)}^2+\frac{67}{4}\left(6^i-2^i\right){u^{\prime }}_0^1{u^{\prime }}_0^2+\frac{13}{7}\left(9^i-2^i\right){\left({u^{\prime }}_0^2\right)}^2+\cdots \\ {u^{\prime }}_i^2=3^i{u^{\prime }}_0^2-212\left(4^i-3^i\right){\left({u^{\prime }}_0^1\right)}^2-\frac{164}{3}\left(6^i-3^i\right){u^{\prime }}_0^1{u^{\prime }}_0^2-\frac{16}{3}\left(9^i-3^i\right){\left({u^{\prime }}_0^2\right)}^2+\cdots \end{array}$ (3.48)

The above solution to the primed recurrence (3.43) can now be transformed back to the original variables, therefore giving a rather complicated expression for the solution to recurrence (3.40):

$\begin{array}{c}{u}_i^1=\left(6\cdot 3^i-5\cdot 2^i\right)u_0^1+10\left(3^i-2^i\right)u_0^2\\ +\left(\frac{87}{7}{9}^i-\frac{307}{4}{6}^i+\frac{413}{2}{4}^i-192\cdot 3^i+\frac{1395}{28}{2}^i\right){\left(u_0^1\right)}^2\\ +\left(\frac{290}{7}{9}^i-\frac{3377}{12}{6}^i+826\cdot 4^i-\frac{2440}{3}{3}^i+\frac{6365}{28}{2}^i\right)u_0^1{u}_0^2\\ +\left(\frac{725}{21}{9}^i-\frac{1535}{6}{6}^i+826\cdot 4^i-\frac{2608}{3}{3}^i+\frac{3705}{14}{2}^i\right){\left(u_0^2\right)}^2+\cdots \end{array}$

$\begin{array}{c}{u}_i^2=-3\left(3^i-2^i\right)u_0^1-\left(5\cdot 3^i-6\cdot 2^i\right)u_0^2\\ +\left(\frac{15}{7}{9}^i+\frac{53}{4}{6}^i-\frac{163}{2}{4}^i+96\cdot 3^i-\frac{837}{28}{2}^i\right){\left(u_0^1\right)}^2\\ +\left(\frac{50}{7}{9}^i+\frac{583}{12}{6}^i-326\cdot 4^i+\frac{1220}{3}{3}^i-\frac{3819}{28}{2}^i\right)u_0^1{u}_0^2\\ +\left(\frac{125}{21}{9}^i+\frac{265}{6}{6}^i-326\cdot 4^i+\frac{1304}{3}{3}^i-\frac{2223}{14}{2}^i\right){\left(u_0^2\right)}^2+\cdots \end{array}$ (3.49)

Comparing (3.49) with the general form of solution (3.39), the coefficient functions can be obtained.

4. Notes on Multi-Variable Recurrences of Arbitrary Depth

Theorem 3.4 is also applicable to arbitrary-depth recursions. For any finite-depth recursion, auxiliary variables can be introduced and the original recurrence can be reduced to multi-variable recurrence of depth one. The method follows closely the one used for the finite-depth linear recurrences [4].

Let the depth-n recurrence be in the form (1.1). Set of $k\times n$ auxiliary variables is defined by:

$u_{i-j}^l=u_i^{l+jk}$ (4.1)

for $0\le j<n$. Equation (1.1) now takes the form:

$\{\begin{array}{l}{u}_i^1=F_1\left(u_{i-1}^1,\cdots ,u_{i-1}^{kn}\right)\\ ⋮\\ u_i^k=F_k\left(u_{i-1}^1,\cdots ,u_{i-1}^{kn}\right)\\ u_i^{k+1}=u_i^1\\ ⋮\\ u_i^{kn}=u_i^{k\left(n-1\right)}\end{array}$ (4.2)

Which is a multi-variable depth one recurrence, already discussed in Section 3.

5. Further Remarks

In Section 2, an explicit formula for the solution of wide range of uni-variable depth-one recurrences is given. Although the final expression is rather complicated, it can be used to derive the coefficient functions $f_j\left(i\right)$ in (2.37). Even if the recurrence does not fulfill the diagonalizability condition, there is a chance that it can be reformulated in terms of new variables. The transformation should be chosen such that the corresponding transition matrix has different eigenvalues and Lemmas 2.4 and 2.5 can be applied. The shift transformation is of particular importance, nevertheless, any transformation that leaves the recurrence polynomial could be used. Transformations by a higher degree invertible polynomial, such us ${u^{\prime }}_i=u_i^3+1$ are just one of many examples.

Although the method allows tackling many different polynomial recurrence relations, there are some notable shortcomings. The diagonalizability conditions in Sections 2 and 3 are only sufficient and not necessary for the transition matrix to be diagonalizable. Therefore, the method may not be applicable for all systems that are solvable by diagonalization approach. In addition, transition matrices which are known to be non-diagonalizabe are yet another group of cases where Theorems 2.6 and 3.4 cannot be used. Possibly, an analogous method could be formulated, where Jordanization is used to calculate the powers of the transition matrix instead of diagonalization.

Theorem 3.4 utilizes the non-redundant system (with $\tilde{T}$ ) instead of the one defined in Theorem 3.1. It is dictated by the fact that the transition matrix defined for the auxiliary vectors 3.4 necessary contains redundant copies of its eigenvalues. This may lead to difficulties with evaluating the expansion (3.39) for large $j_1\mathrm{,}\cdots \mathrm{,}{j}_k$, as only the general rule for constructing $\tilde{T}$ is given as opposed to explicit expression as in the case of T, i.e. (3.9).

It can be noted that examples 2.7 and 3.5 belong to a family of recurrences for which no shift is needed. The logistic map automatically fulfills the assumptions for $r\ne 1$, while the recurrence relation (3.40) needs to be transformed only by matrix A, with $B=0$. Those examples were chosen on purpose, as the solution of the transformed series can be truncated and the calculations simplify significantly. In general, however, if one wants to obtain the exact solution to the original recurrence with non-zero shift parameter, all of the terms of the solution corresponding to the transformed recurrence have to be kept, which often makes the calculations cumbersome.

Even though the Carleman linearization method does not provide a way to systematically approach the topic of polynomial recurrences, a large group of relations (1.1) fulfill the conditions and the method can be applied. Furthermore, it is straightforward to generalize the method to systems of recurrences (Theorem 3.4), as well as systems of arbitrary-depth recurrences, as discussed in Section 4.

Improvements to the Carleman linearization of polynomial recurrences may include research on the diagonalizability of the transition matrices, possibly resulting in less-restrictive diagonalizability conditions. In addition, polynomial transformations may play an important role in further extending the group of recurrences the method is applicable to.

Acknowledgements

The author is grateful to Daniel Gagliardi for valuable remarks, and to the anonymous referees for devoting their time and effort to the manuscript.

Conflicts of Interest

The author declares no conflicts of interest regarding the publication of this paper.

References

[1]	Borovkov, A.A. (2013) Probability Theory. Springer, London. https://doi.org/10.1007/978-1-4471-5201-9
[2]	Mladenović, P. (2019) Combinatorics: A Problem-Based Approach. Springer, Cham. https://doi.org/10.1007/978-3-030-00831-4
[3]	Andrica, D. and Bagdasar, O. (2020) Recurrent Sequences: Key Results, Applications, and Problems. Springer, Cham. https://doi.org/10.1007/978-3-030-51502-7
[4]	Everest, G., Poorten, A., Shparlinski, I. and Ward, T. (2003) Recurrence Sequences. American Mathematical Society, Providence. https://doi.org/10.1090/surv/104
[5]	Zhang, X., Shi, Y. and Chen, G. (2009) Constructing Chaotic Polynomial Maps. International Journal of Bifurcation and Chaos, 19, 531-543. https://doi.org/10.1142/S0218127409023172
[6]	Grosjean, N. and Huillet, T. (2016) Some Combinatorial Aspects of Discrete Non-Linear Population Dynamics. Chaos, Solitons & Fractals, 93, 71-79. https://doi.org/10.1016/j.chaos.2016.10.004
[7]	Han, D.D., Min, L.Q., Zang, H.Y. and Yang, X.P. (2019) Robust Chaos of Cubic Polynomial Discrete Maps with Application to Pseudorandom Number Generators. Mathematical Problems in Engineering, 2019, Article ID: 8250903. https://doi.org/10.1155/2019/8250903
[8]	Wang, C.F. and Ding, Q. (2019) A Class of Quadratic Polynomial Chaotic Maps and Their Fixed Points Analysis. Entropy, 21, 658-671. https://doi.org/10.3390/e21070658
[9]	Rabinovich, S., Berkolaiko, G. and Havlin, S. (1996) Solving Nonlinear Recursions. Journal of Mathematical Physics, 37, Article No. 5828. https://doi.org/10.1063/1.531702
[10]	Shang, Y. (2012) A Brief Note on an Exponential Recursive Sequence. International Journal of Open Problems in Computer Science and Mathematics, 5, 5 p. https://doi.org/10.12816/0006093
[11]	Cadilhac, M., Mazowiecki, F., Paperman, C., Pilipczuk, M. and Sénizergues, G. (2021) On Polynomial Recursive Sequences. Theory of Computing Systems. https://doi.org/10.1007/s00224-021-10046-9
[12]	Cull, P., Flahive, M. and Robson, R. (2005) Difference Equations: From Rabbits to Chaos. Springer, New York.
[13]	Hogben, L. (2014) Handbook of Linear Algebra. 2nd Edition, Chapman and Hall/CRC, New York. https://doi.org/10.1201/b16113
[14]	Berkolaiko, G., Rabinovich, S. and Havlin, S. (1998) Analysis of Carleman Representation of Analytical Recursions. Journal of Mathematical Analysis and Applications, 224, 81-90. https://doi.org/10.1006/jmaa.1998.5986
[15]	Kowalski, K. and Steeb, W.H. (1991) Nonlinear Dynamical Systems and Carleman Linearization. World Scientific, Singapore. https://doi.org/10.1142/1347
[16]	Gralewicz, P. and Kowalski, K. (2002) Continuous Time Evolution from Iterated Maps and Carleman Linearization. Chaos, Solitons & Fractals, 14, 563-572. https://doi.org/10.1016/S0960-0779(01)00222-3
[17]	Rabinovich, S., Berkolaiko, G., Buldyrev, S., Shehter, A. and Havlin, S. (1995) ‘Logistic Map’: An Analytical Solution. Physica A: Statistical Mechanics and its Applications, 218, 457-460. https://doi.org/10.1016/0378-4371(95)00163-2
[18]	Belozyorov, V.Y. and Volkova, S.A. (2016) Role of Logistic and Ricker’s Maps in Appearance of Chaos in Autonomous Quadratic Dynamical Systems. Nonlinear Dynamics, 83, 719-729. https://doi.org/10.1007/s11071-015-2360-2
[19]	Tauber, S. (1963) On Multinomial Coefficients. The American Mathematical Monthly, 70, 1058-1063. https://doi.org/10.1080/00029890.1963.11992172
[20]	Shang, Y. (2011) A Remark on the Chromatic Polynomials of Incomparability Graphs of Posets. International Journal of Pure and Applied Mathematics, 67, 159-164.
[21]	Basu, S. and Velleman, D.J. (2017) On Gauss’s First Proof of the Fundamental Theorem of Algebra. The American Mathematical Monthly, 124, 688-694. https://doi.org/10.4169/amer.math.monthly.124.8.688
[22]	Herrero, D.A. (1991) Triangular Operators. Bulletin of the London Mathematical Society, 23, 513-554. https://doi.org/10.1112/blms/23.6.513
[23]	Bronson, R. and Costa, G.B. (2020) Matrix Methods: Applied Linear Algebra and Sabermetrics. Academic Press, London.
[24]	Ford, W. (2015) Numerical Linear Algebra with Applications: Using MATLAB. Academic Press, London.
[25]	Baliarsingh, P. and Dutta, S. (2015) On an Explicit Formula for Inverse of Triangular Matrices. Journal of the Egyptian Mathematical Society, 23, 297-302. https://doi.org/10.1016/j.joems.2014.06.001
[26]	Strogatz, S.H. (2015) Nonlinear Dynamics and Chaos: With Applications to Physics, Biology, Chemistry, and Engineering. CRC Press, Boca Raton.
[27]	Comon, P., Golub, G., Lim, L. and Mourrain, B. (2008) Symmetric Tensors and Symmetric Tensor Rank. SIAM Journal on Matrix Analysis and Applications, 30, 1254-1279. https://doi.org/10.1137/060661569