Solvability Theory and Iteration Method for One Self-Adjoint Polynomial Matrix Equation

Abstract

The solvability theory of an important self-adjoint polynomial matrix equation is presented, including the boundary of its Hermitian positive definite (HPD) solution and some sufficient conditions under which the (unique or maximal) HPD solution exists. The algebraic perturbation analysis is also given with respect to the perturbation of coefficient matrices. An efficient general iterative algorithm for the maximal or unique HPD solution is designed and tested by numerical experiments.

1. Introduction

In this paper, we consider the following self-adjoint polynomial matrix equation: 𝑋^𝑠−𝐴^∗𝑋^𝑡𝐴=𝑄, (1) where 𝑠,𝑡 are positive integers, 𝐴,𝑄∈ℂ^𝑛×𝑛, and 𝑄>0. As far as we know, the solvability of (1) is not completely solved untill now.

In many fields of applied mathematics, engineering, and economic sciences, (1) plays an important role. The famous discrete-time algebraic Lyapunov equation (DALE) is exactly (1) with 𝑠=𝑡=1. Undoubtedly, DALE is one of the most important mathematical problems in signal processing, system, and control theory and many others (e.g., see the monographs [1, 2]). If 𝐴 is stable (with respect to the unit circle), DALE has a unique Hermitian positive definite (HPD) solution. Such strong relation between the spectral property of 𝐴 and the solvability theory is fortunately owned by (1), which can be considered as a nonlinear DALE if 𝑠≠1 or 𝑡≠1. What about the following algebraic Riccati equation: 𝑌²+𝐵^∗𝑌+𝑌𝐵−𝐴^∗𝑌𝐴−𝑅=0, (2) where 𝐴,𝐵,𝑅∈ℂ^𝑛×𝑛, 𝐵^∗=𝐵≥0, and 𝑅^∗=𝑅>0? Defining 𝑋∶=𝑌+𝐵 and 𝑄∶=𝑅+𝐵²−𝐴^∗𝐵𝐴, we can immediately get (1) with 𝑠=2 and 𝑡=1 as an equivalent form of (2). As we all know, solving algebraic Riccati equations is an important task in the linear-quadratic regulator problem, Kalman filtering, 𝐻_∞-control, model reduction problems, and so forth. See [1, 3–5] and the references therein. Many numerical methods have been proposed, such as invariant subspace methods [6], Schur method [7], doubling algorithm [8], and structure-preserving doubling algorithm [9, 10]. At the same time the perturbation theory was developed in [11–15], as well as the unified methods for the discrete-time and continuous-time algebraic Riccati equations [16, 17]. A general iteration method for (1) given in this paper can be seen as a new algorithm for the algebraic Riccati equation (2), setting 𝑠=2 and 𝑡=1.

Apart from the above applications, (1) is appealing from the mathematical viewpoint since it unifies a large class of systems of polynomial matrix equations. Many nonlinear matrix equations are special cases of (1). For example, nonlinear matrix equations, 𝑋−𝐴^∗𝑋^𝑞𝐴=𝑄 (see, e.g., [18, 19]), are equivalence models of 𝑌^𝑠−𝐴^∗𝑌^𝑡𝐴=𝑄 and 𝑌=𝑋^1/𝑠, where 𝑠,𝑡 are positive integers and 𝑞=𝑡/𝑠. In a rather general form, Ran and Reurings [18] investigated 𝑋+𝐴^∗ℱ(𝑋)𝐴=𝑄 (𝑄>0) for its positive semidefinite solutions under the assumption that the function ℱ(⋅) is monotone and 𝑄−𝐴^∗ℱ(𝑄)𝐴 is positively definite. Besides, Lee and Lim [20] proved that (1) has a unique HPD solution when |𝑠|≥1≥|𝑡| and |𝑡/𝑠|<1. See [21–25] for more recent results on nonlinear matrix equations. To the best of our knowledge, (1) with 𝑠<𝑡 (without monotony in hand) has not been discussed. These facts motivate us to study polynomial matrix equation (1).

This paper is organized as follows. In Section 2 we deduce the existence and uniqueness conditions of HPD solutions of (1); in Section 3 we derive the algebraic perturbation theory for the unique or maximal solution of (1); finally in Section 4, we provide an iterative algorithm and two numerical experiments.

We begin with some notations used throughout this paper. 𝔽^𝑚×𝑛 stands for the set of 𝑚×𝑛 matrices with elements on field 𝔽(𝔽isℝorℂ). If 𝐻 is a Hermitian matrix on 𝔽^𝑛×𝑛, 𝜆_min(𝐻) and 𝜆_max(𝐻) stand for the minimal and the maximal eigenvalues, respectively. Denote the singular values of a matrix 𝐴∈𝔽^𝑚×𝑛 by 𝜎₁(𝐴)≥⋯≥𝜎_𝑙(𝐴)≥0, where 𝑙=min{𝑚,𝑛}. Suppose that 𝑋 and 𝑌 are Hermitian matrices; we write 𝑋≥𝑌(𝑋>𝑌) if 𝑋−𝑌 is positively semidefinite (definite) and denote the matrices set {𝑋∣𝑋−𝛼𝐼≥0and𝛽𝐼−𝑋≥0} by [𝛼𝐼,𝛽𝐼].

2. Solvability of Self-Adjoint Polynomial Matrix Equation

In this section, we study the solvability theory of (1) assuming that 𝐴 is nonsingular; that is, 𝜆_min(𝐴^∗𝐴)>0. To do this, we need two simple but useful functions defined on the positive abscissa axis: 𝑔₁(𝑥)=𝑥^𝑠−𝜆_max𝐴^∗𝐴𝑥^𝑡−𝜆_max(𝑄),𝑔₂(𝑥)=𝑥^𝑠−𝜆_min𝐴^∗𝐴𝑥^𝑡−𝜆_min(𝑄). (3)

The following two famous inequalities will be used frequently in the remaining of this paper.

Lemma 1 (Löwner-Heinz inequality [26, Theorem 1.1]). If 𝐴≥𝐵≥0 and 0≤𝑟≤1, then 𝐴^𝑟≥𝐵^𝑟.

Lemma 2 (see [27, Theorem 2.1]). Let 𝐴 and 𝐵 be positive operators on a Hilbert space 𝐻, such that 𝑀₁𝐼≥𝐴≥𝑚₁𝐼>0,𝑀₂𝐼≥𝐵≥𝑚₂𝐼>0, and 0<𝐴≤𝐵. Then 𝐴^𝑡≤𝑀₁𝑚₁^𝑡−1𝐵^𝑡,𝐴^𝑡≤𝑀₂𝑚₂^𝑡−1𝐵^𝑡 (4) hold for any 𝑡≥1.

2.1. Maximal Solution of (1) with 𝑠<𝑡

Now we derive a necessary condition and a sufficient condition for existence of HPD solutions of (1) with 𝑠<𝑡. With 𝑔₁(𝑥) and 𝑔₂(𝑥) in hand, we can easily get the distribution of eigenvalues of the HPD solution 𝑋 of (1).

Theorem 3. Suppose that 𝜆_max(𝐴^∗𝐴)≤(𝑠/𝑡)((𝑡−𝑠)/𝜆_max(𝑄)𝑡)^{(𝑡−𝑠)/𝑠} and 𝑋∈ℂ^𝑛×𝑛 is an HPD solution of (1); then for any eigenvalue 𝜆(𝑋) of 𝑋, 𝛽₁≤𝜆(𝑋)≤𝛼₁𝑜𝑟𝛼₂≤𝜆(𝑋)≤𝛽₂, (5) where 𝛼₁,𝛼₂ are two positive roots of 𝑔₁(𝑥) and 𝛽₁,𝛽₂ are two positive roots of 𝑔₂(𝑥).

Proof. From Theorem 3.3.16(d) in Horn and Johnson [28], one can see that 𝜎_𝑖𝐴^∗𝑋^𝑡𝐴≤𝜎_𝑖𝑋^𝑡𝜎²₁(𝐴)𝜆,thatis,_𝑖𝐴^∗𝑋^𝑡𝐴≤𝜆_𝑖𝑋^𝑡𝜆_max𝐴^∗𝐴,𝑖=1,…,𝑛. (6)

If 𝐴 is nonsingular, 𝜎_𝑖𝑋^𝑡=𝜎_𝑖𝐴^−1∗𝐴^∗𝑋^𝑡𝐴𝐴⁻¹≤𝜎_𝑛⁻²(𝐴)𝜎_𝑖𝐴^∗𝑋^𝑡𝐴. (7) That means 𝜎_𝑖𝐴^∗𝑋^𝑡𝐴≥𝜎_𝑖𝑋^𝑡𝜎²_𝑛(𝐴)𝜆,thatis,_𝑖𝐴^∗𝑋^𝑡𝐴≥𝜆_𝑖𝑋^𝑡𝜆_min𝐴^∗𝐴,𝑖=1,…,𝑛. (8) The above equations still hold if 𝐴 is singular, since 𝜎_𝑛(𝐴)=0, that is, 𝜆_min(𝐴^∗𝐴)=0, in this case. Applying Weyl theorem in Horn and Johnson [29], 𝑋^𝑠=𝑄+𝐴^∗𝑋^𝑡𝐴 implies 𝜆(𝑋)^𝑠−𝜆_max𝐴^∗𝐴𝜆(𝑋)^𝑡𝜆_max(𝑄)𝜆(𝑋)≤0,_𝑠−𝜆_min𝐴^∗𝐴𝜆(𝑋)^𝑡−𝜆_min(𝑄)≥0. (9)
Define a function 𝑓(𝑥)=𝑥^𝑠−𝑎²𝑥^𝑡−𝑞,𝑎>0,𝑞>0. Then the only positive stationary point of 𝑓(𝑥) is 𝑥₀=((𝑡/𝑠)𝑎²)^{1/(𝑠−𝑡)}. If 𝑎²≤(𝑠/𝑡)((𝑡−𝑠)/𝑞𝑡)^{(𝑡−𝑠)/𝑠}, 𝑓(𝑥) has two positive roots, 𝑥₁ and 𝑥₂, with 𝑞^1/𝑠<𝑥₁≤𝑥₀≤𝑥₂<𝑎^{2/(𝑠−𝑡)}. So 𝜆_max(𝐴^∗𝐴)≤(𝑠/𝑡)(𝑡−𝑠)/𝜆_max(𝑄)𝑡^{(𝑡−𝑠)/𝑠} implies that 𝑔₁(𝑥) has two roots 𝛼₁,𝛼₂>0 and 𝑔₂(𝑥) has two roots 𝛽₁,𝛽₂>0. Since 𝑔₂(𝑥)≥𝑔₁(𝑥), (𝜆_min(𝑄))^1/𝑠≤𝛽₁≤𝛼₁≤𝛼₂≤𝛽₂≤(𝜆_min(𝐴^∗𝐴))^{1/(𝑠−𝑡)}. Then from (9) we obtain (5).

If (1) has an HPD solution, its eigenvalues may skip between [𝛽₁,𝛼₁] and [𝛼₂,𝛽₂]. Next, what we take more attention on is the HPD solution with its eigenvalues distributed only on one interval.

Theorem 4. Suppose that 𝜆_max(𝐴^∗𝐴)≤(𝑠/𝑡)((𝑡−𝑠)/𝜆_max(𝑄)𝑡)^{(𝑡−𝑠)/𝑠}.(1)Equation (1) has an HPD solution, 𝑋∈[𝛽₁𝐼,𝛼₁𝐼], and if 𝜆_min(𝐴^∗𝐴)>𝑠𝛼₁^𝑠−1(𝑡𝛽₁^𝑡−1)⁻¹ such 𝑋 exists uniquely.(2)Equation (1) has an HPD solution, 𝑍∈[𝛼₂𝐼,𝛽₂𝐼], and if 𝜆_min(𝐴^∗𝐴)>𝑠𝛽₂^𝑠−1(𝑡𝛼₂^𝑡−1)⁻¹ such 𝑍 exists uniquely.

Proof. (1) Let ℎ₁(𝑋)=(𝑄+𝐴^∗𝑋^𝑡𝐴)^1/𝑠, where 𝑋∈[(𝜆_min(𝑄))^1/𝑠𝐼,(𝑠/(𝜆_max(𝐴^∗𝐴)𝑡))^{1/(𝑡−𝑠)}𝐼]. Lemmas 1 and 2 and 𝑡−𝑠>0 imply 𝜆_min(𝑄)^1/𝑠𝐼≤ℎ₁(𝑋)≤𝜆_max(𝑄)+𝜆_max𝐴^∗𝐴𝑠𝜆_max(𝐴^∗𝐴)𝑡^{𝑡/(𝑡−𝑠)1/𝑠}𝐼≤𝑠𝜆_max(𝐴^∗𝐴)𝑡^{𝑠/(𝑡−𝑠)×1/𝑠}𝑠𝐼=𝜆_max(𝐴^∗𝐴)𝑡^{1/(𝑡−𝑠)}𝐼. (10) Applying Brouwer’s fixed-point theorem, ℎ₁(𝑋) has a fixed point 𝑋∈[(𝜆_min(𝑄))^1/𝑠𝐼,(𝑠/(𝜆_max(𝐴^∗𝐴)𝑡))^{1/(𝑡−𝑠)}𝐼]. Then from Theorem 3, 𝑋∈[𝛽₁𝐼,𝛼₁𝐼].
We now prove the uniqueness of 𝑋 under the additional condition that 𝜆_min(𝐴^∗𝐴)>𝑠𝛼₁^𝑠−1(𝑡𝛽₁^𝑡−1)⁻¹. Suppose 𝑌∈[(𝜆_min(𝑄))^1/𝑠𝐼,(𝑠/(𝜆_max(𝐴^∗𝐴)𝑡))^{1/(𝑡−𝑠)}𝐼] is another HPD solution of (1) and 𝑌≠𝑋. It has been known that 𝑋^𝑡−𝑌^𝑡_𝐹=𝐴^−1∗(𝑋^𝑠−𝑌^𝑠)𝐴⁻¹_𝐹≤𝜆_min𝐴^∗𝐴⁻¹𝑋^𝑠−𝑌^𝑠_𝐹. (11) Then from ‖𝑋^𝑠−𝑌^𝑠‖_𝐹≤𝑠𝛼₁^𝑠−1‖𝑋−𝑌‖_𝐹 and ‖𝑋^𝑡−𝑌^𝑡‖_𝐹≥𝑡𝛽₁^𝑡−1‖𝑋−𝑌‖_𝐹, ‖𝑋−𝑌‖_𝐹≤𝑠𝛼₁^𝑠−1𝑡𝛽₁^𝑡−1𝜆_min𝐴^∗𝐴⁻¹‖𝑋−𝑌‖_𝐹<‖𝑋−𝑌‖_𝐹, (12) which is impossible. Hence, 𝑋=𝑌.
(2) Let ℎ₂(𝑍)=[(𝐴⁻¹)^∗(𝑍^𝑠−𝑄)𝐴⁻¹]^1/𝑡, where 𝑍∈[𝛼₂𝐼,𝛽₂𝐼]. ℎ₂(𝑍) is continuous, and ℎ₂𝛼₂𝐼≤ℎ₂(𝑍)≤ℎ₂𝛽₂𝐼 (13) because (𝐴⁻¹)^∗(𝛼^𝑠₂𝐼−𝑄)𝐴⁻¹≤(𝐴⁻¹)^∗(𝑍^𝑠−𝑄)𝐴⁻¹≤(𝐴⁻¹)^∗(𝛽^𝑠₂𝐼−𝑄)𝐴⁻¹. By Lemmas 1 and 2 and Brouwer’s fixed-point theorem, it is sufficient to prove ℎ₂(𝛼₂𝐼)≥𝛼₂𝐼 and ℎ₂(𝛽₂𝐼)≤𝛽₂𝐼 in order for an HPD solution 𝑍∈[𝛼₂𝐼,𝛽₂𝐼] to exist. The existence of such 𝑍 follows from inequalities ℎ₂𝛼₂𝐼=𝐴^−1∗𝛼^𝑠₂𝐼−𝑄𝐴^−11/𝑡≥𝐴^−1∗𝛼^𝑠₂𝐼−𝜆_max(𝑄)𝐼𝐴^−11/𝑡≥𝜆_max𝐴^∗𝐴⁻¹𝛼^𝑠₂𝐼−𝜆_max(𝑄)𝐼^1/𝑡=𝛼₂ℎ𝐼,₂𝛽₂𝐼=𝐴^−1∗𝛽^𝑠₂𝐼−𝑄𝐴^−11/𝑡≤𝐴^−1∗𝛽^𝑠₂𝐼−𝜆_min(𝑄)𝐼𝐴^−11/𝑡≤𝜆_min𝐴^∗𝐴⁻¹𝛽^𝑠₂𝐼−𝜆_min(𝑄)𝐼^1/𝑡=𝛽₂𝐼. (14)
Next we prove the uniqueness of 𝑍 under the additional condition that 𝜆_min(𝐴^∗𝐴)>𝑠𝛽₂^𝑠−1(𝑡𝛼₂^𝑡−1)⁻¹. Suppose (1) has two different HPD solutions 𝑍 and 𝑌 on [𝛼₂𝐼,𝛽₂𝐼]. Then 𝑍^𝑡−𝑌^𝑡_𝐹=𝐴^−1∗(𝑍^𝑠−𝑌^𝑠)𝐴⁻¹_𝐹≤𝜆_min𝐴^∗𝐴⁻¹𝑍^𝑠−𝑌^𝑠_𝐹≤𝜆_min𝐴^∗𝐴⁻¹𝑠𝛽₂^𝑠−1‖𝑍−𝑌‖_𝐹. (15) Moreover, if 𝜆_min(𝐴^∗𝐴)>𝑠𝛽₂^𝑠−1(𝑡𝛼₂^𝑡−1)⁻¹, applying the inequality ‖𝑍^𝑡−𝑌^𝑡‖_𝐹≥𝑡𝛼₂^𝑡−1‖𝑍−𝑌‖_𝐹, we have ‖𝑍−𝑌‖_𝐹≤𝑡𝛼₂^𝑡−1𝜆_min𝐴^∗𝐴⁻¹𝑠𝛽₂^𝑠−1‖𝑍−𝑌‖_𝐹<‖𝑍−𝑌‖_𝐹, (16) which is impossible. Hence, 𝑌=𝑍.

The maximal solution (see, e.g., [30, 31]) of (1) is defined as follows.

Definition 5. An HPD solution 𝑋_𝑀∈ℂ^𝑛×𝑛 of (1) is the maximal solution if, for any HPD solution 𝑌∈ℂ^𝑛×𝑛 of (1), there is 𝑋_𝑀≥𝑌.

So the second term of Theorem 4 implies that the maximal solution of (1) is on [𝛼₂𝐼,𝛽₂𝐼].

Theorem 6. Suppose that 𝜆_max(𝐴^∗𝐴)≤(𝑠/𝑡)(𝑡/(𝑡−𝑠)𝜆_max(𝑄))^{(𝑠−𝑡)/𝑠} and 𝜆_min(𝐴^∗𝐴)>𝑠𝛽₂^𝑠−1(𝑡𝛼₂^𝑡−1)⁻¹; then (1) has a maximal solution 𝑋_max∈[𝛼₂𝐼,𝛽₂𝐼] which can be computed by 𝑋_𝑖=𝐴^−1∗𝑋^𝑠_𝑖−1−𝑄𝐴^−11/𝑡,𝑖=1,2,…. (17) with the initial value 𝑋₀=𝛽₂𝐼.

Proof. Let 𝜉=(𝑡𝛼₂^𝑡−1𝜆_min(𝐴^∗𝐴))⁻¹𝑠𝛽₂^𝑠−1; then 𝜉<1. From the proof of Theorem 4 (2), 𝑡𝛼₂^𝑡−1𝑋_𝑖+1−𝑋_𝑖𝐹≤𝑋^𝑡_𝑖+1−𝑋^𝑡_𝑖𝐹=𝐴^−1∗(𝑋^𝑠_𝑖−𝑋^𝑠_𝑖−1)𝐴⁻¹_𝐹≤𝜆_min𝐴^∗𝐴⁻¹𝑋^𝑠_𝑖−𝑋^𝑠_𝑖−1𝐹≤𝜆_min𝐴^∗𝐴⁻¹𝑠𝛽₂^𝑠−1𝑋_𝑖−𝑋_𝑖−1𝐹. (18) Then 𝑋_𝑖+1−𝑋_𝑖𝐹≤𝜉𝑋_𝑖−𝑋_𝑖−1𝐹≤𝜉^𝑖𝑋₁−𝑋_0𝐹, (19) which indicates the convergence of matrix series {𝑋₀,𝑋₁,𝑋₂,…}, generated by (17).
Set 𝑋₀=𝛽₂𝐼. Assuming 𝑋_𝑖∈[𝛼₂𝐼,𝛽₂𝐼], then from inequalities (14) we have 𝛼₂𝐼≤ℎ𝛼₂𝐼≤𝑋_𝑖+1=𝐴^−1∗𝑋^𝑠_𝑖−𝑄𝐴^−11/𝑡≤ℎ𝛽₂𝐼≤𝛽₂𝐼. (20) That means, for any 𝑖=0,1,2,…,𝑋_𝑖∈[𝛼₂𝐼,𝛽₂𝐼]. By Theorem 4 (2), we can see that 𝑋_max=lim_𝑖→+∞𝑋_𝑖 is the unique HPD solution of (1) on [𝛼₂𝐼,𝛽₂𝐼].
Now we prove the maximality of 𝑋_max. Suppose that 𝑋 is an arbitrary HPD solution of (1); then 𝑋₀≥𝑋, and Theorem 3 implies 𝑋^𝑡₀≥𝑋^𝑡 (since 𝑋₀=𝛽₂𝐼). Assuming that 𝑋^𝑡_𝑖≥𝑋^𝑡, Lemma 1 with 𝑠/𝑡<1 implies 𝑋^𝑡_𝑖+1=𝐴^−1∗𝑋^𝑡_𝑖^𝑠/𝑡−𝑄𝐴⁻¹≥𝐴^−1∗𝑋^{𝑡𝑠/𝑡}−𝑄𝐴⁻¹=𝑋^𝑡. (21) Then 𝑋^𝑡_max=lim_𝑖→+∞𝑋^𝑡_𝑖≥𝑋^𝑡, which implies that 𝑋_max≥𝑋 by the Löwner-Heinz inequality.

Note that similar iteration formula ever appeared in some papers such as [20, 21] for other nonlinear matrix equations. Here we firstly proved that the iteration form (17) preserves the maximality of 𝑋_𝑖 over all HPD solutions of (1).

2.2. Unique Solution of (1) with 𝑠≥𝑡

If 𝑠>𝑡, Lee and Lim [20, Theorem 9.4] show that (1) always has a unique HPD solution, denoted by 𝑋_𝑢. Now we give an upper bound and a lower bound of 𝑋_𝑢 and suggest an iteration method for computing 𝑋_𝑢.

As defined in (3), 𝑔₁(𝑥) and 𝑔₂(𝑥) with 𝑠>𝑡 have unique positive roots, denoted by 𝛾₁ and 𝛾₂, respectively.

Since 𝑔₁(𝜆(𝑋_𝑢))≤0 and 𝑔₂(𝜆(𝑋_𝑢))≤0, 𝛾₂≤𝜆(𝑋_𝑢)≤𝛾₁.

Theorem 7. If 𝑠>𝑡, (1) has a unique HPD solution 𝑋_𝑢∈[𝛾₂𝐼,𝛾₁𝐼]. Let 𝑋₀=𝛾₁𝐼 or 𝛾₂𝐼, then matrix series {𝑋₀,𝑋₁,𝑋₂,…} generated by 𝑋_𝑖=𝑄+𝐴^∗𝑋^𝑡_𝑖−1𝐴^1/𝑠,𝑖=0,1,2,… (22) will converge to 𝑋_𝑢.

Proof. We only need to prove the convergence of matrix series {𝑋₀,𝑋₁,𝑋₂,…}. Set 𝑋₀=𝛾₁𝐼. From (22) we have 𝑋₁=𝑄+𝛾^𝑡₁𝐴^∗𝐴^1/𝑠≤𝜆_max(𝑄)+𝛾^𝑡₁𝜆_max𝐴^∗𝐴^1/𝑠𝐼=𝛾₁𝐼, (23) and then 𝑋^𝑠₁≤𝑋^𝑠₀. Assuming that 𝑋^𝑠_𝑖≤𝑋^𝑠_𝑖−1, 𝑋^𝑠_𝑖+1=𝑄+𝐴^∗𝑋^𝑡_𝑖𝐴=𝑄+𝐴^∗𝑋^𝑠_𝑖^𝑡/𝑠𝐴≤𝑄+𝐴^∗𝑋^𝑠_𝑖−1^𝑡/𝑠𝐴=𝑋^𝑠_𝑖. (24) Then for any 𝑖=0,1,2,…, we have 𝑋^𝑠_𝑖+1≤𝑋^𝑠_𝑖 and then 𝑋_𝑖+1≤𝑋_𝑖 by Löwner-Heinz inequality. On the other hand, 𝑋₀≥𝛾₂𝐼 implies 𝑋_𝑖≥𝛾₂𝐼 for any 𝑖=0,1,2,…, because if 𝑋_𝑖−1≥𝛾₂𝐼, then 𝑋_𝑖=𝑄+𝐴^∗𝑋^𝑡_𝑖−1𝐴^1/𝑠≥𝑄+𝛾^𝑡₂𝐴^∗𝐴^1/𝑠≥𝜆_min(𝑄)+𝛾^𝑡₂𝜆_min𝐴^∗𝐴^1/𝑠𝐼=𝛾₂𝐼. (25) Then {𝑋₀,𝑋₁,𝑋₂,…} with 𝑋₀=𝛾₁𝐼 is a decreasingly monotone matrix series with a lower bound 𝛾₂𝐼. Similarly we can prove that {𝑋₀,𝑋₁,𝑋₂,…} generated by (22) with 𝑋₀=𝛾₂𝐼 is an increasingly monotone matrix series with an upper bound 𝛾₁𝐼. Therefore, the convergence of {𝑋₀,𝑋₁,𝑋₂,…} has been proved.

From the above proof, we can see that the iteration form (22) preserves the minimality (𝑋₀=𝛾₁𝐼) or maximality (𝑋₀=𝛾₂𝐼) of 𝑋_𝑖 in process.

If 𝑠=𝑡, (1) can be reduced to a linear matrix equation 𝑌−𝐴^∗𝑌𝐴=𝑄, which is the discrete-time algebraic Lyapunov equation (DALE) or Hermitian Stein equation, [1, Page 5], assuming that 𝑌=𝑋^𝑠. It is well known that if 𝐴 is d-stable (see [1]), 𝑌−𝐴^∗𝑌𝐴=𝑄 has a unique solution, and matrix series {𝑌₀,𝑌₁,𝑌₂,…}, generated by 𝑌_𝑖+1=𝑄+𝐴^∗𝑌_𝑖𝐴 with an initial value 𝑌₀, will converge to the unique solution. Besides, it is not difficult to get an expression of the unique solution 𝑋_𝑢=(∑^∞_𝑗=0(𝐴^∗)^𝑗𝑄𝐴^𝑗)^1/𝑠, applying [32, Theorem 1, Section 13.2), [1, Theorem 1.1.18], and the results in Section 6.4 in [28].

Now we have presented the solvability theory of the self-adjoint polynomial matrix equation (1) in three cases. A general iterative algorithm for its maximal solution (𝑠<𝑡) or unique solution (𝑠≥𝑡) will be given in Section 4. Before it, we study the algebraic perturbation of the maximal or unique solution of (1).

3. Algebraic Perturbation Analysis

In this section, we present the algebraic perturbation analysis of the HPD solution of (1) with respect to the perturbation of its coefficient matrices. Similar to [30], we define the perturbed matrix equation of (1) as 𝑋^𝑠−𝐴^∗𝑋^𝑡𝐴=𝑄, (26) where 𝐴=𝐴+Δ𝐴∈ℂ^𝑛×𝑛 and 𝑄=𝑄+Δ𝑄∈ℂ^𝑛×𝑛. We always suppose that (1) has a maximal (or unique) solution, denoted by 𝑋_𝑀∈[𝛼₂𝐼,𝛽₂𝐼], and (26) has a maximal (or unique) solution, denoted by 𝑋_𝑀∈[𝛼₂𝐼,𝛽₂𝐼].

Now we present the perturbation bound for 𝑋_𝑀 when 𝑠≠𝑡. Define a function 𝜏: 𝜏𝛼,𝛽=𝑠𝛼^𝑠−1−𝑡𝛽^𝑡−1‖𝐴‖²₂,𝛼,𝛽∈ℝ². (27)

Theorem 8. Let 𝜀>0 be an arbitrary real number, and 𝜏(𝛼₂,𝛽₂)≥0. If ‖Δ𝐴‖_𝐹<‖𝐴‖²₂+2𝜀3𝜏𝛼₂,𝛽₂𝑋_𝑀2^−𝑡1/2−‖𝐴‖₂,‖Δ𝑄‖_𝐹<13𝜏𝛼₂,𝛽₂𝜀, (28) then 𝑋_𝑀−𝑋_𝑀𝐹<𝜀. (29)

Proof. It is easy to induce that 𝑋^𝑠_𝑀−𝑋^𝑠_𝑀𝐹≥^𝑠−1_𝑘=0𝛼₂^{𝑠−1−𝑘}𝛼^𝑘₂𝑋_𝑀−𝑋_𝑀𝐹≥𝑠𝛼₂^𝑠−1𝑋_𝑀−𝑋_𝑀𝐹,𝑋^𝑡_𝑀−𝑋^𝑡_𝑀𝐹≤^𝑡−1_𝑘=0𝛽₂^{𝑡−1−𝑘}𝛽^𝑘₂𝑋_𝑀−𝑋_𝑀𝐹≤𝑡𝛽₂^𝑡−1𝑋_𝑀−𝑋_𝑀𝐹. (30) Then from (1) and (26), we have 𝜏𝛼₂,𝛽₂𝑋_𝑀−𝑋_𝑀𝐹≤2‖𝐴‖₂𝑋_𝑀^𝑡₂‖Δ𝐴‖_𝐹+𝑋_𝑀^𝑡₂‖Δ𝐴‖²_𝐹+‖Δ𝑄‖_𝐹. (31) Since 𝜏(𝛼₂,𝛽₂)>0, 𝑋_𝑀−𝑋_𝑀𝐹≤𝜏𝛼₂,𝛽₂⁻¹2‖𝐴‖₂𝑋_𝑀^𝑡₂‖Δ𝐴‖_𝐹+𝑋_𝑀^𝑡₂‖Δ𝐴‖²_𝐹+‖Δ𝑄‖_𝐹. (32) Then for an arbitrary 𝜀>0, if ‖Δ𝐴‖_𝐹<(‖𝐴‖²₂+(2𝜀/3)𝜏(𝛼₂,𝛽₂)‖𝑋_𝑀‖₂^−𝑡)^1/2−‖𝐴‖₂ and ‖Δ𝑄‖_𝐹<(1/3)𝜏(𝛼₂,𝛽₂)𝜀, we have (29).

If 𝑠=𝑡, for an arbitrary 𝜀>0, define 𝜚(𝜀)=‖𝐴‖₂+‖𝐴‖²₂+2𝜀3𝜌^1/2, (33) where 𝜌=𝑋_𝑀^𝑠₂𝑠𝛼₂^𝑠−11−‖𝐴‖²₂⁻¹. (34)

Theorem 9. Let 𝜀>0 be an arbitrary real number, and ‖𝐴‖₂<1. If ‖Δ𝐴‖_𝐹<2𝜀3𝜌𝜚(𝜀)⁻¹,‖Δ𝑄‖_𝐹<𝜀𝑋3𝜌^𝑀𝑠2, (35) then 𝑋_𝑀−𝑋_𝑀𝐹<𝜀. (36)

Proof. Similar to the proof of Theorem 8, we can induce that 1−‖𝐴‖²₂𝑋^𝑠_𝑀−𝑋^𝑠_𝑀𝐹≤2‖𝐴‖₂𝑋_𝑀^𝑠₂‖Δ𝐴‖_𝐹+𝑋_𝑀^𝑠₂‖Δ𝐴‖²_𝐹+‖Δ𝑄‖_𝐹. (37) Then 𝑋^𝑠_𝑀−𝑋^𝑠_𝑀𝐹≤1−‖𝐴‖²₂⁻¹2‖𝐴‖₂𝑋_𝑀^𝑠₂‖Δ𝐴‖_𝐹+𝑋_𝑀^𝑠₂‖Δ𝐴‖²_𝐹+‖Δ𝑄‖_𝐹. (38) With the help of (30) and (34), (38) implies 𝑋_𝑀−𝑋_𝑀𝐹≤𝜌‖Δ𝐴‖_𝐹+2‖𝐴‖₂‖Δ𝐴‖_𝐹+‖Δ𝑄‖_𝐹. (39) Then if ‖Δ𝐴‖_𝐹<(2𝜀/3)(𝜌𝜚(𝜀))⁻¹ and ‖Δ𝑄‖_𝐹<(𝜀/3𝜌)‖𝑋_𝑀‖^𝑠₂, we have (36).

Theorems 8 and 9 make sure that the perturbation of 𝑋_𝑀 can be controlled if Δ𝐴 and Δ𝑄 have a proper upper bound.

4. Algorithm and Numerical Experiments

In this section we give a general iterative algorithm for the maximal or unique solutions of (1) and two numerical experiments. All reported results were obtained using MATLAB-R2012b on a personal computer with 2.4 GHz Intel Core i7 and 8 GB 1600 MHz DDR3.

Example 10. Let matrices 𝐴=rand(100)×10⁻² and 𝑄=eye(100). With tol=10⁻¹² and not more than 200 iterations, we apply Algorithm 1 to compute the maximal or unique HPD solutions of (1) with 𝑠≠𝑡 and compare the results with those by the iteration method from [33] (denoted by MONO in Table 1).
Table 1 shows iterations, CPU times before convergence, and the residues of the computed HPD solution 𝑋, defined by 𝑒()=𝑋𝑠,𝑡^𝑠−𝐴^∗𝑋^𝑡𝐴−𝑄_𝐹‖[𝐴,𝑄]‖_𝐹. (40)

Table 1 

Iteration, CPU time (seconds) and residue for solving (1) with 𝑠≠𝑡.

Step  1. Compute  𝜆max(𝐴∗𝐴), 𝜆min(𝐴∗𝐴), 𝜆max(𝑄), 𝜆min(𝑄). Step  2. Input (3). Step  3. If  𝑠<𝑡, run Steps  4-5; if  𝑡<𝑠, run Steps  6-7; otherwise, run Steps  8-9. Step  4. Compute the roots  𝛼1, 𝛼2  of  𝑔1(𝑥), and  𝛽1, 𝛽2  of  𝑔2(𝑥), respectively. Step  5. Let  𝑋0=𝛽2𝐼, run (17). Step  6. Compute the root  𝛾1  of  𝑔1(𝑥)  and the root  𝛾2  of  𝑔2(𝑥), respectively. Step  7. Let  𝑍0=𝛾1𝐼, run (22). Step  8. Compute the root  𝛿1  of  𝑔1(𝑥)  and the root  𝛿2  of  𝑔2(𝑥), respectively. Step  9. If  𝜆max(𝐴∗𝐴)<1  and  𝛿1≥𝛿2, then let  𝑋0=𝛿1𝐼  and run (22).

Algorithm 1 

Given matrices 𝐴,𝑄∈ℂ^𝑛×𝑛 and positive integers 𝑠,𝑡.

From Table 1, we can see that it takes more iterations and CPU times to solve the maximal solution of (1) with 𝑠<𝑡 than to solve the unique solution of (1) with 𝑠>𝑡. At the same time, the accuracy of the latter is better than the former. MONO can not be used to solve (1) with 𝑠<𝑡, and it costs more iterations and CPU times than Algorithm 1 when solving (1) with 𝑠>𝑡.

Now we use Example 4.1 of [33] to test our method.

Example 11. Let 𝐴=0.5𝐵/‖𝐵‖_∞ with 𝐵=[𝐵_𝑖𝑗]_𝑛×𝑛,𝑏_𝑖𝑗=𝑖+𝑗+1 and let 𝑄=eye(𝑛), with 𝑛=100. We solve (1) with 𝑠=𝑡 and with two different initial solutions. The iterations, CPU times, and the residues of the computation are reported in Table 2.

Table 2 

Iteration, CPU time (seconds) and residue for solving (1) with 𝑠=𝑡 and different initial solutions.

Table 2 shows that for Algorithm 1 the choice 𝑋₀=𝛿₁𝐼_𝑛 is better than 𝑋₀=𝛿₂𝐼_𝑛. When 𝑠 and 𝑡 rise, MONO might lose its efficiency. It seems not proper to apply the iteration method designed for 𝑌−𝐴^∗𝑌^𝑡/𝑠𝐴=𝑄 with 𝑌=𝑋^𝑠 to solve 𝑋^𝑠−𝐴^∗𝑋^𝑡𝐴=𝑄, although they are equivalent to each other in theory.

5. Conclusion

In this paper, we considered the solvability of the self-adjoint polynomial matrix equation (1). Sufficient conditions were given to guarantee the existence of the maximal or unique HPD solutions of (1). The algebraic perturbation analysis including perturbation bounds was also developed for (1) under the perturbation of given coefficient matrices. At last a general iterative algorithm with maximality preserved in process was presented for the maximal or unique solution with two numerical experiments reported.

Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.

Acknowledgments

Zhigang Jia’s research was supported in part by National Natural Science Foundation of China under Grants 11201193 and 11171289 and a project funded by the Priority Academic Program Development of Jiangsu Higher Education Institutions. Minghui Wang’s research was supported in part by the National Natural Science Foundation of China (Grant no. 11001144), the Science and Technology Program of Shandong Universities of China (J11LA04), and the Research Award Fund for Outstanding Young Scientists of Shandong Province in China (BS2012DX009). Sitao Ling’s research was supported in part by National Natural Science Foundations of China under Grant 11301529, Postdoctoral Science Foundation of China under Grant 2013M540472, and Jiangsu Planned Projects for Postdoctoral Research Funds 1302036C. The authors would like to thank three anonymous referees for giving valuable comments and suggestions.