# Decay estimates of Green's matrices for discrete-time linear periodic systems

A. Malyshev, Miloud Sadkane

► To cite this version:

A. Malyshev, Miloud Sadkane. Decay estimates of Green's matrices for discrete-time linear periodic systems. Linear and Multilinear Algebra, 2023, 72 (14), pp.2312-2328. 10.1080/03081087.2023.2256937. hal-04306777

HAL Id: hal-04306777

https://hal.univ-brest.fr/hal-04306777v1

Submitted on 11 Sep 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Decay estimates of Green's matrices for discrete-time linear periodic systems

A. N. Malyshev & M. Sadkane

**To cite this article:** A. N. Malyshev & M. Sadkane (2024) Decay estimates of Green's matrices for discrete-time linear periodic systems, Linear and Multilinear Algebra, 72:14, 2312-2328, DOI: 10.1080/03081087.2023.2256937

**To link to this article:** https://doi.org/10.1080/03081087.2023.2256937

<figure><img src="image_2.png" /></figure>

© 2023 The Author(s). Published by Informa UK Limited, trading as Taylor & Francis Group.

<figure><img src="image_3.png" /></figure>

Published online: 16 Sep 2023.

<figure><img src="image_4.png" /></figure>

Submit your article to this journal

<figure><img src="image_5.png" /></figure>

Article views: 229

<figure><img src="image_6.png" /></figure>

View related articles

<figure><img src="image_7.png" /></figure>

View Crossmark data

CrossMark
# Decay estimates of Green's matrices for discrete-time linear periodic systems

A. N. Malyshevª and M. Sadkaneᵇ

ªDepartment of Mathematics, University of Bergen, Bergen, Norway; ᵇUniv Brest, CNRS – UMR 6205, Laboratoire de Mathématiques de Bretagne Atlantique, Brest Cedex 3, France

## ABSTRACT

We study periodic Lyapunov matrix equations for a general discrete-time linear periodic system $B_p x_p - A_p x_{p-1} = f_p$, where the matrix coefficients $B_p$ and $A_p$ can be singular. The block coefficients of the inverse operator of the system are referred to as the Green matrices. We derive new decay estimates of the Green matrices in terms of the spectral norms of special solutions to the periodic Lyapunov matrix equations. The study is based on the periodic Schur decomposition of matrices.

## ARTICLE HISTORY

Received 3 June 2022

Accepted 3 July 2023

## COMMUNICATED BY

F. Poloni

## KEYWORDS

linear periodic system; Green matrix; discrete-time stability; periodic Schur decomposition; periodic Lyapunov matrix equation; periodic Sylvester matrix equation

# 1. Introduction

Linear periodic systems are frequently used mathematical models of periodic phenomena. A variety of such models is presented in the IFAC proceedings volume [1]. The book [2] provides fundamentals of modelling and analysis of both continuous and discrete-time linear periodic systems with emphasis on control theory. Modern algorithmic and computational issues related to linear periodic systems are discussed in papers [3–5].

We are concerned with the general discrete-time linear periodic systems here, namely with the doubly infinite linear systems of the form

$$B_p x_p - A_p x_{p-1} = f_p, \quad \forall p \in \mathbb{Z}, \qquad (1)$$

where $\mathbb{Z}$ is the set of integer numbers. The $N \times N$ complex matrices $B_p$ and $A_p$ are periodic in $p$ with period $P \ge 1$ so that $B_{p+P_k} = B_p$ and $A_{p+P_k} = A_p$ for all integers $k$. The matrices $A_p$ and $B_p$ are allowed to be singular and $f_p$ is a given complex $N$-vector. The doubly infinite vector sequences $x_p$ and $f_p$ are assumed to be elements of $\ell^2$, i.e. $\sum_{p=-\infty}^{\infty} x_p^* x_p < \infty$ and $\sum_{p=-\infty}^{\infty} f_p^* f_p < \infty$. By $V^*$ we denote the conjugate transpose of a matrix or vector $V$.

Let us first consider the time-invariant case where $B_p = I$ is the identity matrix and $A_p = A$ is constant for all $p \in \mathbb{Z}$, and the matrix $A$ is discrete-time stable, that is, all its eigenvalues lie in the open unit disk $\{z \in \mathbb{C} : |z| < 1\}$ of the complex plane $\mathbb{C}$. Then, the

## CONTACT

A. N. Malyshev

alexander.malyshev@uib.no

Postbox 7803, Bergen 5020, Norway

© 2023 The Author(s). Published by Informa UK Limited, trading as Taylor & Francis Group.

This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited. The terms on which this article has been published allow the posting of the Accepted Manuscript in a repository by the author(s) or with their consent.
<figure><img src="image_8.png" /></figure>

sensitivity of the stability of A to small perturbations can be characterized by the value $\min_{\phi \in \mathbb{R}} \sigma_{\min}(A - e^{i\phi}I)$, where $\sigma_{\min}$ denotes the minimum singular value of a matrix, $i = \sqrt{-1}$ and $\mathbb{R}$ is the set of real numbers. The sensitivity can be alternatively characterized by the parameter $1/\|H\|_2$, where $H = \sum_{n=0}^{\infty} (A^n)^* A^n$ is a positive definite solution of the discrete-time Lyapunov matrix equation $H - A^*HA = I$. Moreover, the magnitude $\|H\|_2$ allows us to derive the decay estimate (see, e.g. [6, chap. 9] or Remark 1.1)

$$ \|A^n\|_2 \le \sqrt{\|H\|_2} (1 - 1/\|H\|_2)^{n/2}, \quad n \ge 0. \qquad (2) $$

The interest of this estimate is not so much its sharpness (although it is very good for sufficiently large $n$), but rather the fact that it is easily implementable as it is based on standard Lyapunov matrix equations and eigenvalues of Hermitian positive-definite matrices. The literature on efficient algorithms dealing with these issues is indeed rich, see, e.g. [7, 8]. As we shall see later in this section and throughout the paper, our goal is to extend this estimate to system (1). This will be achieved again by means of some Lyapunov matrix equations, but unlike the time-invariant case, here the literature is still in its infancy when it comes to implementation issues, see the discussion at the end of this section.

The concept of matrix stability is naturally extended to the concept of spectral dichotomy of matrices and regular matrix pencils, or operators as in [9]. For example, a matrix A is said to possess the circular spectral dichotomy if it has no eigenvalues on the unit circle $\{z \in \mathbb{C} : |z| = 1\}$. A sensitivity theory of the circular dichotomy of matrices is given in [10]. Its extension to regular matrix pencils is developed in [11, 12]. The monograph [6] presents a unified theory of the spectral dichotomy of matrices and regular matrix pencils.

The spectral dichotomy theory developed in [6, 10–12] can be extended to discrete-time periodic systems. A standard way of extension is based on the so called lifting method [13], which transforms the $N \times N$ periodic system (1) to an augmented time-invariant system with matrix coefficients of size $NP \times NP$, for example of the form

$$ \mathcal{B}\mathcal{X}_n - \mathcal{A}\mathcal{X}_{n-1} = \mathcal{F}_n, \quad n \in \mathbb{Z}, \qquad (3) $$

where $\mathcal{X}_n = [x_{np}^T, x_{np+1}^T, \dots, x_{np+P-1}^T]^T$, $\mathcal{F}_n = [f_{np}^T, f_{np+1}^T, \dots, f_{np+P-1}^T]^T$ and

$$ \mathcal{A} = \begin{bmatrix} A_0 & & & & \\ & B_0 & & & \\ & -A_1 & B_1 & & \\ & & -A_2 & B_2 & \\ & & & \ddots & \\ & & & & -A_{P-1} & B_{P-1} \end{bmatrix}, \quad \mathcal{B} = \begin{bmatrix} B_0 & & & & \\ & B_1 & & & \\ & -A_1 & B_2 & & \\ & & -A_2 & B_3 & \\ & & & \ddots & \\ & & & & -A_{P-1} & B_{P-1} \end{bmatrix}. $$

Throughout this section we assume that the matrix $\mathcal{B} - e^{i\phi}\mathcal{A}$ is nonsingular for all $\phi \in \mathbb{R}$ (when this assumption fails, the system (3) may not be solvable, see Proposition 1.1). Consider the Fourier series

$$ (\mathcal{B} - e^{i\phi}\mathcal{A})^{-1} = \sum_{n=-\infty}^{\infty} \mathcal{G}_n e^{in\phi} \qquad (4) $$

with the matrix coefficients $\mathcal{G}_n$, which we refer to as the Green matrices for the system (3). This denomination is used, for example, in [10] for systems of type (3), where $\mathcal{B}$ is the identity matrix. It follows that the unique $l^2$-solution to (3) is represented in the convolutional
<figure><img src="image_9.png" /></figure>

form

$$ \mathcal{X}_n = \sum_{k=-\infty}^{\infty} \mathcal{G}_{n-k} \mathcal{F}_k, \quad n \in \mathbb{Z}. \qquad (5) $$

The matrix of the inverse operator determined by (5) is block Toeplitz and doubly infinite. The Toeplitz blocks are of size $NP \times NP$.

**Proposition 1.1:** If $\mathcal{B} - \mathcal{A} e^{i\phi_0}$ is singular for some $\phi_0 \in \mathbb{R}$, then system (2) has no solution $\mathcal{X}_n$ in $l^2$ for some $\mathcal{F}_n \in l^2$.

**Proof:** We may assume that the matrices $\mathcal{B}$ and $\mathcal{A}$ are upper triangular and that $\mathcal{B}_{NN} - \mathcal{A}_{NN} e^{i\phi_0} = 0$ by using the QZ factorization of the matrix pencil $\mathcal{B} - \mu \mathcal{A}$ with eigenvalue reordering [7]. Then the entries $y_n = (\mathcal{X}_n)_N$ and $g_n = (\mathcal{F}_n)_N$ satisfy the doubly infinite system

$$ \mathcal{B}_{NN} y_n - \mathcal{A}_{NN} y_{n-1} = g_n. \qquad (6) $$

When $\mathcal{B}_{NN} = \mathcal{A}_{NN} = 0$, the right-hand side $g_n$ must be 0 identically and this contradicts to existence of a solution for all sequences $g_n$. Now assume that $\mathcal{B}_{NN} \neq 0$ and $\mathcal{A}_{NN} \neq 0$. Then the solution of (6) with $g_0 = 1$ and $g_n = 0$ for all $n \neq 0$ has the components $y_n = x_0 e^{-in\phi_0}$ for $n \ge 0$ and $y_n = y_{-1} e^{-i(n+1)\phi_0}$ for $n < 0$. The sequence $y_n$ belongs to $l^2$ only if $y_0 = y_{-1} = 0$, i.e. $g_0$ must equal 0, which contradicts the assumption $g_0 = 1$. ■

Proposition 1.1 proves that the operator $\mathcal{L}$ in $l^2$ defined by (2) is bijective under the condition that $\mathcal{B} - \mathcal{A} e^{i\phi}$ is nonsingular for all $\phi \in \mathbb{R}$. Since the inverse operator $\mathcal{L}^{-1}$ is block convolutional as in (5), $\mathcal{G}_n \to 0$ as $|n| \to \infty$. The matrix equation $I = \sum_{n=-\infty}^{\infty} \mathcal{G}_n e^{in\phi} (\mathcal{B} - \mathcal{A} e^{i\phi})$ yields the following identities for the matrix Fourier coefficients $\mathcal{G}_n$:

$$ \mathcal{G}_n \mathcal{B} - \mathcal{G}_{n-1} \mathcal{A} = \delta_n I, \quad n \in \mathbb{Z}, \qquad (7) $$

where $\delta_n$ is the Kronecker symbol such that $\delta_0 = 1$ and $\delta_n = 0$ for $n \neq 0$. Similarly, $I = (\mathcal{B} - \mathcal{A} e^{i\phi}) \sum_{n=-\infty}^{\infty} \mathcal{G}_n e^{in\phi}$ yields $\mathcal{B}\mathcal{G}_n - \mathcal{A}\mathcal{G}_{n-1} = \delta_n I$, $n \in \mathbb{Z}$.

The following theorem shows that the decay rate of $\mathcal{G}_n$ as $|n| \to \infty$ can still be expressed analogously to (2), but the derivation is rather involved (see appendix for more detail).

**Theorem 1.2:**

$$ \|\mathcal{G}_n\|_2 \le \sqrt{\|\mathcal{H}^+\|_2} \left(1 - \frac{1}{\|\Pi_+ \mathcal{B}^* \mathcal{H}^+ \mathcal{B} \Pi_+\|_2}\right)^{n/2}, \quad n \ge 0, \qquad (8) $$

$$ \|\mathcal{G}_{-n}\|_2 \le \sqrt{\|\mathcal{H}^-\|_2} \left(1 - \frac{1}{\|\Pi_- \mathcal{A}^* \mathcal{H}^- \mathcal{A} \Pi_-\|_2}\right)^{(n-1)/2}, \quad n > 0, \qquad (9) $$

where $\Pi_+$ and $\Pi_-$ are the orthogonal projectors onto the right deflating subspaces of the pencil $\lambda \mathcal{B} - \mathcal{A}$ corresponding to the eigenvalues respectively inside and outside the unit circle, and $\mathcal{H}^+ = \sum_{n=0}^{\infty} \mathcal{G}_n^* \mathcal{G}_n$, $\mathcal{H}^- = \sum_{n=-\infty}^{-1} \mathcal{G}_n^* \mathcal{G}_n$ are Hermitian positive semidefinite matrices satisfying the discrete-time Lyapunov matrix equations

$$ \mathcal{B}^* \mathcal{H}^+ \mathcal{B} - \mathcal{A}^* \mathcal{H}^+ \mathcal{A} = (\mathcal{G}_0 \mathcal{B})^* (\mathcal{G}_0 \mathcal{B}), \qquad (10) $$
<figure><img src="image_10.png" /></figure>

$$A^*H^-A - B^*H^-B = (G_{-1}A)^*(G_{-1}A). \quad (11)$$

**Remark 1.1:** In the particular case when $P = 1$, $A_p = A$ and $B_p = I$ for all $p$ and $A$ is discrete-time stable, then $\Pi_+ = I$, $\Pi_- = 0$, $G_n = A^n$ if $n \ge 0$, $G_n = 0$ if $n < 0$, $H^+ = \sum_{n=0}^{\infty} (A^n)^*A^n$, $H^- = 0$, and the estimate (8) reduces to (2).

Given a sequence $\mathcal{F}_n$ such that $\mathcal{F}_n = 0$ whenever $|n| > n_K$ for some positive integer $n_K$, the corresponding solution $\mathcal{X}_n$ to (3) decays as $|n| \to \infty$. Owing to (5) the decay rate of $\mathcal{X}_n$ is bounded by the decay rate of $\mathcal{G}_n$ as $|n| \to \infty$.

A drawback of estimates (8) and (9) is that they involve the dense matrices $H^+$ and $H^-$ of large size $NP \times NP$ (which is the case, e.g. for small $N$ and large $P$ or for medium sizes $N$ and $P$). This drawback is overcome by means of special periodic Lyapunov matrix equations of size $N \times N$ as in [14]. However, this reference restricts the study to the case where $B_p = I$ and nonsingular $A_p$ and furthermore, the periodic Lyapunov matrix equations are constructed in a form which is not suitable for numerical computation. The present paper removes these restrictions and provides periodic Lyapunov matrix equations in the general case, where the matrices $B_p$ and $A_p$ can be singular. In addition, we give tighter decay estimates than those in [14].

More specifically, our approach, where periodicity plays an essential role, is based on the periodic Schur decomposition of matrices with reordering and solutions of generalized periodic Sylvester systems, which allow us to block diagonalize periodic matrix sequences according to the spectrum parts of the monodromy matrix inside and outside the unit circle. The periodic Lyapunov matrix equations are then written for the diagonal blocks thus avoiding projected right-hand sides. The non-projected periodic Lyapunov equations that is without projected right-hand side have advantage that they admit unique solutions, which can be found by efficient algorithms of numerical linear algebra; see e.g. [15]. These issues are developed in Sections 2, 3 and 4. Finally, in Sections 5 and 6, the desired decay estimates are written in terms of solutions to the non-projected periodic Lyapunov matrix equations.

The periodic Schur decomposition is a powerful tool for solving periodic Lyapunov matrix equations and periodic Riccati matrix equations; see e.g. [3, 4, 13, 16] and references therein. Periodic Lyapunov matrix equations in the context of balanced model reduction for periodic linear descriptor systems are treated in [17]. Note however that, so far, convergence failure issues in the periodic Schur decomposition limit its practical utility, but progress is being made; see [18].

The exponential decay of the entries of inverse of band matrices has been of some use in spline approximations as mentioned in [19]. Decay patterns of matrix inverses have also attracted considerable interest in linear system preconditioning, low-rank approximation strategies such as hierarchical matrices, wavelets etc; see [20, 21] and references therein. The exponential decay estimates of the present paper can be useful in the estimation of stability regions near periodic trajectories and in the Smith type methods [22, 23] for periodic matrix equations; see e.g. the analysis in [24] based on the decay bounds for the continuous periodic systems. The decay estimates also provide a quantitative convergence analysis of the so called doubling iterations for solving certain matrix equations [25].

Part IV of [26] contains other applications of exponential decay estimates and describes alternative decay estimates based on the Kreiss matrix theorem. However, the constants in
<figure><img src="image_11.png" /></figure>

the Kreiss matrix theorem are too expensive for computing [27]. Moreover, these estimates depend on the matrix size.

## 2. Green's matrices and periodic Lyapunov equations

The doubly infinite system (1) determines a bounded linear transform $\mathcal{L}: x = (x_p) \mapsto f = (f_p)$ in the space $l^2$ of square-summable doubly infinite vector sequences. We assume that $\mathcal{L}$ is invertible and represent its bounded inverse $\mathcal{L}^{-1}$ in the matrix form

$$x_l = \sum_{p \in \mathbb{Z}} G_{l,p} f_p, \quad l \in \mathbb{Z}, \qquad (12)$$

where we call the coefficient matrices $G_{l,p}$ by Green's matrices for (1) in $l^2$. In the control theory literature, the Green matrices are usually called the fundamental matrices.

Formula (5) determines the block structure of $G_{n-k}$:

$$G_{n-k} = \begin{bmatrix} G_{nP,kP} & G_{nP,kP+1} & \dots & G_{nP,kP+P-1} \\ G_{nP+1,kP} & G_{nP+1,kP+1} & \dots & G_{nP+1,kP+P-1} \\ \vdots & \vdots & \ddots & \vdots \\ G_{nP+P-1,kP} & G_{nP+P-1,kP+1} & \dots & G_{nP+P-1,kP+P-1} \end{bmatrix}. \qquad (13)$$

It follows from (13) that $G_{l,p}$ is $P$-periodic:

$$G_{l+Pm,p+Pm} = G_{l,p} \quad \text{for all } m \in \mathbb{Z}. \qquad (14)$$

The structure (13) also implies that the positive semidefinite Hermitian matrices

$$H_p^+ = \sum_{l=p}^{\infty} G_{l,p}^* G_{l,p}, \quad H_p^- = \sum_{l=-\infty}^{p-1} G_{l,p}^* G_{l,p}, \quad p \in \mathbb{Z} \qquad (15)$$

are the diagonal $N \times N$ blocks of the matrices $H^+$ and $H^-$ used in (10) and (11). The matrices $H_p^+$ and $H_p^-$ are periodic in $p$ with period $P$ on account of (14).

It follows from (12) that Green's matrices satisfy the identity

$$B_l G_{l,p} = A_l G_{l-1,p} + \delta_{l,p} I \qquad (16)$$

for all $l,p \in \mathbb{Z}$. The Kronecker symbol $\delta_{l,p}$ equals 1 if $l=p$ or 0 otherwise. The identity (16) is simply the elementwise expression of equation $\mathcal{L}\mathcal{L}^{-1} = \mathcal{I}$, where $\mathcal{I}$ stands for the identity operator in $l^2$. Inserting the equality $f_p = B_p x_p - A_p x_{p-1}$ into (12) gives the elementwise expression of equation $\mathcal{L}^{-1}\mathcal{L} = \mathcal{I}$:

$$G_{l,p} B_p = G_{l,p+1} A_{p+1} + \delta_{l,p} I. \qquad (17)$$

The matrices $H_p^+$ and $H_p^-$ are solutions to the following pair of periodic Lyapunov matrix equations

$$B_p^* H_p^+ B_p - A_{p+1}^* H_{p+1}^+ A_{p+1} = (G_{p,p} B_p)^* (G_{p,p} B_p), \qquad (18)$$
<figure><img src="image_12.png" /></figure>

$$A_{p+1}^* H_{p+1}^- A_{p+1} - B_p^* H_p^- B_p = (G_{p,p+1} A_{p+1})^* (G_{p,p+1} A_{p+1}). \quad (19)$$

These equations are direct implications of the identity (17). For instance, since $G_{l,p}B_p = G_{l,p+1}A_{p+1}$ for $l \neq p$,

$$A_{p+1}^* H_{p+1}^+ A_{p+1} = \sum_{l=p+1}^{\infty} (G_{l,p+1} A_{p+1})^* G_{l,p+1} A_{p+1} = \sum_{l=p+1}^{\infty} (G_{l,p} B_p)^* (G_{l,p} B_p) \\ = B_p^* H_p^+ B_p - (G_{p,p} B_p)^* (G_{p,p} B_p).$$

Equation (19) is derived analogously. We show in Section 4 that $G_{p,p}B_p$ and $-G_{p,p+1}A_{p+1} = I - G_{p,p}B_p$ are projectors for all $p \in \mathbb{Z}$; see Corollary 4.2.

Derivations in Section 4 show that the solution $H_p^+$ of (18) and solution $H_p^-$ of (19) are unique only under the additional conditions

$$(B_p G_{p,p})^* H_p^+ B_p G_{p,p} = H_p^+ \quad (20)$$

$$(A_p G_{p-1,p})^* H_p^- A_p G_{p-1,p} = H_p^- \quad (21)$$

Note that computation of $H_p^+$ and $H_p^-$ from the systems of equations (18), (20) and (19), (21) is complicated by the fact that they are overdetermined. Below we get rid of the projected right-hand sides in (18), (19) and conditions (20), (21).

### 3. Periodic Schur decomposition and block diagonalization

By a slight abuse of notation we write the monodromy matrix formally as

$$M = B_p^{-1} A_p B_{p-1}^{-1} A_{p-1} \dots B_1^{-1} A_1 \quad (22)$$

even if some matrices $B_p$ are singular. Such usage can be justified for our purposes by small perturbation of the matrices $B_p$. A strict but more involved consideration is based on the lifting method. If some $B_p$ are singular then the matrix (22) may have infinite eigenvalues and indefinite eigenvalues. The eigenvalue structure of $M$ is easily revealed by means of the periodic Schur decomposition, which is an extension of the QZ factorization of linear matrix pencils to periodic matrix pencils; see [28, 29] for more details.

**Theorem 3.1 (Periodic Schur decomposition):** There exist unitary matrices $Q_p, Z_p, p = 1, \dots, P$, and $Q_{P+1} = Q_1$ such that all the transformed matrices

$$\hat{B}_p = Z_p^* B_p Q_{p+1}, \quad \hat{A}_p = Z_p^* A_p Q_p, \quad p = 1, \dots, P,$$

are upper triangular.
<figure><img src="image_13.png" /></figure>

The diagonal entries of $\hat{B}_p$, $\hat{A}_p$ provide spectral information about system (1). By the aid of Theorem 3.1 the monodromy matrix (20) reads

$$M = Q_1 \hat{M} Q_1^*, \quad \text{where} \quad \hat{M} = \hat{B}_p^{-1} \hat{A}_p \hat{B}_{p-1}^{-1} \hat{A}_{p-1} \dots \hat{B}_1^{-1} \hat{A}_1. \qquad (23)$$

The eigenvalues $\lambda_i$ of $M$ are the diagonal entries of $\hat{M}$, namely

$$\lambda_i = a_i/b_i, \quad b_i = \prod_{p=1}^{P} (\hat{B}_p)_{ii}, \quad a_i = \prod_{p=1}^{P} (\hat{A}_p)_{ii}.$$

When $(a_i, b_i) = (0, 0)$, we say that $\lambda_i$ is indefinite. From now on, we restrict ourselves to the regular case, when the monodromy matrix $M$ has no indefinite eigenvalues. If $b_i \neq 0$ ($b_i = 0$), then $\lambda_i$ is a finite (infinite) eigenvalue of $M$.

We also assume below that the monodromy matrix $M$ in (22) has no eigenvalues on the unit circle $\{\lambda \in \mathbb{C} : |\lambda| = 1\}$. Otherwise, it is easy to construct a solution $x_p$ of (1) such that the solution $x_p \notin l^2$ and the corresponding right-hand side $f_p \in l^2$.

It is possible to select an arbitrary order of eigenvalues in the periodic Schur decomposition of Theorem 3.1; see [30]. Let the eigenvalues be ordered such that $\lambda_1, \dots, \lambda_{N_0}$ lie in the open unit disk, i.e. $|\lambda_i| < 1$ for $i = 1, \dots, N_0$, and $\lambda_{N_0+1}, \dots, \lambda_{N_0+N_\infty}$ lie outside the closed unit disk, so that $N_0 + N_\infty = N$ and $|\lambda_i| > 1$ for $i = N_0 + 1, \dots, N$.

Let us block partition the transformed matrices subject to the above ordering

$$\hat{B}_p = \begin{bmatrix} B_p^0 & B_p^{0\infty} \\ 0 & B_p^{\infty} \end{bmatrix}, \quad \hat{A}_p = \begin{bmatrix} A_p^0 & A_p^{0\infty} \\ 0 & A_p^{\infty} \end{bmatrix}. \qquad (24)$$

The matrices $\hat{B}_p$ and $\hat{A}_p$ are then simultaneously block-diagonalized as

$$\begin{bmatrix} I & -L_p \\ 0 & I \end{bmatrix} \hat{B}_p \begin{bmatrix} I & K_{p+1} \\ 0 & I \end{bmatrix} = \begin{bmatrix} B_p^0 & 0 \\ 0 & B_p^{\infty} \end{bmatrix} = \tilde{B}_p, \qquad (25)$$

$$\begin{bmatrix} I & -L_p \\ 0 & I \end{bmatrix} \hat{A}_p \begin{bmatrix} I & K_p \\ 0 & I \end{bmatrix} = \begin{bmatrix} A_p^0 & 0 \\ 0 & A_p^{\infty} \end{bmatrix} = \tilde{A}_p, \qquad (26)$$

where the blocks $K_p$ and $L_p$ satisfy the generalized P-periodic Sylvester system (an analog of the generalized Sylvester matrix equations; see, e.g. [31])

$$B_p^0 K_{p+1} - L_p B_p^{0\infty} = -B_p^{0\infty}, \quad A_p^0 K_p - L_p A_p^{0\infty} = -A_p^{0\infty}, \quad p = 1, \dots, P, \qquad (27)$$

$$K_{p+1} = K_1. \qquad (28)$$

Since $B_p^0$ and $A_p^\infty$ are nonsingular due to the ordering, the Sylvester matrix equations can be rewritten in the form $K_{p+1} = (B_p^0)^{-1} L_p B_p^\infty - (B_p^0)^{-1} B_p^{0\infty}$ and $L_p = A_p^0 K_p (A_p^\infty)^{-1} + A_p^{0\infty} (A_p^\infty)^{-1}$. After subsequent elimination of $L_p$ we obtain the periodic Lyapunov matrix system with respect to the set of matrices $K_p, p = 1, \dots, P$,

$$K_{p+1} = (B_p^0)^{-1} A_p^0 K_p (A_p^\infty)^{-1} B_p^\infty - F_p, \qquad (29)$$

where $K_{p+1} = K_1$, and $F_p = (B_p^0)^{-1} B_p^{0\infty} - (B_p^0)^{-1} A_p^{0\infty} (A_p^\infty)^{-1} B_p^\infty$.
<figure><img src="image_14.png" /></figure>

The system of P equations (29) can be further reduced to a single nonsymmetric discrete-time Lyapunov matrix equation

$$K_p - M_p^0 K_p M_{p-1}^\infty = C_p, \quad (30)$$

where

$$M_p^0 = (B_{p+P-1}^0)^{-1} A_{p+P-1}^0 \cdots (B_p^0)^{-1} A_p^0, \\ M_{p-1}^\infty = (A_p^\infty)^{-1} B_p^\infty \cdots (A_{p+P-1}^\infty)^{-1} B_{p+P-1}^\infty,$$

and

$$C_p = \sum_{j=1}^{P-1} \left( \left( \prod_{i=P-1}^{j} (B_{p+i}^0)^{-1} A_{p+i}^0 \right) F_{p+j-1} \prod_{i=j}^{P-1} (A_{p+i}^\infty)^{-1} B_{p+i}^\infty \right) - F_{p+P-1}.$$

Equation (30) has a unique solution because both matrices $M_p^0$ and $M_{p-1}^\infty$ are discrete-time stable; see, e.g. [6]. This leads to the following

**Theorem 3.2 (Periodic block diagonalization with eigenvalue ordering):** Assume that the monodromy matrix M has no indefinite eigenvalues and has no eigenvalues on the unit circle. Then, there exist nonsingular matrices $\tilde{Q}_p, \tilde{Z}_p, p = 1, \dots, P$, and $\tilde{Q}_{P+1} = \tilde{Q}_1$ such that

$$\tilde{B}_p = \tilde{Z}_p^{-1} B_p \tilde{Q}_{p+1}, \quad \tilde{A}_p = \tilde{Z}_p^{-1} A_p \tilde{Q}_p,$$

where all the transformed matrices $\tilde{B}_p, \tilde{A}_p, p = 1, \dots, P$, are upper triangular and block diagonal. The diagonal blocks $B_p^0, A_p^0, B_p^\infty, A_p^\infty$ of the transformed matrices

$$\tilde{B}_p = \begin{bmatrix} B_p^0 & 0 \\ 0 & B_p^\infty \end{bmatrix}, \quad \tilde{A}_p = \begin{bmatrix} A_p^0 & 0 \\ 0 & A_p^\infty \end{bmatrix}$$

are such that the eigenvalues of the monodromy matrices $M_p^0 = (B_{p+P-1}^0)^{-1} A_{p+P-1}^0 \cdots (B_p^0)^{-1} A_p^0$ and $M_p^\infty = (A_{p-P+1}^\infty)^{-1} B_{p-P+1}^\infty \cdots (A_p^\infty)^{-1} B_p^\infty$ lie inside the open unit disk for all integer p.

**Proof:** The matrices

$$\tilde{Q}_p = Q_p \begin{bmatrix} I & K_p \\ 0 & I \end{bmatrix}, \quad \tilde{Z}_p = Z_p \begin{bmatrix} I & L_p \\ 0 & I \end{bmatrix}$$

provide the desired block diagonalization; see (25), (26).

<figure><img src="image_15.png" /></figure>

## 4. Block diagonal representations

Recall that the matrix sequences $B_p, A_p, H_p^+, H_p^-$ are P-periodic in p. The Green matrices $G_{l,p}$ are also periodic as $G_{l,p} = G_{l+P_k, p+P_k}$ for all $k \in \mathbb{Z}$. If necessary, the matrices $Q_p, Z_p, K_p, L_p$ are extended P-periodically for all $p \in \mathbb{Z}$.
<figure><img src="image_16.png" /></figure>

The transformed Green matrices

$$ \tilde{G}_{l,p} = \tilde{Q}_{l+1}^{-1} G_{l,p} \tilde{Z}_p \quad (31) $$

satisfy the analogue of (7), namely

$$ \tilde{G}_{l,p} \tilde{B}_p = \tilde{G}_{l,p+1} \tilde{A}_{p+1} + \delta_{l,p} I \quad (32) $$

and also the similar equation $\tilde{B}_l \tilde{G}_{l,p} = \tilde{A}_l \tilde{G}_{l-1,p} + \delta_{l,p} I$.

**Theorem 4.1:** *The transformed Green matrices $\tilde{G}_{l,p}$ are block diagonal*

$$ \tilde{G}_{l,p} = \begin{bmatrix} G_{l,p}^0 & 0 \\ 0 & G_{l,p}^\infty \end{bmatrix} \quad (33) $$

with the blocks

$$ G_{l,p}^0 = \begin{cases} 0, & l < p, \\ (B_p^0)^{-1}, & l = p, \\ ((B_l^0)^{-1} A_l^0) \dots ((B_{p+1}^0)^{-1} A_{p+1}^0) (B_p^0)^{-1}, & l > p. \end{cases} \quad (34) $$

$$ G_{l,p}^{\infty} = \begin{cases} 0, & l > p-1, \\ -(A_p^{\infty})^{-1}, & l = p-1, \\ -((A_{l+1}^{\infty})^{-1} B_{l+1}^{\infty}) \dots ((A_{p-1}^{\infty})^{-1} B_{p-1}^{\infty}) (A_p^{\infty})^{-1}, & l < p-1. \end{cases} \quad (35) $$

**Proof:** It is easy to verify that matrices (33) with the blocks (34) and (35) satisfy Equation (32). Moreover,

$$ \lim_{l \to +\infty} \|G_{l,p}^0\| = 0, \quad \lim_{l \to -\infty} \|G_{l,p}^\infty\| = 0. \quad (36) $$

Convergence in (36) is linear. Therefore, the operator (12) constructed from the matrices (33) is bounded in $l^2$ and coincides with $L^{-1}$ owing to uniqueness of the inverse operator in $l^2$. ■

**Corollary 4.2:** *The identities $G_{p,p}^0 B_p^0 = I$ and $G_{p-1,p}^\infty A_p^\infty = -I$ are valid and*

$$ G_{p,p} B_p = \tilde{Q}_{p+1} \begin{bmatrix} I & 0 \end{bmatrix} \tilde{Q}_{p+1}^{-1}, \quad -G_{p-1,p} A_p = \tilde{Q}_p \begin{bmatrix} 0 & I \end{bmatrix} \tilde{Q}_p^{-1} $$

are projectors.

Let us denote the upper triangular factor of the Cholesky factorization of $I + L_p L_p^*$ by $(I + L_p L_p^*)^{1/2}$ and that of $I + K_p^* K_p$ by $(I + K_p^* K_p)^{1/2}$.

**Proposition 4.3:**

$$ \|G_{l,p}\|_2 = \|G_{l,p}^0 (I + L_p L_p^*)^{1/2}\|_2, \quad l \ge p, \quad (37) $$

$$ \|G_{l,p}\|_2 = \|(I + K_{l+1}^* K_{l+1})^{1/2} G_{l,p}^\infty\|_2, \quad l \le p-1. \quad (38) $$
<figure><img src="image_17.png" /></figure>

**Proof:** If $l \ge p$, then $G_{l,p}G_{l,p}^* = Q_{l+1} \begin{bmatrix} G_{l,p}^0(I+L_pL_p^*)(G_{l,p}^0)^* & 0 \\ 0 & 0 \end{bmatrix} Q_{l+1}^*$. If $l \le p-1$, then
$$G_{l,p}^*G_{l,p} = Z_p \begin{bmatrix} 0 & 0 \\ 0 & (G_{l,p}^*)^*(I+K_{l+1}^*K_{l+1})G_{l,p}^* \end{bmatrix} Z_p^*.$$

<figure><img src="image_18.png" /></figure>

The transformed matrices

$$\tilde{H}_p^+ = \tilde{Z}_p^* H_p^+ \tilde{Z}_p, \quad \tilde{H}_p^- = \tilde{Z}_p^* H_p^- \tilde{Z}_p. \qquad (39)$$

are block diagonal

$$\tilde{H}_p^+ = \begin{bmatrix} H_p^0 & 0 \\ 0 & 0 \end{bmatrix}, \quad \tilde{H}_p^- = \begin{bmatrix} 0 & 0 \\ 0 & H_p^\infty \end{bmatrix} \qquad (40)$$

with the diagonal blocks

$$H_p^0 = \sum_{l \ge p} (G_{l,p}^0)^* G_{l,p}^0 \quad \text{and} \quad H_p^\infty = \sum_{l \le p-1} (G_{l,p}^\infty)^* (I + K_{l+1}^* K_{l+1}) G_{l,p}^\infty. \qquad (41)$$

The P-periodic matrix sequences $H_p^0$ and $H_p^\infty$ satisfy the periodic discrete-time Lyapunov matrix systems

$$(B_p^0)^* H_p^0 B_p^0 - (A_{p+1}^0)^* H_{p+1}^0 A_{p+1}^0 = I, \qquad (42)$$

$$(A_{p+1}^{\infty})^* H_{p+1}^{\infty} A_{p+1}^{\infty} - (B_p^{\infty})^* H_p^{\infty} B_p^{\infty} = (I + K_{p+1}^{*} K_{p+1}). \qquad (43)$$

To estimate the norms in (37) and (38) we introduce the following matrices:

$$\tilde{B}_p^0 = (I + L_p L_p^*)^{-1/2} \tilde{B}_p^0, \quad \tilde{B}_p^\infty = B_p^\infty (I + K_{p+1}^* K_{p+1})^{-1/2}, \qquad (44)$$

$$\tilde{A}_p^0 = (I + L_p L_p^*)^{-1/2} \tilde{A}_p^0, \quad \tilde{A}_p^\infty = A_p^\infty (I + K_p^* K_p)^{-1/2}, \qquad (45)$$

$$\check{G}_{l,p}^0 = G_{l,p}^0 (I + L_p L_p^*)^{1/2}, \quad \check{G}_{l,p}^\infty = (I + K_{l+1}^* K_{l+1})^{1/2} G_{l,p}^\infty, \qquad (46)$$

$$\check{H}_p^0 = \sum_{l \ge p} (\check{G}_{l,p}^0)^* \check{G}_{l,p}^0, \quad \check{H}_p^\infty = \sum_{l \le p-1} (\check{G}_{l,p}^\infty)^* \check{G}_{l,p}^\infty = H_p^\infty. \qquad (47)$$

Note that the monodromy matrices $\check{M}_p^0 = (\check{B}_{p+p-1}^0)^{-1} \check{A}_{p+p-1}^0 \cdots (\check{B}_p^0)^{-1} \check{A}_p^0 = M_p^0$ and
$$\check{M}_p^\infty = (\check{A}_{p-p+1}^\infty)^{-1} \check{B}_{p-p+1}^\infty \cdots (\check{A}_p^\infty)^{-1} \check{B}_p^\infty = (I + K_{p+1}^* K_{p+1})^{1/2} M_p^\infty (I + K_{p+1}^* K_{p+1})^{-1/2}$$

are discrete-time stable. We also obtain $\|G_{l,p}\|_2 = \|\check{G}_{l,p}^0\|_2$ if $l \ge p$ and $\|G_{l,p}\|_2 = \|\check{G}_{l,p}^\infty\|_2$ if $l \le p-1$.

The matrices $\check{H}_p^0$ and $\check{H}_p^\infty$ satisfy the periodic discrete-time Lyapunov matrix systems

$$(\check{B}_p^0)^* \check{H}_p^0 \check{B}_p^0 - (\check{A}_{p+1}^0)^* \check{H}_{p+1}^0 \check{A}_{p+1}^0 = I \qquad (48)$$

$$(\check{A}_{p+1}^{\infty})^* \check{H}_{p+1}^{\infty} \check{A}_{p+1}^{\infty} - (\check{B}_p^{\infty})^* \check{H}_p^{\infty} \check{B}_p^{\infty} = I \qquad (49)$$

which are uniquely solvable. Moreover, all coefficient matrices $\check{B}_p^0$, $\check{A}_{p+1}^0$, $\check{B}_p^\infty$ and $\check{A}_{p+1}^\infty$ in systems (48) and (49) are upper triangular. Therefore, (48) and (49) are triangular linear systems when writing them with the Kronecker product of matrices.
<figure><img src="image_19.png" /></figure>

## 5. Decay estimates I

**Theorem 5.1:** The following decay estimates hold

$$ \|G_{l,p}^0\|_2 \le \sqrt{\|\check{H}_p^0\|_2} \prod_{i=p}^{l-1} \left(1 - 1/\|(\check{B}_i^0)^* \check{H}_i^0 \check{B}_i^0\|_2\right)^{\frac{1}{2}}, \quad l \ge p, \quad (50) $$

$$ \|G_{l,p}^{\infty}\|_2 \le \sqrt{\|\check{H}_p^{\infty}\|_2} \prod_{i=l+2}^{p} \left(1 - 1/\|(\check{A}_i^{\infty})^* \check{H}_i^{\infty} \check{A}_i^{\infty}\|\right)^{\frac{1}{2}}, \quad l \le p-1, \quad (51) $$

where the product over an empty set of indices is equal to 1.

**Proof:** The bound (50) is valid for $l = p$ as a consequence of the inequality $(\check{G}_{p,p}^0)^* \check{G}_{p,p}^0 \le \check{H}_p^0$. Further we need the variables $x_l^0 = \check{G}_{l,p}^0 f^0$ for $l \ge p$ such that $\check{B}_l^0 x_l^0 = \check{A}_l^0 x_{l-1}^0$ when $l > p$ and $\check{B}_p^0 x_p^0 = f^0$. If $l > p$, then owing to (48) and $\check{B}_l^0 x_l^0 = \check{A}_l^0 x_{l-1}^0$ we obtain the estimates

$$ \begin{aligned} (x_l^0)^*(\check{B}_l^0)^* \check{H}_l^0 \check{B}_l^0 x_l^0 &= (x_{l-1}^0)^*(\check{A}_l^0)^* \check{H}_l^0 \check{A}_l^0 x_{l-1}^0 \\ &= (x_{l-1}^0)^*(\check{B}_{l-1}^0)^* \check{H}_{l-1}^0 \check{B}_{l-1}^0 x_{l-1}^0 - (x_{l-1}^0)^* x_{l-1}^0 \\ &\le (1 - 1/\|(\check{B}_{l-1}^0)^* \check{H}_{l-1}^0 \check{B}_{l-1}^0\|_2) \cdot (x_{l-1}^0)^* (\check{B}_{l-1}^0)^* \check{H}_{l-1}^0 \check{B}_{l-1}^0 x_{l-1}^0 \\ &\le \left[ \prod_{i=p}^{l-1} (1 - 1/\|(\check{B}_i^0)^* \check{H}_i^0 \check{B}_i^0\|_2) \right] (x_p^0)^* (\check{B}_p^0)^* \check{H}_p^0 \check{B}_p^0 x_p^0. \end{aligned} $$

Corollary 4.2 and (44)-(46) imply that $\check{G}_{l,l}^0 \check{B}_l^0 = I$. Since

$$ (x_l^0)^*(x_l^0) = (x_l^0)^*(\check{B}_l^0)^* (\check{G}_{l,l}^0)^* \check{G}_{l,l}^0 \check{B}_l^0 x_l^0 \le (x_l^0)^*(\check{B}_l^0)^* \check{H}_l^0 \check{B}_l^0 x_l^0 $$

we arrive at the inequality

$$ \|\check{G}_{l,p}^0 f^0\|_2^2 \le \left[ \prod_{i=p}^{l-1} \left(1 - 1/\|(\check{B}_i^0)^* \check{H}_i^0 \check{B}_i^0\|_2\right) \right] (f^0)^* \check{H}_p^0 f^0, $$

which yields (50). The estimate (51) is derived similarly using (47), (49), the variables $x_l^\infty = \check{G}_{l,p}^\infty f^\infty$ for $l \le p-1$ satisfying $\check{B}_l^\infty x_l^\infty = \check{A}_l^\infty x_{l-1}^\infty$ when $l < p-1$ and $\check{A}_p^\infty x_{p-1}^\infty = -f^\infty$, the identities $G_{p,p+1}^\infty A_{p+1}^\infty = -I$ and

$$ \check{G}_{p,p+1}^{\infty} \check{A}_{p+1}^{\infty} = (I + K_{p+1}^{*} K_{p+1})^{1/2} G_{p,p+1}^{\infty} A_{p+1}^{\infty} (I + K_{p+1}^{*} K_{p+1})^{-1/2} = -I. $$

<figure><img src="image_20.png" /></figure>

**Remark 5.1:** Let $\Pi_p^+$ and $\Pi_p^-$ be the orthogonal projectors corresponding to the spectral projector pair $P_p = G_{p,p}B_p$ and $I - P_p = -G_{p,p+1}A_{p+1}$. It can be shown that for all $p = 1, 2, \dots, P$

$$ \begin{aligned} \|\check{H}_p^0\|_2 &= \|H_p^+\|_2, \quad \|\check{H}_p^\infty\|_2 = \|H_p^-\|_2, \\ \|(\check{B}_p^0)^* \check{H}_p^0 \check{B}_p^0\|_2 &= \|\Pi_p^+ B_p^* H_p^+ B_p \Pi_p^+\|_2, \quad \|(\check{A}_p^0)^* \check{H}_p^0 \check{A}_p^0\|_2 = \|\Pi_p^- A_p^* H_p^- A_p \Pi_p^-\|_2. \end{aligned} $$
<figure><img src="image_21.png" /></figure>

## 6. Decay estimates II

In this section we propose decay estimates which are based on the P-periodic sequences of the discrete-time stable monodromy matrices $M_p^0$ and $M_p^\infty$ from Theorem 3.2. The discrete-time Lyapunov matrix equations $\mathfrak{H}_p^0 - (M_p^0)^*\mathfrak{H}_p^0 M_p^0 = I$ and $\mathfrak{H}_p^\infty - (M_p^\infty)^*\mathfrak{H}_p^\infty M_p^\infty = I$ have unique solutions

$$ \mathfrak{H}_p^0 = \sum_{k \ge 0} ((M_p^0)^k)^k (M_p^0)^k, \quad \mathfrak{H}_p^\infty = \sum_{k \ge 0} ((M_p^\infty)^k)^k (M_p^\infty)^k, \qquad (52) $$

and (2) yields the following estimates for all integer $k \ge 0$:

$$ \|(M_p^0)^k\|_2 \le \sqrt{\|\mathfrak{H}_p^0\|_2 (1 - 1/\|\mathfrak{H}_p^0\|_2)^{k/2}}, \quad \|(M_p^\infty)^k\|_2 \le \sqrt{\|\mathfrak{H}_p^\infty\|_2 (1 - 1/\|\mathfrak{H}_p^\infty\|_2)^{k/2}}. \qquad (53) $$

**Theorem 6.1:** For all integer $k \ge 0$

$$ \|G_{l+Pk,p}^0\|_2 \le \sqrt{\|\mathfrak{H}_{l+1}^0\|_2 (1 - 1/\|\mathfrak{H}_{l+1}^0\|_2)^{k/2}} \|G_{l,p}^0\|_2, \quad l = p, p+1, \dots, p+P-1, \qquad (54) $$

$$ \|G_{l-Pk,p}^\infty\|_2 \le \sqrt{\|\mathfrak{H}_l^\infty\|_2 (1 - 1/\|\mathfrak{H}_l^\infty\|_2)^{k/2}} \|G_{l,p}^\infty\|_2, \quad l = p-1, \dots, p-P. \qquad (55) $$

*Proof:* It is easy to see that $G_{l+Pk,p}^0 = (M_{l+1}^0)^k G_{l,p}^0$ for $l \ge p$ and $G_{l-Pk,p}^\infty = (M_l^\infty)^k G_{l,p}^\infty$ for $l \le p-1$. The estimates (53) yield the desired estimates. ■

Alternatively, one can apply the decay estimates from [32] formulated in the following

**Theorem 6.2:** If A is a discrete-time stable matrix, then

$$ \|A^n\| \le \begin{cases} e(n+1)(1-d_A)^n, & n > (1-d_A)/d_A, \\ 1/d_A, & \text{otherwise,} \end{cases} \qquad (56) $$

where $d_A = \min_{\phi \in \mathbb{R}} \sigma_{\min}(e^{i\phi}I - A)$.

**Corollary 6.3:** Determine the parameters

$$ d_{M_p^0} = \min_{\phi \in \mathbb{R}} \sigma_{\min}(e^{i\phi}I - M_p^0) \quad \text{and} \quad d_{M_p^\infty} = \min_{\phi \in \mathbb{R}} \sigma_{\min}(e^{i\phi}I - M_p^\infty). $$

For all integer $k \ge (1 - d_{M_p^0})/d_{M_p^0}$,

$$ \|G_{l+Pk,p}^0\|_2 \le e(k+1)(1-d_{M_p^0})^k, \quad l=p, p+1, \dots, p+P-1, \qquad (57) $$

For all integer $k \ge (1 - d_{M_p^\infty})/d_{M_p^\infty}$,

$$ \|G_{l-Pk,p}^\infty\|_2 \le e(k+1)(1-d_{M_p^\infty})^k, \quad l=p-1, \dots, p-P. \qquad (58) $$
<figure><img src="image_22.png" /></figure>

**Example:** This simple illustrative example shows that the decay estimates of Theorem 6.1 can be better than the decay estimates of Theorem 5.1.

Let $P = 2$ and $\alpha$ be a very large real number, say, $\alpha = 2^{50}$. The coefficients of a stable 2-periodic linear system are $A_1 = \begin{pmatrix} 2^{-1} & 0 \end{pmatrix}$, $A_2 = \begin{pmatrix} 1 & \alpha \\ 0 & 1 \end{pmatrix}$, $B_1 = B_2 = I$. The monodromy matrices are $M_1^0 = A_2 A_1 = \begin{pmatrix} 2^{-1} & 0 \end{pmatrix}$ and $M_2^0 = A_1 A_2 = 2^{-1} \begin{pmatrix} 1 & \alpha \\ 0 & 0 \end{pmatrix}$. Green's matrices $G_{l,p}$ for $p = 0$ are $G_{0,0} = I$, $G_{2k-1,0} = G_{2k,0} = \begin{pmatrix} 2^{-k} & 0 \end{pmatrix}$, $k = 1, 2, \dots$, and $G_{l,0} = 0$ for $l < 0$.

Solution to the system (42) is $H_0 = \begin{pmatrix} 8/3 & 5\alpha/3 \\ 5\alpha/3 & 5\alpha^2/3 + 2 \end{pmatrix}$, $H_1 = \begin{pmatrix} 5/3 & 0 \\ 0 & 1 \end{pmatrix}$. Since $\|H_0\|_2 \approx 5\alpha^2/3$ and $\|H_1\|_2 = 5/3$, the bound (50) gives the weak estimate

$$ \begin{aligned} \|G_{2k,0}\|_2 &\le \sqrt{\|H_0\|_2} (1 - 1/\|H_0\|_2)^{k/2} (1 - 1/\|H_1\|_2)^{k/2} \\ &\approx \sqrt{\frac{5}{3}\alpha} \left(1 - \frac{3}{5\alpha^2}\right)^{k/2} \left(1 - \frac{3}{5}\right)^{k/2}. \end{aligned} $$

On the other hand, from (52) we obtain $\mathfrak{H}_1^0 = \begin{pmatrix} 4/3 & 0 \\ 0 & 1 \end{pmatrix}$, so $\|\mathfrak{H}_1^0\|_2 = 4/3$, and (54) yields the tight estimate

$$ \|G_{2k,0}\|_2 \le \sqrt{\|\mathfrak{H}_1^0\|_2} (1 - 1/\|\mathfrak{H}_1^0\|)^{k/2} = \frac{2}{\sqrt{3}} 2^{-k}. $$

We also have $d_{M_1^0} = 1/2$. Corollary 6.3 provides the weaker estimate $\|G_{2k,0}\|_2 \le e(k + 1)2^{-k}$ for $k > 1$.

When the matrices $A_p$ and $B_p$, $p = 1, 2$, are interchanged, we easily obtain an illustrative example with singular $B_p$. The Green matrices $G_{l,p}$ for $p = 0$ are $G_{-1,0} = -I$, $G_{-2k-1,0} = G_{-2k,0} = \begin{pmatrix} 2^{-k} & 0 \end{pmatrix}$, $k = 1, 2, \dots$, and $G_{l,0} = 0$ for $l \ge 0$.

Solution to the system (43) is $H_0 = \begin{pmatrix} 5/3 & 0 \\ 0 & 1 \end{pmatrix}$, $H_1 = \begin{pmatrix} 8/3 & 5\alpha/3 \\ 5\alpha/3 & 5\alpha^2/3 + 2 \end{pmatrix}$. Since $\|H_0\|_2 = 5/3$ and $\|H_1\|_2 \approx 5\alpha^2/3$, the bound (51) gives the estimate

$$ \begin{aligned} \|G_{-2k-1,0}\|_2 &\le \sqrt{\|H_0\|_2} (1 - 1/\|H_0\|_2)^{k/2} (1 - 1/\|H_1\|_2)^{k/2} \\ &\le \sqrt{\frac{5}{3}} \left(1 - \frac{3}{5\alpha^2}\right)^{k/2} \left(1 - \frac{3}{5}\right)^{k/2}. \end{aligned} $$

On the other hand, from (54) we obtain $\mathfrak{H}_1^\infty = \begin{pmatrix} 4/3 & 0 \\ 0 & 1 \end{pmatrix}$, so $\|\mathfrak{H}_1^\infty\|_2 = 4/3$, and (54) yields the estimate

$$ \|G_{-2k-1,0}\|_2 \le \sqrt{\|\mathfrak{H}_1^\infty\|_2} (1 - 1/\|\mathfrak{H}_1^\infty\|)^{k/2} = \frac{2}{\sqrt{3}} 2^{-k}. $$
<figure><img src="image_23.png" /></figure>

# Acknowledgments

The authors thank the anonymous referees for their criticisms and suggestions that helped to improve the presentation of this paper.

# Disclosure statement

No potential conflict of interest was reported by the author(s).

# Funding

This work was supported by the L. Meltzers Høyskolefond Research Fund.

# References

[1] Bittanti S, Colaneri P, editors. Periodic control systems. Amsterdam: Elsevier Science & Technology; 2002. (IFAC Proceedings; vol. 2001).

[2] Bittanti S, Colaneri P. Periodic systems: filtering and control. London: Springer; 2009.

[3] Varga A. Periodic Lyapunov equations: some applications and new algorithms. Int J Control. 1997;67:69–87. doi: 10.1080/002071797224360

[4] Varga A. Computational issues for linear periodic systems: paradigms, algorithms, open problems. Int J Control. 2013;86:1227–1239. doi: 10.1080/00207179.2013.773088

[5] Gusev S, Johansson S, Kågström B, et al. A numerical evaluation of solvers for the periodic Riccati differential equation. BIT Numer Math. 2010;50:301–329. doi: 10.1007/s10543-010-0257-5

[6] Godunov SK. Modern aspects of linear algebra. Providence (RI): American Math. Soc; 1998. (Translations of Mathematical Monographs, vol. 175).

[7] Golub GH, Van Loan CF. Matrix computations. 3rd ed., Baltimore (MD): The Johns Hopkins Univ Press; 1996.

[8] Datta BN. Numerical methods for linear control systems: design and analysis. San Diego (CA): Elsevier Academic Press; 2004.

[9] Daleckii JL, Krein MG. Stability of solutions to differential equations in Banach space. Providence (RI): American Math. Soc; 1974.

[10] Bulgakov AJ, Godunov SK. Circular dichotomy of a matrix spectrum. Sib Math J. 1988;29:59–70.

[11] Malyshev AN. Guaranteed accuracy in spectral problems of linear algebra I. Sib Adv Math. 1992;2(1):144–197.

[12] Malyshev AN. Guaranteed accuracy in spectral problems of linear algebra II. Sib Adv Math. 1992;2(2):153–204.

[13] Varga A, Van Dooren P. Computational methods for periodic systems - an overview. IFAC Proceedings Volumes (IFAC Papers-OnLine). 2001;34(12):167–172. doi: 10.1016/S1474-6670(17)34079-X

[14] Demidenko GV, Bondar LA. Exponential dichotomy of systems of linear difference equations with periodic coefficients. Sib Math J. 2016;57:969–980. doi: 10.1134/S0037446616060045

[15] Penzl T. Numerical solution of generalized Lyapunov equations. Adv Comput Math. 1998;8:33–48. doi: 10.1023/A:1018979826766

[16] Sreedhar J, Van Dooren P. Periodic Schur form and some matrix equations. In: Helmke U, Mennicken R, Saurer J, editors. Proceedings of Symposium on Mathematical Theory of Networks and Systems; Regensburg, Germany; 1993. vol. I. p. 339–362.

[17] Chu EKW, Fan HY, Lin WW. Projected generalized discrete-time periodic Lyapunov equations and balanced realization of periodic descriptor systems. SIAM J Matrix Anal Appl. 2007;29(3):982–1006. doi: 10.1137/040606715

[18] Sima V, Gahinet P. Using semi-implicit iterations in the periodic QZ algorithm. Proceedings of the 17th International Conference on Informatics in Control, Automation and Robotics; 2020, p. 35–46. doi: 10.5220/0009823300350046.
<figure><img src="image_24.png" /></figure>

[19] Demko S, Moss WF, Smith PW. Decay rates for inverses of band matrices. Math Comput. 1984;43(168):491-499. doi: 10.1090/mcom/1984-43-168

[20] Canuto C, Simoncini V, Verani M. On the decay of the inverse of matrices that are sum of Kronecker products. Linear Algebra Appl. 2014;452:21-39. doi: 10.1016/j.laa.2014.03.029

[21] Frommer A, Schimmel C, Schweitzer M. Non-Toeplitz decay bounds for inverses of Hermitean positive definite tridiagonal matrices. Electron Trans Numer Anal. 2018;48:362-372. doi: 10.1553/etna

[22] Smith R. Matrix equation XA + BX = C. SIAM J Appl Math. 1968;16:198-201. doi: 10.1137/0116017

[23] Davison EJ, Man FT. The numerical solution of $A'Q + QA = -C$. IEEE Trans Automat Control. 1968;13(4):448-449. doi: 10.1109/TAC.1968.1098954

[24] Klevtsova YY. A numerical algorithm for analyzing the asymptotic stability of solutions to linear systems with periodic coefficients. J Appl Industrial Math. 2009;3(2):234-245. doi: 10.1134/S1990478909020094

[25] Malyshev AN. Parallel algorithm for solving some spectral problems of linear algebra. Linear Algebra Appl. 1993;188/189:489-520. doi: 10.1016/0024-3795(93)90477-6

[26] Trefethen LN, Embree M. Spectra and pseudospectra: the behavior of nonnormal matrices and operators. Princeton (NJ): Princeton Univ. Press; 2005.

[27] Mitchell T. Fast interpolation-based globality certificates for computing Kreiss constants and the distance to uncontrollability. SIAM J Matrix Anal Appl. 2021;42(2):578-607. doi: 10.1137/20M1358955

[28] Bojanczyk A, Golub G, Van Dooren P. Periodic Schur decomposition: algorithms and applications. Proceedings of SPIE 1770, Advanced Signal Processing Algorithms, Architectures, and Implementations III; 1992. p. 31-42. doi: 10.1117/12.130915.

[29] Hench JJ, Laub AJ. Numerical solution of the discrete-time periodic Riccati equation. IEEE Trans Automat Control. 1994;39:1197-1210. doi: 10.1109/9.293179

[30] Granat R, Kågström B, Kressner D. Computing periodic deflating subspaces associated with a specified set of eigenvalues. BIT. 2007;47:763-791. doi: 10.1007/s10543-007-0143-y

[31] Stewart GW, Sun JG. Matrix perturbation theory. San Diego (CA): Academic Press; 1990.

[32] Bai Z, Demmel J, Gu M. An inverse free parallel spectral divide and conquer algorithm for nonsymmetric eigenproblems. Numer Math. 1997;76:279-308. doi: 10.1007/s002110050264

[33] Lancaster P, Rodman L. Algebraic Riccati equations. Oxford: Clarendon Press; 1995.

[34] Godunov SK. Problem of the dichotomy of the spectrum of a matrix. Sib Mat J. 1986;27:649-660. doi: 10.1007/BF00969193

[35] Bulgakov AY. Generalization of a matrix Lyapunov equation. Sib Math J. 1989;30:525-532. doi: 10.1007/BF00971750

[36] Stykel T. Analysis and numerical solution of generalized Lyapunov equations [PhD thesis]. Berlin (Germany): Technical University; 2002.

## Appendix. Proof of Theorem 1.2

The proof requires some intermediate results which we give in the form of lemmas in order to facilitate the reading.

**Lemma A.1: The Hermitian positive semidefinite matrices**

$$ \mathcal{H}^{+} = \sum_{n=0}^{\infty} \mathcal{G}_{n}^{*} \mathcal{G}_{n}, \quad \mathcal{H}^{-} = \sum_{n=-\infty}^{-1} \mathcal{G}_{n}^{*} \mathcal{G}_{n} \qquad (A1) $$

satisfy the discrete-time Lyapunov matrix equations

$$ \mathcal{B}^{*} \mathcal{H}^{+} \mathcal{B} - \mathcal{A}^{*} \mathcal{H}^{+} \mathcal{A} = (\mathcal{G}_{0} \mathcal{B})^{*} (\mathcal{G}_{0} \mathcal{B}), \qquad (A2) $$

$$ \mathcal{A}^{*} \mathcal{H}^{-} \mathcal{A} - \mathcal{B}^{*} \mathcal{H}^{-} \mathcal{B} = (\mathcal{G}_{-1} \mathcal{A})^{*} (\mathcal{G}_{-1} \mathcal{A}). \qquad (A3) $$
<figure><img src="image_25.png" /></figure>

**Proof:** The proof is based on Equation (7). For $n \ge 1$, $G_{n-1}A = G_nB$, and so

$$A^*H^+A = \sum_{n=1}^{\infty} (G_{n-1}A)^*G_{n-1}A = \sum_{n=1}^{\infty} (G_nB)^*G_nB = B^*H^+B - (G_0B)^*G_0B.$$

Similarly, for $n \le -1$, $G_nB = G_{n-1}A$, and so

$$B^*H^-B = \sum_{n=-\infty}^{-1} (G_nB)^*G_nB = \sum_{n=-\infty}^{-1} (G_{n-1}A)^*G_{n-1}A = A^*H^-A - (G_{-1}A)^*G_{-1}A.$$

<figure><img src="image_26.png" /></figure>

**Lemma A.2:** (1) The matrices $P_+ = G_0B$ and $P_- = I - P_+ = -G_{-1}A$ are the spectral projectors onto the right deflating subspaces of the pencil $\lambda B - A$ corresponding to the eigenvalues respectively inside and outside the unit circle.

(2) The orthogonal projectors $\Pi_+$ and $\Pi_-$ corresponding to $P_+$ and $P_-$ satisfy

$$P_+ = \Pi_+P_+, \quad P_- = \Pi_-P_- \quad (A4)$$

(3) The matrices $G_n$ satisfy the identities

$$G_0BG_n = G_n = G_nBG_0 \quad \text{if } n \ge 0, \quad (A5)$$

$$-G_{-1}AG_n = G_n = -G_nAG_{-1} \quad \text{if } n < 0. \quad (A6)$$

**Proof:** (1) The spectral projectors onto the right deflating subspace of the pencil $\lambda B - A$ corresponding to the eigenvalues inside the unit circle is given by (see, e.g. [33, chap. 1] or [6, chap. 10]) $P_+ = \frac{1}{2\pi i} \int_\gamma (zB - A)^{-1}B dz$, where $\gamma$ is any closed contour enclosing the eigenvalues in the open unit disk and excluding the other eigenvalues. Letting $\gamma$ be the unit circle $z = e^{-i\phi}$, $0 \le \phi < 2\pi$ leads to $P_+ = \frac{1}{2\pi} \int_0^{2\pi} (B - e^{i\phi}A)^{-1}B d\phi$.

On the other hand, it follows from (4) that $G_n = \frac{1}{2\pi} \int_0^{2\pi} (B - e^{i\phi}A)^{-1}e^{-in\phi} d\phi = \frac{1}{2\pi i} \int_\gamma (zB - A)^{-1}z^n dz$. Therefore $G_0B = P_+$, and equation (7) gives $P_- = I - G_0B = -G_{-1}A$.

(2) The orthogonal projectors associated with $P_+$ and $P_-$ are given by $\Pi_+ = P_+(P_+^*P_+ + P_+^*P_-)^{-1}P_+$ and $\Pi_- = P_-(P_+^*P_+ + P_+^*P_-)^{-1}P_-$.$\Pi_+ = P_+(P_+^*P_+ + P_+^*P_-)^{-1}P_-$ and hence $P_+ = \Pi_+P_+$. Similarly, we have $P_- = \Pi_-P_-$.

In fact, $P_+ = \tilde{P}_+P_+$ for any projector $\tilde{P}_+$ onto the same subspace as $P_+$, i.e. when range $(\tilde{P}_+)$ = range $(P_+)$. Indeed, for each vector $x$, $\tilde{P}_+(P_+x) = P_+x$. Analogously, $P_- = \tilde{P}_-P_-$ for any projector $\tilde{P}_-$ onto the same subspace as $P_-$.

(3) Let $\gamma_1$ and $\gamma_2$ be closed contours enclosing the eigenvalues in the open unit disk and excluding the other eigenvalues and suppose that $\gamma_1$ is inside $\gamma_2$. Then

$$\begin{align*}
G_0BG_n &= \frac{1}{(2\pi i)^2} \int_{\gamma_2} \int_{\gamma_1} (z_1B - A)^{-1}B(z_2B - A)^{-1}z_2^n dz_1 dz_2 \\
&= \frac{1}{(2\pi i)^2} \int_{\gamma_2} \int_{\gamma_1} \frac{(z_1B - A)^{-1} - (z_2B - A)^{-1}}{z_2 - z_1} z_2^n dz_1 dz_2 \\
&= \frac{1}{2\pi i} \int_{\gamma_1} \left( \frac{1}{2\pi i} \int_{\gamma_2} \frac{z_2^n}{z_2 - z_1} dz_2 \right) (z_1B - A)^{-1} dz_1 \\
&\quad - \frac{1}{2\pi i} \int_{\gamma_2} \left( \frac{1}{2\pi i} \int_{\gamma_1} \frac{dz_1}{z_2 - z_1} \right) (z_2B - A)^{-1} z_2^n dz_2 = G_n,
\end{align*}$$

where we have used the Cauchy formulas

$$\frac{1}{2\pi i} \int_{\gamma_2} \frac{z_2^n}{z_2 - z_1} dz_2 = z_1^n \quad \text{and} \quad \frac{1}{2\pi i} \int_{\gamma_1} \frac{dz_1}{z_2 - z_1} = 0.$$

The other equalities in (A5) and (A6) can be derived similarly.

<figure><img src="image_27.png" /></figure>
<figure><img src="image_28.png" /></figure>

**Proof of Theorem 1.2.:** Note first that the equalities $\mathcal{G}_n = \mathcal{G}_n \mathcal{B} \mathcal{G}_0$ for $n \ge 0$ and $\mathcal{G}_n = -\mathcal{G}_n \mathcal{A} \mathcal{G}_{-1}$ for $n < 0$ imply the equations

$$ \mathcal{H}^{+} = \mathcal{G}_{0}^{*} \mathcal{B}^{*} \mathcal{H}^{+} \mathcal{B} \mathcal{G}_{0}, \qquad (A7) $$

$$ \mathcal{H}^{-} = \mathcal{G}_{-1}^{*} \mathcal{A}^{*} \mathcal{H}^{-} \mathcal{A} \mathcal{G}_{-1}, \qquad (A8) $$

which should be added to the Lyapunov Equations (A2) and (A3) in order to provide uniqueness of their solutions $\mathcal{H}^+$ and $\mathcal{H}^-$. Now, owing to (A2), the equalities

$$ \mathcal{G}_{n}^{*}(\mathcal{B}^{*}\mathcal{H}^{+}\mathcal{B})\mathcal{G}_{n} = \mathcal{G}_{n-1}^{*}\mathcal{A}^{*}\mathcal{H}^{+}\mathcal{A}\mathcal{G}_{n-1} = \mathcal{G}_{n-1}^{*}[\mathcal{B}^{*}\mathcal{H}^{+}\mathcal{B} - (\mathcal{G}_{0}\mathcal{B})^{*}\mathcal{G}_{0}\mathcal{B}]\mathcal{G}_{n-1} $$

hold if $n \ge 1$. In view of (A4), $\mathcal{G}_0 \mathcal{B} = \Pi_+ \mathcal{G}_0 \mathcal{B}$. As a consequence of (A7)

$$ \begin{aligned} \mathcal{G}_{n-1}^{*} \mathcal{B}^{*} \mathcal{H}^{+} \mathcal{B} \mathcal{G}_{n-1} &= \mathcal{G}_{n-1}^{*} \mathcal{B}^{*} \mathcal{G}_{0}^{*} \mathcal{B}^{*} \mathcal{H}^{+} \mathcal{B} \mathcal{G}_{0} \mathcal{B} \mathcal{G}_{n-1} \\ &= \mathcal{G}_{n-1}^{*} \mathcal{B}^{*} \mathcal{G}_{0}^{*} \Pi_{+} \mathcal{B}^{*} \mathcal{H}^{+} \mathcal{B} \Pi_{+} \mathcal{G}_{0} \mathcal{B} \mathcal{G}_{n-1} \\ &\le \| \Pi_{+} \mathcal{B}^{*} \mathcal{H}^{+} \mathcal{B} \Pi_{+} \|_{2} (\mathcal{G}_{0} \mathcal{B} \mathcal{G}_{n-1})^{*} \mathcal{G}_{0} \mathcal{B} \mathcal{G}_{n-1} \end{aligned} $$

and

$$ \mathcal{G}_{n}^{*}(\mathcal{B}^{*}\mathcal{H}^{+}\mathcal{B})\mathcal{G}_{n} \le (1 - 1/\|\Pi_{+}\mathcal{B}^{*}\mathcal{H}^{+}\mathcal{B}\Pi_{+}\|_{2})\mathcal{G}_{n-1}^{*}(\mathcal{B}^{*}\mathcal{H}^{+}\mathcal{B})\mathcal{G}_{n-1}. $$

Hence

$$ \mathcal{G}_{n}^{*}(\mathcal{B}^{*}\mathcal{H}^{+}\mathcal{B})\mathcal{G}_{n} \le (1 - 1/\|\Pi_{+}\mathcal{B}^{*}\mathcal{H}^{+}\mathcal{B}\Pi_{+}\|_{2})^{n}\mathcal{G}_{0}^{*}(\mathcal{B}^{*}\mathcal{H}^{+}\mathcal{B})\mathcal{G}_{0}, $$

and using again (A7) leads to the estimate

$$ \mathcal{G}_{n}^{*}(\mathcal{B}^{*}\mathcal{H}^{+}\mathcal{B})\mathcal{G}_{n} \le (1 - 1/\|\Pi_{+}\mathcal{B}^{*}\mathcal{H}^{+}\mathcal{B}\Pi_{+}\|_{2})^{n}\mathcal{H}^{+}. $$

Since $\mathcal{G}_0 \mathcal{B} \mathcal{G}_n = \mathcal{G}_n$, we can derive the estimate

$$ \mathcal{G}_{n}^{*}(\mathcal{B}^{*}\mathcal{H}^{+}\mathcal{B})\mathcal{G}_{n} \ge \mathcal{G}_{n}^{*}\mathcal{B}^{*}\mathcal{G}_{0}^{*}\mathcal{G}_{0}\mathcal{B}\mathcal{G}_{n} = \mathcal{G}_{n}^{*}\mathcal{G}_{n}. \qquad (A9) $$

The estimate (8) of Theorem 1.2 then follows from the inequality $\mathcal{G}_n^* \mathcal{G}_n \le (1 - 1/\|\Pi_+ \mathcal{B}^* \mathcal{H}^+ \mathcal{B} \Pi_+ \|_2)^n \mathcal{H}^+$. The estimate (9) is derived analogously. ■

## Comments

(1) The Lyapunov equations (A2) and (A3) are often called generalized Lyapunov equations [6, chap. 10]) because the matrices $P_+ = \mathcal{G}_0 \mathcal{B}$ and $P_- = -\mathcal{G}_{-1} \mathcal{A}$ are projectors. The generalized Lyapunov equations were first introduced by Godunov for continuous-time linear systems in [34]. Since the equations from [34] have nonunique solutions, Bulgakov supplemented them by additional equations in [35] in order to get a matrix system, which has a unique solution. Generalized Lyapunov equations with unique solution in the general discrete-time case, including singular $\mathcal{B}$, were derived by Malyshev in [11, 12]. Generalized Lyapunov equations for descriptor continuous-time case were introduced by Stykel in [36].

(2) To our best knowledge, the earliest publication which contains decay estimates of the form (2) for the continuous-time Lyapunov matrix equations, is [9, proof of Theorem 5.1]. Decay estimates for discrete-time linear systems $\mathcal{B} \mathcal{X}_n - \mathcal{A} \mathcal{X}_{n-1} = \mathcal{F}_n$ with an arbitrary matrix $\mathcal{B}$ first appeared in [11, 12]. Decay estimates for descriptor continuous-time linear systems are derived in [36].

(3) In previous publications, the denominators in (8) and (9) are $\|\mathcal{B}^* \mathcal{H}^+ \mathcal{B}\|_2$ instead of $\|\Pi_+ \mathcal{B}^* \mathcal{H}^+ \mathcal{B} \Pi_+\|_2$ and $\|\mathcal{A}^* \mathcal{H}^- \mathcal{A}\|_2$ instead of $\|\Pi_- \mathcal{A}^* \mathcal{H}^- \mathcal{A} \Pi_-\|_2$. The variant in (8) and (9) gives tighter estimates, which, as far as we know, are not available in the literature.