Denominator Bounds and Polynomial Solutions for Systems of q-Recurrences over K(t) for Constant K

TL;DR

This paper extends methods for bounding the denominators and degrees of polynomial solutions to systems of higher order q-recurrence equations with rational coefficients, without relying on uncoupling or reduction techniques.

Contribution

It generalizes scalar case techniques to systems of q-recurrences, providing direct methods for bounds without uncoupling or reduction.

Findings

Provides bounds on the maximal power of t in denominators of solutions

Establishes degree bounds for polynomial solutions

Offers a direct approach that avoids uncoupling or reduction

Abstract

We consider systems A_\ell(t) y(q^\ell t) + ... + A_0(t) y(t) = b(t) of higher order q-recurrence equations with rational coefficients. We extend a method for finding a bound on the maximal power of t in the denominator of arbitrary rational solutions y(t) as well as a method for bounding the degree of polynomial solutions from the scalar case to the systems case. The approach is direct and does not rely on uncoupling or reduction to a first order system. Unlike in the scalar case this usually requires an initial transformation of the system.

Equations

$$A_s{\sigma }^s(y)+\dots +A_1\sigma (y)+A_{0}y=b$$

$$A_s{\sigma }^s(y)+\dots +A_1\sigma (y)+A_{0}y=t^{-\nu }b$$

$$\sum _{j=0}^s\frac{A_j{\sigma }^j(u)}{{\sigma }^j(g)}+\sum _{j=0}^s\frac{A_j{\sigma }^j(v)}{q^{nj}{t}^n}=t^{-\nu }b$$

$$t^{\nu-n}\Bigl{(}\prod_{j=0}^{s}\sigma^{j}(g)\Bigr{)}\sum_{j=0}^{s}q^{-nj}A_{j}\sigma^{j}(v)\\ =\Bigl{(}\prod_{j=0}^{s}\sigma^{j}(g)\Bigr{)}b-\sum_{j=0}^{s}\Bigl{(}\prod_{k\neq j}\sigma^{k}(g)\Bigr{)}A_{j}\sigma^{j}(u).$$

$$t\Bigl{.}\;\Bigm{|}\;\Bigr{.}\sum_{j=0}^{s}q^{-nj}A_{j}\sigma^{j}(v).$$

$$0=\Bigl{(}\sum_{j=0}^{s}q^{-nj}A_{j0}\Bigr{)}v_{0}.$$

$$\begin{gathered} \left(\begin{array}{cc}8& 0\\ 8& 0\end{array}\right)\left(\begin{array}{c}{y}_1(4t)\\ y_2(4t)\end{array}\right)+\left(\begin{array}{cc}-16t+4& 8\\ -16t^2+16t-12& 8\end{array}\right)\left(\begin{array}{c}{y}_1(2t)\\ y_2(2t)\end{array}\right) \\ +\left(\begin{array}{cc}16t-4& -8t^3-1\\ 16t^2-8t+4& -8t^4-1\end{array}\right)\left(\begin{array}{c}{y}_1(t)\\ y_2(t)\end{array}\right)=0. \end{gathered}$$

$$\lambda =\det \left(\begin{array}{cc}8x^2+4x-4& 8x-1\\ 8x^2-12x+4& 8x-1\end{array}\right)=128x^2-80x+8.$$

$$\left(\begin{array}{c}1\\ t^{-3}\end{array}\right)\text{and}\left(\begin{array}{c}{t}^{-1}\\ t^{-3}\end{array}\right).$$

$$y=y_n{t}^n+\dots +y_{1}t+y_0$$

$$\sum_{j=0}^{s}\Bigl{(}\sum_{i=0}^{\ell}A_{ji}t^{i}\Bigr{)}\sigma^{j}\Bigl{(}\sum_{k=0}^{n}y_{k}t^{k}\Bigr{)}\\ =\sum_{j=0}^{s}\sum_{i=0}^{\ell}\sum_{k=0}^{n}q^{jk}t^{i+k}A_{ji}y_{k}=t^{-\nu}\sum_{j=0}^{\varkappa}b_{j}t^{j}.$$

$$\Bigl{(}\sum_{j=0}^{s}A_{j\ell}q^{nj}\Bigr{)}y_{n}=b_{\ell+n+\nu}.$$

$$\begin{gathered} \left(\begin{array}{cc}1& 0\\ 1& 0\end{array}\right)\left(\begin{array}{c}{z}_1(4t)\\ z_2(4t)\end{array}\right)+\left(\begin{array}{cc}4-16t& 8\\ 12+16t-16t^2& 8\end{array}\right)\left(\begin{array}{c}{z}_1(2t)\\ z_2(2t)\end{array}\right) \\ +\left(\begin{array}{cc}128t-32& -64t^3-8\\ 128t^2-64t+32& -64t^4-8\end{array}\right)\left(\begin{array}{c}{z}_1(t)\\ z_2(t)\end{array}\right)=0. \end{gathered}$$

$$\varrho =\det \left(\begin{array}{cc}0& 0\\ 0& -64\end{array}\right)=0.$$

$$\text{K}[t,\sigma ]=\{{\alpha }_d(t){\sigma }^d+\dots {\alpha }_1(t)\sigma +{\alpha }_0\mid {\alpha }_0,\dots ,{\alpha }_d\in \text{K}[t]\}$$

$$\sigma t=qt\sigma$$

$$A\bullet y=t^{-\nu }b.$$

$$(QA)\bullet y=Q\bullet t^{-\nu }b$$

$$QM=\left(\begin{array}{c}P\\ 0\end{array}\right)$$

$$AU=\left(\begin{array}{cc}\tilde{A}& 0\end{array}\right)$$

$$QAU=\left(\begin{array}{cc}\text{accentset}\approx A& 0\\ 0& 0\end{array}\right)$$

$$x=Q\left(\begin{array}{c}y\\ z\end{array}\right)\text{and}Q\bullet b=\left(\begin{array}{c}c\\ w\end{array}\right)$$

$$A=t^{\ell }{\tilde{A}}_{\ell }+\dots +t{\tilde{A}}_1+{\tilde{A}}_0$$

$${\tilde{A}}_0=\sum _{j=0}^s{A}_{j0}{\sigma }^j\text{and}{\tilde{A}}_{\ell }=\sum _{j=0}^s{A}_{j\ell }{\sigma }^j.$$

$${\tilde{A}}_0=\left(\begin{array}{c}{\tilde{B}}_0\\ 0\end{array}\right)$$

$${\tilde{A}}_0=\left(\begin{array}{c}{\tilde{B}}_0\\ {\tilde{C}}_0\end{array}\right)$$

$$Q=\left(\begin{array}{cc}{1}_r& 0\\ -X& 1_{m-r}\end{array}\right)$$

$$A_0=\left(\begin{array}{c}{B}_0\\ C_0\end{array}\right)$$

$$\begin{gathered} A=\left(\begin{array}{c}B\\ C\end{array}\right)\stackrel{\text{reduction}}{\to }\left(\begin{array}{cc}1& 0\\ -X& 1\end{array}\right)\left(\begin{array}{c}B\\ C\end{array}\right) \\ \stackrel{\text{shift}}{\to }\left(\begin{array}{cc}1& 0\\ 0& t^{-1}1\end{array}\right)\left(\begin{array}{cc}1& 0\\ -X& 1\end{array}\right)\left(\begin{array}{c}B\\ C\end{array}\right)=\left(\begin{array}{cc}1& 0\\ -t^{-1}X& t^{-1}1\end{array}\right)\left(\begin{array}{c}B\\ C\end{array}\right). \end{gathered}$$

$$\left(\begin{array}{cc}1& 0\\ -t^{-\nu }X& t^{-\nu }\end{array}\right)$$

Full text

recurrence systems, $q$-recurrences, direct solving, denominator bound, polynomial solutions, Popov normal form, Ore operators

This work was supported by the Austrian Science Fund (FWF) grant SFB50 (F5009-N15).

††conference: International Symposium on Symbolic and Algebraic Computation; July 25-28, 2017; Kaiserslautern

1. Introduction

Let $\K$ be a field of characteristic [math], and let $q\in\K\setminus\{0\}$ be not a root of unity. Moreover, let $t$ be an indeterminate over $\K$. We define an automorphism $\sigma$ on the rational function field $\K(t)$ by setting $\sigma(t)=qt$. We will refer to $\sigma$ as the $q$-shift. Note that with this definition we have $\sigma(\alpha)=\alpha$ for all $\alpha\in\K$. This is an easy case of a homogeneous or $\Pi$-extension in the sense of (Karr81, ). In the following we will often let $\sigma$ act on (column) vectors in $\K(t)^{m}$. This is to be understood as component-wise application of $\sigma$.

We consider the recurrence equation

[TABLE]

where $A_{0},\ldots,A_{s}\in\K[t]^{m\times m}$ are square polynomial matrices and $b\in\K[t]^{m}$ is a polynomial vector. We will assume that the matrix $A=A_{s}\sigma^{s}+\ldots+A_{0}$ has full (left row) rank over the operator ring $\K(t)[\sigma]$ (see Section 3 for details). Note that any $q$-recurrence system can be brought into this form using unimodular transformations (see, for example, (myPhDthesis, , Lem. 7.5); it is easy to verify that the method discussed there works for general Ore operators – see Section 3 for the definition – and not just differential operators). We will briefly explain the idea at the end of Section 3.

Generally, one major area where recurrence equations arise is symbolic summation. See, for example, (Karr81, ), (SchneiderThesis, ), (AeqlB, ), or (ChyzakSalvy98, ). In particular, note that symbolic simplification of products in (Karr81, ) and (SchneiderThesis, ) is done within so-called homogeneous or $\Pi$-extension of which our field $\K(t)$ is a special case. Examples which deal explicitly with $q$-hypergeometric sums are (PauleRiese:97, ) or (Riese:03, ).

Our motivating goal is to find new algorithms to compute all rational solutions of (rec); that is, we want to know all $y\in\K(t)^{m}$ such that the equation holds for that $y$.

A first step will be to determine a so-called denominator bound: we are looking for a polynomial $d\in\K[t]\setminus\{0\}$ such that $dy\in\K[t]^{m}$ for every possible rational solution $y\in\K(t)^{m}$. There are comparatively easy ways to determine all the factors in $d$ except for powers of $t$; cf., for example, (AbramovBarkatou1998, ), (AbramovKhmelnov2012, ), (CPS:08, ), or (SchneiderMiddeke2017, ). See (Karr81, ), (Bron:00, ) or (SchneiderThesis, ) for a discussion on why this distinction is necessary. Thus, in this work we will only concentrate on finding a bound on the powers of $t$ in $d$ in a computationally cheap way. The next step will then be to compute polynomial solutions of the system. For this, we are computing a bound on the degree of polynomial solutions $\tilde{y}\in\K[t]^{m}$. Once such a bound is determined, a solution to (rec) can be found by making an ansatz with unknown coefficients and solving the resulting linear system over $\K$.

Methods for determining such a denominator bound (or universal denominator as some authors prefer to call it), have already been discussed in the literature. See, for example, (AbramovBarkatou1998, ), (Barkatou1999, ), (CPS:08, ), or (AbramovKhmelnov2012, ). Most of these approaches rely on the transformation of the system into first order; the method in (AbramovKhmelnov2012, ) seems to be the only one which deals with systems directly (albeit for shift recurrences and not the $q$-case; although it might be possible to adapt the method using ideas from (Abramov2002, )). Another way of determining a denominator bound would be to decouple the system using cyclic vector methods, the Danilevski-Barkatou-Zücher algorithm (see (Danilevski37, ), (Barkatou1993, ), (Zurcher94, ), or (BCP13, )), or Smith–Jacobson form computations.

In this paper, however, we want to avoid uncoupling and conversion to first order systems as these often introduce an additional cost and also tend to blur the structure of the input system. Instead we aim at computing the denominator bound and degree bound directly from the matrices of the higher order system (rec). Our hope here is that this will lead to more efficient algorithms.

This paper is structured as follows: In Section 2, we discuss how the powers of $t$ in the denominator can be bounded and how we may determine a bound for the degree of polynomial solutions. The results will require certain determinants not to be identically zero. Since this of course cannot be guaranteed for an arbitrary system, we need to develop an algorithm which transforms a given system into an equivalent one where the method of Section 2 is applicable. This is done in Section 4 which constitutes the major result of this paper. The section also includes a comparison with competing methods. In between, Section 3 recalls some basic facts about operators and matrix normal forms which are needed for the transformation.

2. Bounding the Power of the Indeterminate and the Degree of

Polynomial Solutions

Consider (rec) where $A_{0},\ldots,A_{s}$ are matrices in $\K[t]^{m\times m}$ and where $b\in\K[t]^{m}$ is a column vector. For technical reasons the usefulness of which will become apparent in Section 4, we consider the slightly modified system

[TABLE]

where $\nu\geqslant 0$, that is, where we allow the right hand side to be a rational vector but only with powers of $t$ in the denominator.

We start by considering possible denominators. The following method is folklore. Similar ideas can be found, for example, in (Karr81, ). Let $y\in\K(t)^{m}$ be an arbitrary solution of (rec’) and let $gt^{n}\in\K[t]$ be the denominator of $y$ where $t\nmid g$. Using an (entry-wise) partial fraction decomposition, we may write $y=u/g+v/t^{n}$ where $u,v\in\K[t]^{m}$ and where $v=0$ or $t\nmid v$ (that is, $t$ does not divide every entry of $v$). In the first case, $n=0$ is a candidate for the power of $t$ in the degree bound. In the other case, substituting this into (rec’) yields

[TABLE]

since $\sigma^{j}(t^{n})=q^{nj}t^{n}$ for all $j$. We may transform the equation into

[TABLE]

Note that the right hand side is in $\K[x]^{m}$, that is, a polynomial vector. That means that also the left hand side must be polynomial, and this in turn implies that either $\nu\geqslant n$ or else that $t$ divides \bigl{(}\prod_{j=0}^{s}\sigma^{j}(g)\bigr{)}\sum_{j=0}^{s}q^{-nj}A_{j}\sigma^{j}(v). Assume that the latter case holds. Then, since $t\nmid g$, we also have $t\nmid\sigma^{j}(g)$ for all $j$ (or else $t\sim q^{-j}t=\sigma^{-j}(t)\mid g$). Hence, $t$ cannot divide $\prod_{j=0}^{s}\sigma^{j}(g)$. We obtain

[TABLE]

In particular, the coefficient vector of $t^{0}$ in the sum must vanish.

Write now $A_{j}=A_{j\ell}t^{\ell}+\ldots+A_{j1}t+A_{j0}$ with $A_{j0},\ldots,A_{j\ell}\in\K^{m\times m}$ being constant matrices for all $j$; and write $v=v_{d}t^{d}+\ldots+v_{1}t+v_{0}$ with $v_{0},\ldots,v_{d}\in\K^{m}$. From the assumption that $t\nmid v$ we know that in particular $v_{0}\neq 0$. Using these representations in the sum above, we find that the constant term is

[TABLE]

Now, $\Lambda=\sum_{j=0}^{s}q^{-nj}A_{j0}\in\K^{m\times m}$ is just a matrix and since $v_{0}\neq 0$ we see that $\Lambda$ must have a non-trivial kernel. This implies that the determinant of $\Lambda$ vanishes. We may regard $\det\Lambda$ as a polynomial expression in $q^{-n}$. Therefore, assuming that $\det\Lambda$ does not identically vanish, we can obtain candidates for $n$ by checking if it has roots of the form $q^{-n}$.

This observation motivates the following definition.

Definition 2.1 ($t$-tail regular).

We call the system (rec’) $t$-tail regular if \lambda=\det\bigl{(}\sum_{j=0}^{s}A_{j0}x^{j}\bigr{)}\in\K[x] is not the zero polynomial.

Summarising the argumentation so far, we obtain Theorem 2.2 below. It remains to deal with the case of systems which are not $t$-tail regular. We will do so in Section 4.

Theorem 2.2.

Assume that the system (rec’) is $t$-tail regular. If $y\in\K(t)^{m}$ is a rational solution of (rec’) and if $t^{n}$ divides the denominator of $y$, then

1. (1)

$n=0$, 2. (2)

or $n\leqslant\nu$, 3. (3)

or $q^{-n}$ is a root of \lambda=\det\bigl{(}\sum_{j=0}^{s}A_{j0}x^{j}\bigr{)}\in\K[x].

Consequently, we obtain a bound on $n$ by taking the maximum of those values. ∎

If $\K$ happens to be of the form $\mathbb{F}(q)$ where $\mathbb{F}$ is a computable field and $q$ is transcendental over $\mathbb{F}$, then the method of (AbramovPaulePetkovsec, , Ex. 1) can be applied to find all roots of $\lambda$ of the required form. In the case that $\K=\mathbb{Q}$, we can use the method of (BauerPetkovsek1999, , Sect. 3.3) for that.

Example 2.3.

Let $\K=\mathbb{Q}$ and $q=2$. We consider the system

[TABLE]

Here, $\nu=0$. Computing the polynomial $\lambda$ yields

[TABLE]

The roots of $\lambda$ are $1/2=q^{-1}$ and $1/8=q^{-3}$. That yields the upper bound $3$ for $n$. Note that this fits the actual solutions well which are generated by

[TABLE]

(We can easily check that the two vectors are indeed solutions; and – after row reduction – (AbramovBarkatou2014, , Thm. 6) gives the dimension of the solution space as $2$.)

We turn our attention to polynomial solutions now. The approach is rather similar to the degree bounding and has been discussed in the literature before, for example, in (AbramovPaulePetkovsec, , Sect. 2). Assume that $y\in\K[t]^{m}$ is a solution of (rec’). We make an ansatz

[TABLE]

where $y_{0},\ldots,y_{n}\in\K^{m}$. We assume that $y_{n}\neq 0$, that is, that $n$ is the degree of $y$. In addition, let $b=b_{\varkappa}t^{\varkappa}+\ldots+b_{1}t+b_{0}$ with $b_{0},\ldots,b_{\varkappa}\in\K^{m}$; and write once more $A_{j}=A_{j\ell}t^{\ell}+\ldots+A_{j1}t+A_{j0}$ with $A_{j0},\ldots,A_{j\ell}\in\K^{m\times m}$ for all $j$. We assume that $b_{\varkappa}\neq 0$ and that $A_{j\ell}\neq\mathbf{0}$ for at least one $j$. Substituting all of this into (rec’) gives

[TABLE]

If the right hand side contains negative powers of $t$, then there cannot be a polynomial solution. Else, the degree (in $t$) of the left hand side is at most $\ell+n$. We compute the coefficient of $t^{\ell+n}$ on both sides. This yields (since $i=\ell$ and $k=n$)

[TABLE]

There are two cases:

1. (1)

Either $\varkappa\geqslant\ell+n+\nu$ (that is, $n\leqslant\varkappa-\nu-\ell$), 2. (2)

or $\varkappa<\ell+n+\nu$.

In the first case, we do already have a bound on $n$ (namely $\varkappa-\nu-\ell$). In the second case, the right hand side of the equation for the coefficient of $t^{\ell+n}$ is zero. Just as for the denominator bound, we see that the left hand side contains a matrix with a non-trivial kernel. Again, its determinant is a polynomial expression—this time in $q^{n}$. We make an analogous definition.

Definition 2.4 ($t$-head regular).

We call the system (rec’) $t$-head regular if the polynomial \varrho=\det\bigl{(}\sum_{j=0}^{s}A_{j\ell}x^{j}\bigr{)}\in\K[x] is not identically zero.

The computations above have show that for a $t$-head regular system $q^{n}$ is a root of $\varrho$. Thus, we can again obtain all the candidates for the degree $n$ of $y$ by computing the roots of $\varrho$. The following theorem summarises these findings.

Theorem 2.5.

Assume that the system (rec’) is $t$-head regular. Let $\ell$ be the maximum of $\deg_{t}A_{0},\ldots,\deg_{t}A_{s}$ and let $\varkappa=\deg_{t}b$. If $y\in\K[t]^{m}$ is a polynomial solution of degree $n$, then

1. (1)

either $n\leqslant\varkappa-\nu-\ell$, 2. (2)

or $q^{n}$ is a root of \varrho=\det\bigl{(}\sum_{j=0}^{s}A_{j\ell}x^{j}\bigr{)}\in\K[x].

As in Theorem 2.2, we obtain a finite number of candidates for $n$.

Example 2.6.

We continue Example 2.3. We have already seen that rational solutions have the shape $y=t^{-3}z$ for some $z\in\K[t]^{m}$.111We did not explicitly show that the denominator of $y$ does not have any other factors apart from $t$ – however, this is the case here since we do know the complete solution space. Substituting this into the equation and clearing denominators yields

[TABLE]

The maximal degree in $t$ is $\ell=4$ and thus, we have

[TABLE]

Thus, Theorem 2.5 is not applicable here. We will revisit the example in Section 4.

3. Ore Operators and Matrix Normal Forms

This section is devoted to introducing some concepts which are necessary in order to deal with the cases $\lambda=0$ or $\varrho=0$. This section introduces the algebraic model which we are going to employ as well as some necessary background information and other results.

In Section 4, we will try to transform the system (rec’) into an equivalent form such that the new $\lambda$ or $\varrho$ becomes non-zero. During this, we are allowed to use the following transformations:

1. (1)

multiply a row of the system by a rational function in $\K(t)$; 2. (2)

add a $\K(t)$-linear combination of shifted versions of one row to another row.

We will apply the transformations in such a way that the coefficients on the left hand side of (rec’) will always be polynomial matrices. However, the right hand side might contain rational functions in $t$.

It will be convenient to express the transformations using the language of Ore operators. Roughly speaking, these are a class of non-commutative polynomials. They are named after their first introduction by Ø. Ore in (Ore33, ); a first interpretation of them as operators seems to be (JacPLT, ). Ore operators have been used to model differential or difference equations (see (BP96, ) for an introduction) and have been applied to problems in symbolic summation (cf. (ChyzakSalvy98, )). Implementations exist for Maple (for example, (AbramovLeLi2005, )) and Mathematica (for example, (KoutschanPhD, )).

We will only require a special case where the derivative is the zero map. These rings are sometimes also referred to as skew polynomials. For us, this is simply the (non-commutative) ring

[TABLE]

which consist of polynomial expressions in $\sigma$ with coefficients being polynomials in $t$. The addition in this ring is the same as for regular polynomials. The multiplication, however, is governed by the commutation rule

[TABLE]

(and consequently $\sigma t^{-1}=q^{-1}t^{-1}\sigma$ in the extensions which we introduce below). We will usually refer to the elements of $R$ as operators. These operators act on $\K(t)$ in the usual way: Denoting the action of $a\in\K[t,\sigma]$ on $y\in\K(t)$ by $a\bullet y$ we have $t\bullet y(t)=ty(t)$ and $\sigma\bullet y(t)=y(qt)$.

We will also need to consider a related ring where we allow rational functions in $t$ in the coefficients of our operators. The corresponding ring $\K(t)[\sigma]$ may be constructed from $\K[t,\sigma]$ since the non-zero polynomials in $t$ form a so-called (left or right) Ore set in $\K[t,\sigma]$ (for a definition see, for example, (cohnRT, , p. 177)). We have $\K[t,\sigma]\subseteq\K(t)[\sigma]$. The action of $\K[t,\sigma]$ on $\K(t)$ can easily be extended to an action of $\K(t)[\sigma]$.

We may now associate the system (rec’) with a matrix $A\in\K[t,\sigma]^{m\times m}$ writing it as

[TABLE]

Doing so, we may identify the transformations above as multiplication of the system by a certain matrix from the left. More precisely, each of the allowed transformations correspond to multiplication by a unimodular matrix $Q\in\K(t)[\sigma]^{m\times m}$, that is, a matrix which has an inverse $Q^{-1}\in\K(t)[\sigma]^{m\times m}$ in the same matrix ring. We denote the set of unimodular matrices by $\operatorname{GL}_{m}(\K(t)[\sigma])$. Note that $Q$ transforms the system (rec’) into

[TABLE]

where both the left and the right hand side are modified. It is easy to see that the solutions of the original and the transformed system are exactly the same if $Q$ is unimodular. We will call the two systems (and by extension also the matrices $A$ and $QA$) equivalent.

We would like to point out that being unimodular is a stronger property than merely having full (left row) rank where we understand the rank as the maximal number of linearly independent rows over $\K(t)[\sigma]$ – see also (BeckermannChengLabahn2006, , Def. 2.1, Thm. A.2). For instance, the $1$-by-$1$ matrix $(\sigma)\in\K(t)[\sigma]^{1\times 1}$ has full rank, but no inverse in $\K(t)[\sigma]^{1\times 1}$. We will call a matrix of full rank regular in order to distinguish this case.

Before we are able to propose an algorithm for transforming the system (rec’) into an equivalent system which is $t$-tail regular or $t$-head regular, we need to introduce one more concept. First, we would like to point out that the subring $\K[\sigma]$ of $\K[t,\sigma]$ may actually be regarded as a commutative polynomial ring because $\sigma$ commutes with all elements of $\K$. Therefore, it makes sense to look at matrix normal forms for polynomial matrices. We will again use the terms unimodular and regular, this time referring to matrices which have an inverse in $\K[\sigma]^{m\times m}$ or which have full rank over $\K[\sigma]$, respectively.

The normal form which we are going to employ here is the so-called Popov normal form. The most succinct way to describe it is using Gröbner bases for modules. See (myPhDthesis, ) for a detailed explanation; or see (Villard96, ) or Popov’s original paper (Popov72, ) for a more traditional definition. See, for example, (AdamsLoustaunau1994, ) for more details on Gröbner bases over modules. We introduce the term over position ordering on the polynomial module $\K[\sigma]^{1\times m}$ by saying that $e_{i}\sigma^{a}\leqslant_{\text{top}}e_{j}\sigma^{b}$ if $a<b$ or $a=b$ and $i\geqslant j$. Here, $e_{i}$ denotes the $i$th unit vector. We say that a matrix is in Popov normal form if its rows form a minimal Gröbner basis with respect to $\leqslant_{\text{top}}$. It is possible to prove that for every matrix $M\in\K[\sigma]^{m\times m}$ there exists a unimodular matrix $Q\in\operatorname{GL}_{m}(\K[\sigma])$ over $\K[\sigma]$ such that222That is, here unimodularity means that the inverse of $Q$ is a polynomial matrix $Q^{-1}\in\K[\sigma]^{m\times m}$.

[TABLE]

where $P\in\K[\sigma]^{\ell\times m}$ is in Popov normal form. It can further be shown that a matrix in Popov normal form has maximal row rank. Also, using Gröbner basis division (for each row of $C$) we can write every matrix $C\in\K[\sigma]^{\ell\times m}$ as $C=XP+Y$ where $X\in\K[\sigma]^{\ell\times\ell}$ and where $Y\in\K[\sigma]^{\ell\times m}$ has no rows which are reducible by the Popov normal form $P$ (again in the Gröbner basis sense). This offers an algorithmic way to determine whether $C$ is a left multiple of $P$. We are going to need this property in Algorithm 4.1 and Algorithm 4.5 below.

We would like to point out, that despite its definition in terms of Gröbner bases, the Popov normal form can easily be computed in polynomial time. In fact, even the naive method which is based on row reduction will be polynomial (cf. (myPhDthesis, , Lem. 5.14)). See, for example, (BeckermannLabahnVillard2006, , Cor. 6.1) for a faster method.

It is possible to use other row equivalent forms instead of the Popov normal form; as long as they have a similar division property. One possible example is the Hermite normal form which is a Gröbner basis with respect to the position over term ordering (see, for example, (myPhDthesis, ) for a details on why the Hermite normal form is a Gröbner basis). However, to our best knowledge, with current methods computing a Popov normal form can be done more efficiently than computing a Hermite normal form.

We briefly recall the reasoning of (myPhDthesis, , Lem 7.5) to explain how the Popov normal form can also be used to remove redundancies from a system: Starting with a general matrix $A\in\K[t,\sigma]^{m\times n}$ – that is, with a matrix which is not necessarily square or of full row rank – we first compute the column Popov normal form of $A$. It is defined analogously but for right modules of column vectors instead of left modules of row vectors. Moreover, algorithms to compute the (row) Popov normal form and the corresponding transformation matrix can easily be translated to the column Hermite normal form. We obtain a representation

[TABLE]

where $U\in\operatorname{GL}_{n}(\K(t)[\sigma])$ and $\tilde{A}$ has full column rank. Next, we compute the row Popov normal form of $AU$ which yields

[TABLE]

where $Q\in\operatorname{GL}_{m}(\K(t)[\sigma])$ and $\accentset{\approx}{A}$ has full row rank. Since the row transformations do not change the column rank and since (left) row and (right) column rank coincide by (cohn, , Thm. 8.1.1), we conclude that $\accentset{\approx}{A}$ must be square and of full row rank.

Now, as $Q$ and $U$ are both unimodular, the system $A\bullet x=b$ has a solution if and only if the system $\accentset{\approx}{A}\bullet y=c$ has a solution and $w=0$ where

[TABLE]

(with the blocks matching those of $QAU$). Thus, either the compatibility condition $w=0$ does not hold and we know that the system is not solvable; or it suffices to concentrate on the system $\accentset{\approx}{A}\bullet y=c$ which is of the correct form for the method presented in this paper. After solutions $y$ have been found, we can easily translate them into solutions $x$ of the original system using the matrix $Q$. (The vector $z$ does not have any conditions imposed on it and contributes free variables to the solution $x$.)

It is not necessary to employ the Popov normal form in order to remove the redundancies of the system. Any pair of rank revealing column and row transformation would be sufficient. For instance, we could use the Hermite normal form. Since we do not need the division property, also row and column reduction (see (BeckermannChengLabahn2006, ) for a definition and an efficient algorithm) or EG-eliminations (see (Abramov:99, )) could be used.

4. Desingularising the Determinants

With the notation and concepts from Section 3, we are finally prepared to deal with the case of identically vanishing determinants in Section 2. This section here will introduce an algorithm for transforming any system into an equivalent one (using elementary row transformations over $\K(t)[\sigma]$) where $\lambda$ or $\varrho$ do not vanish. This constitutes the main result of this paper.

We write system (rec’) in operator form $A\bullet y=t^{-\nu}b$ where $A\in\K[t,\sigma]^{m\times m}$ is an operator matrix as explained in Section 3. We express $A$ as

[TABLE]

where $\tilde{A}_{0},\ldots,\tilde{A}_{\ell}\in\K[\sigma]^{m\times m}$ are matrices in $\sigma$ alone. We will call $\tilde{A}_{0}$ the $t$-trailing matrix of $A$ and – assuming that $\tilde{A}_{\ell}\neq\mathbf{0}$ – we will call $\tilde{A}_{\ell}$ the $t$-leading matrix of $A$. It is easy to see that (using the notation of Section 2) we have

[TABLE]

In other words, $\lambda$ and $\varrho$ are the same as the determinants of $\tilde{A}_{0}$ and $\tilde{A}_{\ell}$ – except that we used $x$ instead of $\sigma$ as the name for the variable. (As explained in Section 3, $\K[\sigma]$ can be interpreted as univariate commutative polynomial ring.) Thus, the task of transforming the system into an equivalent one (using a unimodular multiplier) with $\lambda$ or $\varrho$ non-zero can be equivalently described as the task of having the $t$-trailing or the $t$-leading matrix being regular.

The definition of the $t$-trailing matrix in particular is also the reason why we allowed a denominator for the right hand side of (rec’). Should $\tilde{A}_{0}$ be the zero matrix, then we simply divide the entire equation by a suitable power of $t$ in order rectify that problem.

There are similar works which consider recurrence systems and place requirements on certain leading or trailing matrices. For instance, see the method in (AbramovKhmelnov2012, ) or the concept of strong row-reduction in (AbramovBarkatou2014, , Def. 4). There is, however, an essential difference between those approaches and the one presented here: We consider the leading and trailing matrices with respect to $t$ as the main variable while those other approaches consider $\sigma$ to be the main variable (and are consequently dealing with leading and trailing matrices that are in $\K(t)$).

The choice of variable does make a big difference. While we can easily do a non-commutative version of row-reduction (that is, working over $\K(t)[\sigma]$) in order to work on the leading coefficient with respect to $\sigma$, say; simply switching the variables would catapult us into $\K(\sigma)[t]$. The later ring, however, contains arbitrary denominators in $\sigma$ and therefore does not have a natural action on $\K(t)$ which extends the action of $\K[t,\sigma]$. (To see that let $p\in\K[\sigma]$ and $f\in\K(t)$. If we want to compute $p^{-1}\bullet f$, then we have to find $g\in\K(t)$ such that $p\bullet g=f$. However, not all equations of that form have rational solutions.)

We will thus need to develop a different approach which does not require division by expressions in $\sigma$.

Also note that our algorithm bears some similarities to EG-eliminations as described in (Abramov:99, ). In particular, the idea to reduce the trailing matrix $\tilde{A}_{0}$ and to shift down coefficients from the higher matrices whenever a row of $\tilde{A}_{0}$ becomes zero are the same. The difference is again the choice of the main variable. The transformations used during EG-eliminations are unimodular over $\K(t)[\sigma,\sigma^{-1}]$ (the Laurent skew polynomials333These are well-defined since the powers of $\sigma$ form an Ore set in $\K(t)[\sigma]$ – see, for example, (cohnRT, , Chptr. 5)) and do thus not alter the solutions; however switching the main variable does again expose us to arbitrary fractions in $\sigma$.

The need to keep all transformations unimodular over $\K(t)[\sigma]$ also means that we cannot proceed row by row as (Abramov:99, ) does where the $\sigma$-trailing or -leading matrix is brought into trapezoidal form one row at a time. Instead we always have to consider the entire $t$-trailing or -leading matrix for each elimination step. In particular, we cannot force the “$t$-width” (as defined in (Abramov:99, )) of the lower rows to decrease.

We are going to deal with the $t$-trailing matrix first. Below we state an algorithm which we claim will transform the system into the desired shape. Note that while we do not explicitly compute the transformation matrix $Q$, this could be easily done by applying all the transformations done to $A$ in parallel to the identity matrix $\mathbf{1}_{m}$. Alternatively, we may also apply them directly to the right hand side of (rec’).

Algorithm 4.1.

**Input:: **

A regular matrix $A\in\K[t,\sigma]^{m\times m}$.

**Output:: **

An equivalent matrix $B\in\K[t,\sigma]^{m\times m}$ such that the $t$-trailing matrix of $B$ is regular.

**Procedure:: **

1. **(1): **

Compute the $t$-trailing matrix $\tilde{A}_{0}$ of $A$. If $\det\tilde{A}_{0}\neq 0$, then return $A$. 2. **(2): **

Else, compute the Popov normal form of $\tilde{A}_{0}$ and let $X\in\operatorname{GL}_{m}(\K[\sigma])$ be the corresponding transformation matrix and let $r$ be the rank of $\tilde{A}_{0}$.444That is, the number of nonzero rows. Set $A\leftarrow XA$ and write the new trailing matrix of $A$ as

[TABLE]

where $\tilde{B}_{0}\in\K[\sigma]^{r\times m}$ is in Popov normal form. 3. **(3): **

Let $T=\operatorname{diag}(\mathbf{1}_{r},t^{-1}\mathbf{1}_{m-r})$ and set $A\leftarrow TA$.555This shifts the lower rows of $A$ by $t^{-1}$. The new trailing matrix is now

[TABLE]

with the same $\tilde{B}_{0}$ as before and some $\tilde{C}_{0}\in\K[\sigma]^{(m-r)\times m}$. 4. **(4): **

Use Gröbner bases reduction of $\tilde{C}_{0}$ by $\tilde{B}_{0}$ trying to eliminate all the rows in $\tilde{C}_{0}$. Write the result as $\tilde{C}_{0}=X\tilde{B}_{0}+Y$ with a matrix $X\in\K[\sigma]^{r\times r}$ and a matrix $Y\in\K[\sigma]^{(m-r)\times m}$. Let

[TABLE]

and set $A\leftarrow QA$. 5. **(5): **

If $Y=\mathbf{0}$, then continue the inner loop with Step 3. Else, go back to the outer loop in Step 1.

We would like to point out that for any matrix $X\in\K[\sigma]^{m\times m}$ we have $Xt=t\overline{X}$ for some $\overline{X}\in\K[\sigma]^{m\times m}$ since $\sigma$ is the $q$-shift. Thus, in Step 2 of Algorithm 4.1 the product $XA$ has the form $t^{\ell}X^{(\ell)}\tilde{A}_{\ell}+\ldots+tX^{(1)}\tilde{A}_{1}+X\tilde{A}_{0}$ for some matrices $X^{(1)},\ldots,X^{(\ell)}\in\K[\sigma]^{m\times m}$. In other words, $X$ acts on the $t$-trailing matrix of $A$ without interference from the other parts. Consequently, if after applying the transformation by $X$ the lower rows of the $t$-trailing matrix $\tilde{A}_{0}$ become zero, this means that the lower rows of the entire matrix $A$ must be divisible by $t$. This explains why the transformation by $T$ in Step 3 will still result in a matrix $A\in\K[t,\sigma]^{m\times m}$. A similar reasoning as in Step 2 also holds true for Step 4: again the transformation of $A$ by $Q$ acts on the $t$-trailing matrix without disturbance from the other coefficients. Its effect is to replace the lower block $\tilde{C}_{0}$ of $\tilde{A}_{0}$ by $Y$.

Note that this reasoning only works because of our special choice of $\sigma$ and would not be possible with more general automorphisms. In particular, the algorithm is not applicable for the normal shift case.

Theorem 4.2.

Algorithm 4.1* is correct and terminates.*

Before we can prove the theorem, we state two simple remarks:

Remark 4.3.

As long as the inner loop runs, we always consider the trailing matrix

[TABLE]

and compute a matrix $X$ (if possible) such that $C_{0}=XB_{0}$. (The inner loop terminates if that is not possible any more.) We then apply the following transformations to $A$

[TABLE]

Thus, we see that the transformation matrix for a single step has the form

[TABLE]

for some $\nu\geqslant 1$ (where we write $t^{-\nu}$ in the lower right block instead of $t^{-\nu}\mathbf{1}$). Moreover,

[TABLE]

for all $\mu,\nu\geqslant 0$ and matrices $X,\tilde{X}\in\K[t,\sigma]^{(m-r)\times r}$. Thus, we easily see that all transformation matrices in the inner loop must be of this shape.

Remark 4.4.

Let $A\in\K[t,\sigma]^{m\times m}$ have full (left) row rank. We want to prove that $A$ has a (two-sided) inverse in the quotient skew field $\K(t,\sigma)$. We first remark that we can embed $\K[t,\sigma]$ into the non-commutative Euclidean domain $\K(t)[\sigma]$. Over $\K(t)[\sigma]$, we can compute the Jacobson normal form

[TABLE]

of $A$ where $S,T\in\operatorname{GL}_{m}(\K(t)[\sigma])$ are unimodular (see, for example, (cohn, , Thm. 8.1.1)). The diagonal entries $e_{1},\ldots,e_{m}$ cannot be zero since otherwise $A$ did not have full row rank. Thus, the inverse $N^{-1}=\operatorname{diag}(e_{1}^{-1},\ldots,e_{m}^{-1})\in\K(t,\sigma)^{m\times m}$ of $N$ exists. We obtain

[TABLE]

Since $S^{-1}NT^{-1}$ is the product of (from both sides) invertible matrices, we conclude that $A$ is invertible.

Proof of Theorem 4.2.

If the algorithm terminates, then the trailing matrix $\tilde{A}_{0}$ is in Popov normal form which implies that it has full rank (over $\K[\sigma]$) which in turn implies that $\det\tilde{A}_{0}\neq 0$ as a polynomial in $\K[\sigma]$.

It remains to reason about the termination of the algorithm. First, we note that the outer loop starting in Step 1 cannot be reached infinitely often: A new Popov normal form is only computed if the lower rows of $\tilde{A}_{0}$ are not in the row space of the upper block. This implies that we either gain new nonzero rows in the Popov normal form or that the degrees or positions of the leading monomials decrease. In both cases, the module generated by the rows of $\tilde{A}_{0}$ becomes strictly larger which can happen only finitely often since $\K(t)[\sigma]^{1\times m}$ is noetherian.

Second, we have to show that the inner loop starting in Step 3 cannot be run infinitely often. For this, assume that the matrix

[TABLE]

with $B\in\K[t,\sigma]^{r\times m}$ and $C\in\K[t,\sigma]^{(m-r)\times m}$ has full (row) rank and assume that the inner loop repeats infinitely, that is, assume that for every $\nu\geqslant 0$ there are $C_{\nu}\in\K[t,\sigma]^{(m-r)\times m}$ and $X_{\nu}\in\K[t,\sigma]^{(m-r)\times r}$ such that

[TABLE]

(where we write $t^{-\nu}$ instead of $t^{-\nu}\mathbf{1}_{m-r}$ for the lower left block). We explained why the transformation matrices always look like this in Remark 4.3 above.

We now form the quotient (skew) field $\K(t,\sigma)$ of $\K[t,\sigma]$. Since $A$ has full row rank, it does have a (two-sided) inverse by Remark 4.4. Thus, the equation above is equivalent to the identity

[TABLE]

where we write $A^{-1}=\begin{pmatrix}P&Q\end{pmatrix}$ with $P\in\K(t,\sigma)^{m\times r}$ and $Q\in\K(t,\sigma)^{m\times(m-r)}$. In particular, we have the identity $C_{\nu}Q=t^{-\nu}.$ Note that $Q$ does not depend on $\nu$, that is, the denominator in $t$ of $Q$ is the same for every $\nu\geqslant 0$. In addition, $C_{\nu}$ is a polynomial matrix, that is, its denominator is always $1$. Thus, the denominator on the left hand side is bounded. However, the denominator of the right hand side is not. This is a contradiction. Hence, there cannot be infinitely many steps where $C_{\nu}$ is in the row space of $B$. ∎

An analogous method works for the $t$-leading matrix. The only differences to Algorithm 4.1 are that we have to work with the Popov normal form of the $t$-leading matrix $\tilde{A}_{\ell}$ in Step 3 and that we shift the lower rows by $t$ instead of $t^{-1}$ in Step 3. Again, we state the algorithm without explicit computation of the transformation matrix.

Algorithm 4.5.

**Input:: **

A regular matrix $A\in\K[t,\sigma]^{m\times m}$.

**Output:: **

An equivalent matrix $B\in\K[t,\sigma]^{m\times m}$ such that the $t$-leading matrix of $B$ is regular.

**Procedure:: **

1. **(1): **

Compute the $t$-leading matrix $\tilde{A}_{\ell}$ of $A$. If $\det\tilde{A}_{\ell}\neq 0$, then return $A$. 2. **(2): **

Else, compute the Popov normal form of $\tilde{A}_{\ell}$ and let $X\in\operatorname{GL}_{m}(\K[\sigma])$ be the corresponding transformation matrix and let $r$ be the rank of $\tilde{A}_{\ell}$. Set $A\leftarrow XA$ and write the new trailing matrix of $A$ as

[TABLE]

where $\tilde{B}_{\ell}\in\K[\sigma]^{r\times m}$ is in Popov normal form. 3. **(3): **

Let $T=\operatorname{diag}(\mathbf{1}_{r},t\mathbf{1}_{m-r})$ and set $A\leftarrow TA$.666This shifts the lower $m-r$ rows of $A$ by $t$. The new trailing matrix is now

[TABLE]

with the same $\tilde{B}_{\ell}$ as before and some matrix $\tilde{C}_{\ell}\in\K[\sigma]^{(m-r)\times m}$. 4. **(4): **

Use Gröbner bases reduction of $\tilde{C}_{\ell}$ by $\tilde{B}_{\ell}$ trying to eliminate all the rows in $\tilde{C}_{\ell}$. Write the result as $\tilde{C}_{\ell}=X\tilde{B}_{\ell}+Y$ with a matrix $X\in\K[\sigma]^{r\times r}$ and a matrix $Y\in\K[\sigma]^{(m-r)\times m}$. Let

[TABLE]

and set $A\leftarrow QA$. 5. **(5): **

If $Y=\mathbf{0}$, then continue the inner loop with Step 3. Else, go back to the outer loop in Step 1.

Theorem 4.6.

Algorithm 4.5* is correct and terminates.*

Proof.

The proof is mostly analogous to that of Theorem 4.2. The only noticeable change is that for the termination of the inner loop we have to deal with transformation matrices of the shape

[TABLE]

(with some $X_{\nu}\in\K[t,\sigma]^{(m-r)\times r}$) instead of having fractions in $t$. However, we come to an analogous equation $C_{\nu}Q=t^{\nu}$ where the degree of the left hand side is bounded while the degree of the right hand side is not. (It is easy to check that the degree in $t$ of $A$ during execution of Algorithm 4.5 is always equal to $\ell$). Thus, we arrive at a similar contradiction as for the proof of Theorem 4.2. ∎

Example 4.7.

We continue Example 2.6. The $t$-leading matrix of the system was

[TABLE]

which is not regular. Thus, we will apply Algorithm 4.5. The Popov normal form of the $t$-leading matrix is simply

[TABLE]

(Actually, for the proper Popov normal form we would need to divide the first row by $64$; however, we want to avoid fractions in order to save some space.) We obtain the new system matrix

[TABLE]

Next, we multiply the lower row by $t$ which leads to the new $t$-leading matrix

[TABLE]

which does not have full rank yet. Gröbner basis reduction amounts to adding the upper row to the lower. We obtain the new system matrix

[TABLE]

It turns out that we have to shift the lower row by $t$ three times before the $t$-leading matrix changes again (since the lower row is only linear in $t$). The final result is the system matrix

[TABLE]

with $t$-leading matrix

[TABLE]

This yields $\varrho=-64\sigma^{2}+768\sigma-2048$ with roots $4=q^{2}$ and $8=q^{3}$. Thus, the degree of polynomial solutions is at most $3$. From Example 2.3 we do know that the solution space is spanned by $(t^{3},1)$ and $(t^{2},1)$; so the bound is actually sharp in this case.

5. Conclusions

In this paper we have presented a method for determining degree bounds for the polynomial solutions and for determining the power of $t$ in the denominator bound for rational solutions of $q$-recurrence systems. Although the translation from the scalar case seems to be straight-forward at first, we quickly discovered a problem when certain determinants are vanishing. This required us to develop a new method for transforming the system into an equivalent form where those determinants are non-zero.

There exist other ways of obtaining the same information, as for instance uncoupling or transformation of the system into first order. However, those approaches are often computationally costly (since uncoupling is and since conversion to first order usually introduces a lot of new variables). The method presented in this paper does avoid those problems.

We do not have strict bounds on the complexity of the transformation algorithm yet; we hope to deliver those in the near future. However, in all the examples which we have computed, termination usually occurred within a very small number of steps. Moreover, the degree bounds found were usually tight. Although we did not yet carry out extensive comparisons with the existing methods mentioned above, we are therefore confident that this method will prove to be useful in practical applications.