Perturbation Analysis of the Nonlinear Matrix Equation 𝑋−∑^𝑚_𝑖=1𝐴^∗_𝑖𝑋^𝑝𝑖𝐴_𝑖=𝑄

Abstract

Consider the nonlinear matrix equation 𝑋−∑^𝑚_𝑖=1𝐴^∗_𝑖𝑋^𝑝𝑖𝐴_𝑖=𝑄 with 0<𝑝_𝑖<1. Two perturbation bounds and the backward error of an approximate solution to the equation are derived. Explicit expressions of the condition number for the equation are obtained. The theoretical results are illustrated by numerical examples.

1. Introduction

In this paper we consider the Hermitian positive definite solution of the nonlinear matrix equation 𝑋−^𝑚_𝑖=1𝐴^∗_𝑖𝑋^𝑝𝑖𝐴_𝑖=𝑄, (1) where 0<𝑝_𝑖<1(𝑖=1,2,…,𝑚), 𝐴₁,𝐴₂,…,𝐴_𝑚 are 𝑛×𝑛 complex matrices, 𝑚 is a positive integer, and 𝑄 is a positive definite matrix. Here, 𝐴^∗_𝑖 denotes the conjugate transpose of the matrix 𝐴_𝑖.

When 𝑚>1, (1) is recognized as playing an important role in solving a system of linear equations. For example, in many physical calculations, one must solve the system of linear equation 𝑀𝑥=𝑓, (2) where 𝑀=𝐼0⋯0𝐴₁0𝐼⋯0𝐴₂⋮⋮⋱⋮⋮00⋯𝐼𝐴_𝑚𝐴^∗₁𝐴^∗₂⋯𝐴^∗_𝑚−𝑄 (3) arises in a finite difference approximation to an elliptic partial differential equation (for more information, refer to [1]). We can rewrite 𝑀 as 𝑀=𝑀+𝐷, where 𝑋𝑀=^−𝑝10⋯0𝐴₁0𝑋^−𝑝2⋯0𝐴₂⋮⋮⋱⋮⋮00⋯𝑋^−𝑝𝑚𝐴_𝑚𝐴^∗₁𝐴^∗₂⋯𝐴^∗_𝑚,−𝑄𝐷=𝐼−𝑋^−𝑝10⋯000𝐼−𝑋^−𝑝2⋯00⋮⋮⋱⋮⋮00⋯𝐼−𝑋^−𝑝𝑚0.00⋯00 (4)𝑀 can be factored as 𝑀=−𝐼0⋯000−𝐼⋯00⋮⋮⋱⋮⋮00⋯−𝐼0−𝐴^∗₁𝑋^𝑝1−𝐴^∗₂𝑋^𝑝2⋯−𝐴^∗_𝑚𝑋^𝑝𝑚×−𝐼−𝑋^−𝑝10⋯0−𝐴₁0−𝑋^−𝑝2⋯0−𝐴₂⋮⋮⋱⋮⋮00⋯−𝑋^−𝑝𝑚−𝐴_𝑚00⋯0𝑋 (5) if and only if 𝑋 is a solution of equation 𝑋−∑^𝑚_𝑖=1𝐴^∗_𝑖𝑋^𝑝𝑖𝐴_𝑖=𝑄. When   𝑚=1, this type of nonlinear matrix equations arises in ladder networks, dynamic programming, control theory, stochastic filtering, statistics, and so forth [2–7].

For the similar equations 𝑋±𝐴^∗𝑋^−𝑝𝐴=𝑄, 𝑋^𝑠±𝐴^∗𝑋^−𝑡𝐴=𝑄, 𝑋+∑^𝑚_𝑖=1𝐴^∗_𝑖𝑋⁻¹𝐴_𝑖=𝐼, and 𝑋=𝑄−𝐴^∗𝑋⁻¹𝐴+𝐵^∗𝑋⁻¹𝐵, there were many contributions in the literature to the theory, numerical solutions, and perturbation analysis [8–32]. Jia and Gao [33] derived two perturbation estimates for the solution of the equation 𝑋−𝐴^∗𝑋^𝑞𝐴=𝑄 with 0<𝑞<1. In addition, Duan et al. [34] proved that the equation 𝑋−∑^𝑚_𝑖=1𝐴^∗_𝑖𝑋^𝛿𝑖𝐴_𝑖=𝑄(0<|𝛿_𝑖|<1) has a unique positive definite solution. They also proposed an iterative method for obtaining the unique positive definite solution. However, to our best knowledge, there has been no perturbation analysis for (1) with 𝑚>1 in the known literature.

As a continuation of the previous results, the rest of the paper is organized as follows. In Section 2, some preliminary lemmas are given. In Section 3, two perturbation bounds for the unique solution to (1) are derived. Furthermore, in Section 4, we obtain the backward error of an approximate solution to (1). In Section 5, we also discuss the condition number of the unique solution to (1). Finally, several numerical examples are presented in Section 6.

We denote by 𝒞^𝑛×𝑛 the set of 𝑛×𝑛 complex matrices, by ℋ^𝑛×𝑛 the set of 𝑛×𝑛 Hermitian matrices, by 𝐼 the identity matrix, by 𝐢 the imaginary unit, by ‖⋅‖ the spectral norm, by ‖⋅‖_𝐹 the Frobenius norm, and by 𝜆_max(𝑀) and 𝜆_min(𝑀) the maximal and minimal eigenvalues of 𝑀, respectively. For 𝐴=(𝑎₁,…,𝑎_𝑛)=(𝑎_𝑖𝑗)∈𝒞^𝑛×𝑛 and a matrix 𝐵, 𝐴⊗𝐵=(𝑎_𝑖𝑗𝐵) is a Kronecker product, and vec𝐴 is a vector defined by vec𝐴=(𝑎^𝑇₁,…,𝑎^𝑇_𝑛)^𝑇. For 𝑋,𝑌∈ℋ^𝑛×𝑛, we write 𝑋≥𝑌 (𝑋>𝑌,resp.) if 𝑋−𝑌 is Hermitian positive semidefinite (definite, resp.).

2. Preliminaries

Lemma 1 (see [35]). If 𝐴≥𝐵>0 and 0≤𝛾≤1, then 𝐴^𝛾≥𝐵^𝛾.

Lemma 2 (see [33]). For any Hermitian positive definite matrix 𝑋 and Hermitian matrix Δ𝑋, one has (i)𝑋^𝑞=(sin𝑞𝜋/𝜋)∫^∞₀𝑋^1/2(𝜆𝐼+𝑋)⁻¹𝑋^1/2𝜆^𝑞−1𝑑𝜆,  0<𝑞<1;(ii)𝑋^𝑞=(sin𝑞𝜋/(1−𝑞)𝜋)∫^∞₀𝑋^1/2(𝜆𝐼+𝑋)⁻¹𝑋(𝜆𝐼+𝑋)⁻¹𝑋^1/2𝜆^𝑞−1𝑑𝜆,0<𝑞<1. In addition, if 𝑋+Δ𝑋≥(1/𝜈)𝑋>0 and 0<𝑞<1, then (iii)‖𝑋^−1/2𝐴^∗((𝑋+Δ𝑋)^𝑞−𝑋^𝑞)𝐴𝑋^−1/2‖≤𝑞(‖𝑋^−1/2Δ𝑋𝑋^−1/2‖)+𝜈‖𝑋^−1/2Δ𝑋𝑋^−1/2‖²)‖𝑋^𝑞/2𝐴𝑋^−1/2‖².

Lemma 3 (see [34]). The matrix equation 𝑋−∑^𝑚_𝑖=1𝐴^∗_𝑖𝑋^𝛿𝑖𝐴_𝑖=𝑄(0<|𝛿_𝑖|<1) always has a unique positive definite solution 𝑋. The matrix sequence 𝑋_𝑘: 𝑋_𝑠+𝑚+1=𝑄+^𝑚_𝑖=1𝐴^∗_𝑖𝑋^𝛿𝑖_𝑠+𝑖𝐴_𝑖,𝑠=0,1,2,…, (6) converges to the unique positive definite solution 𝑋 for arbitrary initial positive definite matrices 𝑋₁,𝑋₂,…,𝑋_𝑚.

3. Perturbation Bounds

Here the perturbed equation 𝑋−^𝑚_𝑖=1𝐴^∗_𝑖𝑋^𝑝𝑖𝐴_𝑖=𝑄. (7) is considered, where 0<𝑝_𝑖<1  and  𝐴_𝑖 and 𝑄 are small perturbations of 𝐴_𝑖 and 𝑄 in (1), respectively. We assume that 𝑋 and 𝑋 are the solutions of (1) and (7), respectively. Let Δ𝑋=𝑋−𝑋, Δ𝑄=𝑄−𝑄 and Δ𝐴_𝑖=𝐴_𝑖−𝐴_𝑖.

By Lemma 3, we know that (1) always has a unique positive definite solution 𝑋; then in this section two perturbation bounds for the unique positive definite solution of (1) are developed. The relative perturbation bound in Theorem 5 does not depend on any knowledge of the actual solution 𝑋 of (1). Furthermore, a sharper perturbation bound in Theorem 8 is derived.

To prove the next theorem, we first verify the following lemma.

Lemma 4. If 𝑋 is a solution of (1), then 𝑋≥𝜆_min(𝑄)+^𝑚_𝑖=1𝜆_min𝐴^∗_𝑖𝐴_𝑖𝜆^𝑝𝑖_min(𝑄)𝐼=𝛽𝐼. (8)

Proof. By Lemma 3, (1) with 0<𝑝_𝑖<1 always has a unique positive definite solution 𝑋. Then 𝑋>0, and it follows that 𝑋^𝑝𝑖>0. Therefore 𝑋≥𝑄. By Lemma 1 and (1), we have 𝑋≥𝑄+∑^𝑚_𝑖=1𝐴^∗_𝑖𝑄^𝑝𝑖𝐴_𝑖≥  (𝜆_min(𝑄)+∑^𝑚_𝑖=1𝜆_min(𝐴^∗_𝑖𝐴_𝑖)𝜆^𝑝𝑖_min(𝑄))𝐼=𝛽𝐼.

The next theorem generalizes [33, Theorem 4] with 𝑚=1,  ‖Δ𝑄‖=0 to arbitrary integer 𝑚≥1,  ‖Δ𝑄‖>0.

Theorem 5. Let 𝑏=𝛽+‖Δ𝑄‖−∑^𝑚_𝑖=1𝑝_𝑖𝛽^𝑝𝑖‖𝐴_𝑖‖²,𝑠=∑^𝑚_𝑖=1𝛽^𝑝𝑖‖Δ𝐴_𝑖‖(2‖𝐴_𝑖‖+‖Δ𝐴_𝑖‖). If ,𝑏0<𝑏<2𝛽−𝑠₂−4𝛽−𝑠(𝑠+‖Δ𝑄‖)≥0, (9) then 𝑋−𝑋‖𝑋‖≤𝜚^𝑚_𝑖=1Δ𝐴_𝑖+𝜔‖Δ𝑄‖≜𝜉₁, (10) where 𝜚=2𝑠∑^𝑚_𝑖=1Δ𝐴_𝑖𝑏+𝑏²(‖,2−4𝛽−𝑠𝑠+Δ𝑄‖)𝜔=𝑏+𝑏²(‖.−4𝛽−𝑠𝑠+Δ𝑄‖) (11)

Proof. Let Ω=Δ𝑋∈ℋ^𝑛×𝑛∶𝑋^−1/2Δ𝑋𝑋^−1/2≤𝜚^𝑚_𝑖=1Δ𝐴_𝑖+𝜔‖Δ𝑄‖. (12) Obviously, Ω is a nonempty bounded convex closed set. Let 𝑓(Δ𝑋)=^𝑚_𝑖=1𝐴_𝑖^∗(𝑋+Δ𝑋)^𝑝𝑖𝐴_𝑖−𝐴^∗_𝑖𝑋^𝑝𝑖𝐴_𝑖+Δ𝑄,Δ𝑋∈Ω. (13) Evidently, 𝑓∶Ω↦ℋ^𝑛×𝑛 is continuous. We will prove that 𝑓(Ω)⊆Ω.
For every Δ𝑋∈Ω, it follows that ‖𝑋^−1/2Δ𝑋𝑋^−1/2‖≤𝜚∑^𝑚_𝑖=1‖Δ𝐴_𝑖‖+𝜔‖Δ𝑄‖. Thus 𝜚^𝑚_𝑖=1Δ𝐴_𝑖+𝜔‖Δ𝑄‖𝐼≥𝑋^−1/2Δ𝑋𝑋^−1/2≥−𝜚^𝑚_𝑖=1Δ𝐴_𝑖−𝜔‖Δ𝑄‖𝐼,1+𝜚^𝑚_𝑖=1Δ𝐴_𝑖‖‖𝑋≥+𝜔Δ𝑄≥𝑋+Δ𝑋1−𝜚^𝑚_𝑖=1Δ𝐴_𝑖−𝜔‖Δ𝑄‖𝑋. (14) According to (9), we have 𝜚^𝑚_𝑖=1Δ𝐴_𝑖‖‖=2‖)+𝜔Δ𝑄(‖Δ𝑄+𝑠𝑏+𝑏²‖≤2‖)−4𝛽−𝑠(𝑠+Δ𝑄‖)(‖Δ𝑄+𝑠𝑏≤𝑏2𝛽−𝑠<1. (15) Therefore 1−𝜚^𝑚_𝑖=1Δ𝐴_𝑖−𝜔‖Δ𝑄‖𝑋>0. (16) From Lemmas 2 and 4, it follows that 𝑋^−1/2𝑚_𝑖=1𝐴^∗_𝑖(𝑋+Δ𝑋)^𝑝𝑖−𝑋^𝑝𝑖𝐴_𝑖𝑋^−1/2≤𝑋^−1/2Δ𝑋𝑋^−1/2+𝑋^−1/2Δ𝑋𝑋^−1/221−𝜚∑^𝑚_𝑖=1Δ𝐴_𝑖‖‖×−𝜔Δ𝑄_𝑚𝑖=1𝑝_𝑖𝑋^𝑝𝑖/2𝐴_𝑖𝑋^−1/22≤𝑋^−1/2Δ𝑋𝑋^−1/2+𝑋^−1/2Δ𝑋𝑋^−1/221−𝜚∑^𝑚_𝑖=1Δ𝐴_𝑖‖‖×−𝜔Δ𝑄_𝑚𝑖=1𝑝_𝑖𝛽^1−𝑝𝑖𝐴_𝑖². (17) Therefore 𝑋^−1/2𝑓(Δ𝑋)𝑋^−1/2=𝑋^−1/2𝑚_𝑖=1𝐴^∗_𝑖(𝑋+Δ𝑋)^𝑝𝑖𝐴_𝑖−𝐴^∗_𝑖𝑋^𝑝𝑖𝐴_𝑖×𝑋^−1/2+𝑋^−1/2Δ𝑄𝑋^−1/2≤^𝑚_𝑖=1𝑋^−1/2𝐴^∗_𝑖(𝑋+Δ𝑋)^𝑝𝑖−𝑋^𝑝𝑖𝐴_𝑖𝑋^−1/2+𝑋^−1/2Δ𝑄𝑋^−1/2+^𝑚_𝑖=1𝑋^−1/2Δ𝐴^∗_𝑖(𝑋+Δ𝑋)^𝑝𝑖𝐴_𝑖+Δ𝐴_𝑖+𝐴^∗_𝑖(𝑋+Δ𝑋)^𝑝𝑖Δ𝐴_𝑖𝑋^−1/2≤𝑋^−1/2Δ𝑋𝑋^−1/2+𝑋^−1/2Δ𝑋𝑋^−1/221−𝜚∑^𝑚_𝑖=1Δ𝐴_𝑖‖‖×−𝜔Δ𝑄_𝑚𝑖=1𝑝_𝑖𝛽^1−𝑝𝑖𝐴_𝑖²+^𝑚_𝑖=1Δ𝐴_𝑖2𝐴_𝑖+Δ𝐴_𝑖𝛽^1−𝑝𝑖×1+𝜚^𝑚_𝑖=1Δ𝐴_𝑖‖‖+‖‖+𝜔Δ𝑄Δ𝑄𝛽≤𝜉₁+𝜉²₁1−𝜉₁^𝑚_𝑖=1𝑝_𝑖𝛽^1−𝑝𝑖𝐴_𝑖²+𝑠𝛽1+𝜉₁+‖Δ𝑄‖𝛽=𝜉₁. (18) That is, 𝑓(Ω)⊆Ω. By Brouwer's fixed point theorem, there exists a Δ𝑋∈Ω such that 𝑓(Δ𝑋)=Δ𝑋. Moreover, by Lemma 3, we know that 𝑋 and 𝑋 are the unique solutions to (1) and (7), respectively. Then 𝑋−𝑋‖𝑋‖=‖Δ𝑋‖‖𝑋‖=𝑋^1/2𝑋^−1/2Δ𝑋𝑋^−1/2𝑋^1/2‖𝑋‖≤𝑋^−1/2Δ𝑋𝑋^−1/2≤𝜚^𝑚_𝑖=1Δ𝐴_𝑖+𝜔‖Δ𝑄‖. (19)

Remark 6. According to 𝜚^𝑚_𝑖=1Δ𝐴_𝑖‖‖=2∑+𝜔Δ𝑄_𝑚𝑖=1Δ𝐴_𝑖2𝐴_𝑖+Δ𝐴_𝑖+‖Δ𝑄‖𝑏+𝑏²(‖,−4𝛽−𝑠𝑠+Δ𝑄‖) (20) we get 𝜚∑^𝑚_𝑖=1‖Δ𝐴_𝑖‖+𝜔‖Δ𝑄‖→0 for ‖Δ𝑄‖→0 and ‖Δ𝐴_𝑖‖→0(𝑖=1,2,…,𝑚). Therefore (1) is well posed.

Next, a sharper perturbation estimate is derived.

Subtracting (1) from (7), we have Δ𝑋+^𝑚_𝑖=1sin𝑝_𝑖𝜋𝜋^∞₀(𝜆𝐼+𝑋)⁻¹𝑋^1/2𝐴_𝑖^∗×Δ𝑋(𝜆𝐼+𝑋)⁻¹𝑋^1/2𝐴_𝑖𝜆^𝑝𝑖−1(),𝑑𝜆=𝐸+ℎΔ𝑋 (21) where 𝐵_𝑖=𝑋^𝑝𝑖𝐴_𝑖,𝑖=1,2,…,𝑚,𝐸=^𝑚_𝑖=1𝐵^∗_𝑖Δ𝐴_𝑖+Δ𝐴^∗_𝑖𝐵_𝑖+^𝑚_𝑖=1Δ𝐴^∗_𝑖𝑋^𝑝𝑖Δ𝐴_𝑖𝑍+Δ𝑄,_𝑖(Δ𝑋)=sin𝑝_𝑖𝜋𝜋^∞₀𝑋^1/2(𝜆𝐼+𝑋)⁻¹Δ𝑋(𝜆𝐼+𝑋+Δ𝑋)⁻¹×Δ𝑋(𝜆𝐼+𝑋)⁻¹𝑋^1/2𝜆^𝑝𝑖−1𝑉𝑑𝜆,_𝑖(Δ𝑋)=sin𝑝_𝑖𝜋𝜋^∞₀𝑋^1/2(𝜆𝐼+𝑋)⁻¹×()Δ𝑋𝜆𝐼+𝑋+Δ𝑋⁻¹𝑋^1/2𝜆^𝑝𝑖−1ℎ()=𝑑𝜆,Δ𝑋^𝑚_𝑖=1𝐴^∗_𝑖𝑍_𝑖(Δ𝑋)𝐴_𝑖−𝐴^∗_𝑖𝑉_𝑖(Δ𝑋)Δ𝐴_𝑖−Δ𝐴^∗_𝑖𝑉_𝑖(Δ𝑋)𝐴_𝑖. (22)

Lemma 7. If ∑^𝑚_𝑖=1(‖𝐴_𝑖‖²/𝛽^1−𝑝𝑖)<1, then the linear operator 𝐋∶ℋ^𝑛×𝑛→ℋ^𝑛×𝑛 defined by 𝐋𝑊=𝑊+^𝑚_𝑖=1sin𝑝_𝑖𝜋𝜋×^∞₀(𝜆𝐼+𝑋)⁻¹𝑋^1/2𝐴_𝑖^∗×𝑊(𝜆𝐼+𝑋)⁻¹𝑋^1/2𝐴_𝑖×𝜆^𝑝𝑖−1𝑑𝜆,𝑊∈ℋ^𝑛×𝑛 (23) is invertible.

Proof. It suffices to show that the following equation: 𝐋𝑊=𝑉 (24) has a unique solution for every  𝑉∈ℋ^𝑛×𝑛. Define the operator 𝐌∶ℋ^𝑛×𝑛→ℋ^𝑛×𝑛 by 𝐌𝑍=^𝑚_𝑖=1sin𝑝_𝑖𝜋𝜋×^∞₀𝑋^−1/2𝐴^∗_𝑖𝑋^1/2×(𝜆𝐼+𝑋)⁻¹𝑋^1/2𝑍𝑋^1/2×(𝜆𝐼+𝑋)⁻¹𝑋^1/2𝐴_𝑖𝑋^−1/2×𝜆^𝑝𝑖−1𝑑𝜆,𝑍∈ℋ^𝑛×𝑛. (25) Let 𝑌=𝑋^−1/2𝑊𝑋^−1/2. Thus (21) is equivalent to 𝑌+𝐌𝑌=𝑋^−1/2𝑉𝑋^−1/2. (26) According to Lemma 2, we have ‖𝐌𝑌‖≤^𝑚_𝑖=1sin𝑝_𝑖𝜋𝜋×^∞₀𝑋^−1/2𝐴^∗_𝑖𝑋^1/2×(𝜆𝐼+𝑋)⁻¹𝑋(𝜆𝐼+𝑋)⁻¹×𝑋^1/2𝐴_𝑖𝑋^−1/2𝜆^𝑝𝑖−1𝑌‖≤𝑑𝜆‖_𝑚𝑖=11−𝑝_𝑖𝑋^𝑝𝑖/2𝐴_𝑖𝑋^−1/22‖𝑌‖≤^𝑚_𝑖=1𝐴_𝑖²𝛽^1−𝑝𝑖‖𝑌‖<‖𝑌‖, (27) which implies that ||𝐌||<1 and 𝐼+𝐌 is invertible. Therefore, the operator 𝐋 is invertible.

Furthermore, we define operators 𝐏_𝑖∶𝒞^𝑛×𝑛→ℋ^𝑛×𝑛 by 𝐏_𝑖𝑍_𝑖=𝐋⁻¹𝐵^∗_𝑖𝑍_𝑖+𝑍^∗_𝑖𝐵_𝑖,𝑍_𝑖∈𝒞^𝑛×𝑛,𝑖=1,2,…,𝑚. (28) Thus, we can rewrite (21) as Δ𝑋=𝐋⁻¹+Δ𝑄_𝑚𝑖=1𝐏_𝑖Δ𝐴_𝑖+𝐋^−1𝑚_𝑖=1Δ𝐴^∗_𝑖𝑋^𝑝𝑖Δ𝐴_𝑖+𝐋⁻¹(ℎ(Δ𝑋)). (29) Define 𝐋⁻¹=max_{𝑊∈ℋ𝑛×𝑛‖𝑊‖=1}𝐋⁻¹𝑊,𝐏_𝑖=max_{𝑍∈𝒞𝑛×𝑛‖𝑍‖=1}𝐏_𝑖𝑍,𝑖=1,2,…,𝑚. (30) Now we denote that 𝑙=𝐋⁻¹⁻¹,𝜁=𝑋⁻¹,𝜉_𝑖=𝑋^𝑝𝑖,𝑛_𝑖=𝐏_𝑖,𝜁𝜃=²𝑙^𝑚_𝑖=1𝜉_𝑖𝐴_𝑖²,1𝑖=1,2,…,𝑚,𝜖=𝑙‖Δ𝑄‖+^𝑚_𝑖=1𝑛_𝑖Δ𝐴_𝑖+𝜉_𝑖𝑙Δ𝐴_𝑖²,𝜁𝜎=𝑙^𝑚_𝑖=1𝜉_𝑖2𝐴_𝑖+Δ𝐴_𝑖Δ𝐴_𝑖. (31)

Theorem 8. Suppose that 𝑋 and 𝑋 are the solutions of (1) and (7), respectively. If ∑^𝑚_𝑖=1‖𝐴_𝑖‖²/𝛽^1−𝑝𝑖<1, ()𝜎<1,𝜖<1−𝜎²𝜁+𝜎𝜁+2𝜃+2,((𝜁+𝜃)(𝜎𝜁+𝜃)) (32) then 𝑋−𝑋≤2𝜖1+𝜖𝜁−𝜎+(1+𝜁𝜖−𝜎)²−4𝜖(𝜁+𝜃)≡𝜈. (33)

Proof. Let 𝑓(Δ𝑋)=𝐋⁻¹Δ𝑄+^𝑚_𝑖=1𝐏_𝑖Δ𝐴_𝑖+𝐋^−1𝑚_𝑖=1Δ𝐴^∗_𝑖𝑋^𝑝𝑖Δ𝐴_𝑖+𝐋⁻¹(ℎ(Δ𝑋)). (34) Obviously, 𝑓∶ℋ^𝑛×𝑛→ℋ^𝑛×𝑛 is continuous. The condition (32) ensures that the quadratic equation (𝜁+𝜃)𝑥²−(1+𝜁𝜖−𝜎)𝑥+𝜖=0 with respect to the variable 𝑥 has two positive real roots. The smaller one is 𝜈=2𝜖1+𝜖𝜁−𝜎+(1+𝜁𝜖−𝜎)²().−4𝜖𝜁+𝜃 (35) Define Ω={Δ𝑋∈ℋ^𝑛×𝑛∶‖Δ𝑋‖≤𝜈}. Then for any Δ𝑋∈Ω, by (32), we have 𝑋⁻¹Δ𝑋≤𝑋⁻¹‖Δ𝑋‖≤𝜁𝜈≤𝜁⋅2𝜖1+𝜁𝜖−𝜎=1+𝜁𝜖+𝜎−1()1+𝜁𝜖−𝜎≤1+−21−𝜎)(𝜁𝜎+𝜃(𝜁+𝜎𝜁+2𝜃)(1+𝜁𝜖−𝜎)<1. (36) It follows that 𝐼−𝑋⁻¹Δ𝑋 is nonsingular and 𝐼−𝑋⁻¹≤1Δ𝑋1−𝑋⁻¹≤1Δ𝑋‖‖.1−𝜁Δ𝑋 (37) Using (22) and Lemma 2, we have 𝑍_𝑖(Δ𝑋)≤1−𝑝_𝑖‖Δ𝑋‖²𝑋⁻¹²×𝐼+𝑋⁻¹Δ𝑋⁻¹𝑋^𝑝𝑖≤𝜉_𝑖𝜁²‖Δ𝑋‖²‖‖,𝑉1−𝜁Δ𝑋_𝑖(Δ𝑋)≤𝑋^𝑝𝑖‖Δ𝑋‖𝑋⁻¹𝐼+𝑋⁻¹Δ𝑋⁻¹≤𝜉_𝑖𝜁‖Δ𝑋‖‖‖,1−𝜁Δ𝑋‖ℎ(Δ𝑋)‖≤^𝑚_𝑖=1𝐴_𝑖²𝑍_𝑖(+2𝐴Δ𝑋)_𝑖+Δ𝐴_𝑖×Δ𝐴_𝑖𝑉_𝑖(≤Δ𝑋)_𝑚𝑖=1𝜉_𝑖𝜁²𝐴_𝑖²‖Δ𝑋‖²‖‖+2𝐴1−𝜁Δ𝑋_𝑖+Δ𝐴_𝑖×Δ𝐴_𝑖𝜉_𝑖𝜁‖Δ𝑋‖‖‖.1−𝜁Δ𝑋 (38)(39) Noting (31) and (34), it follows that 𝑓(≤1Δ𝑋)𝑙‖Δ𝑄‖+^𝑚_𝑖=1𝑛_𝑖Δ𝐴_𝑖+𝜁_𝑖𝑙Δ𝐴_𝑖²+1𝑙‖ℎ(𝜎‖‖Δ𝑋)‖≤𝜖+Δ𝑋‖‖+𝜃‖‖1−𝜁Δ𝑋Δ𝑋²1−𝜁‖Δ𝑋‖≤𝜖+𝜎𝜈+1−𝜁𝜈𝜃𝜈²1−𝜁𝜈=𝜈, (40) for Δ𝑋∈Ω. That is, 𝑓(Ω)⊆Ω. According to Schauder fixed point theorem, there exists Δ𝑋_∗∈Ω such that 𝑓(Δ𝑋_∗)=𝑋_∗. It follows that 𝑋+Δ𝑋_∗ is a Hermitian solution of (7). By Lemma 3, we know that the solution of (7) is unique. Then ‖Δ𝑋_∗‖=‖𝑋−𝑋‖≤𝜈.

Remark 9. From Theorem 8, we get the first order perturbation bound for the solution as follows: ≤1𝑋−𝑋𝑙‖Δ𝑄‖+^𝑚_𝑖=1𝑛_𝑖Δ𝐴_𝑖+𝑂(Δ𝐴₁,Δ𝐴₂,…,Δ𝐴_𝑚,Δ𝑄)²_𝐹, (41) as(Δ𝐴₁,Δ𝐴₂,…,Δ𝐴_𝑚,Δ𝑄)→0.
Combining this with (29) gives Δ𝑋=𝐋⁻¹Δ𝑄+𝐋^𝑚−1_𝑖=1𝐵^∗_𝑖Δ𝐴_𝑖+Δ𝐴^∗_𝑖𝐵_𝑖+𝑂Δ𝐴₁,Δ𝐴₂,…,Δ𝐴_𝑚,Δ𝑄²_𝐹. (42) as (Δ𝐴₁,Δ𝐴₂,…,Δ𝐴_𝑚,Δ𝑄)→0.

4. Backward Error

In this section, a backward error of an approximate solution for the unique solution to (1) is obtained.

Theorem 10. Let 𝑋>0 be an approximation to the solution 𝑋 of (1). If Σ=∑^𝑚_𝑖=1𝑝_𝑖‖𝑋^𝑝𝑖/2𝐴_𝑖𝑋^−1/2‖²<1 and the residual 𝑅(𝑋)≡𝑄+∑^𝑚_𝑖=1𝐴^∗_𝑖𝑋^𝑝𝑖𝐴_𝑖−𝑋 satisfies 𝑅𝑋<𝜃₁2𝑋⁻¹𝜃min1,¹2,where𝜃₁≡1+𝑋⁻¹𝑅𝑋−Σ>0, (43) then 𝑋−𝑋≤𝜇𝑅𝑋,where2𝑋𝑋𝜇=⁻¹𝜃₁+𝜃²₁−4𝑋⁻¹𝑅𝑋. (44)

Proof. Let Ψ=Δ𝑋∈ℋ^𝑛×𝑛∶𝑋^−1/2Δ𝑋𝑋^−1/2≤𝜃₂𝑅𝑋, (45) where 𝜃₂=𝜇/‖𝑋‖. Obviously, Ψ is a nonempty bounded convex closed set. Let 𝑔(Δ𝑋)=^𝑚_𝑖=1𝐴^∗_𝑖𝑋+Δ𝑋^𝑝𝑖−𝑋^𝑝𝑖𝐴_𝑖+𝑅𝑋. (46) Evidently 𝑔∶Ψ↦ℋ^𝑛×𝑛 is continuous. We will prove that 𝑔(Ψ)⊆Ψ. For every Δ𝑋∈Ψ, we have 𝑋^−1/2Δ𝑋𝑋^−1/2≤𝜃₂𝑅𝑋. (47) Hence 𝑋^−1/2Δ𝑋𝑋^−1/2≥−𝜃₂𝑅𝑋𝐼. (48) That is, 𝑋+Δ𝑋≥1−𝜃₂𝑅𝑋𝑋. (49) Using (43), one sees that 𝜃₂𝑅𝑋=2𝑋⁻¹𝑅𝑋𝜃₁+𝜃²₁−4𝑋⁻¹𝑅𝑋<2𝑋⁻¹𝑅𝑋𝜃₁<1. (50) Therefore (1−𝜃₂‖𝑅(𝑋)‖)𝑋>0.
According to (17), we obtain 𝑋^−1/2𝑔(Δ𝑋)𝑋^−1/2≤𝑋^−1/2Δ𝑋𝑋^−1/2+𝑋^−1/2Δ𝑋𝑋^−1/221−𝜃₂𝑅𝑋Σ+𝑋^−1/2𝑅𝑋𝑋^−1/2≤𝜃₂𝑅𝑋+𝜃₂𝑅𝑋²1−𝜃₂𝑅𝑋Σ+𝑋⁻¹𝑅𝑋=𝜃₂𝑅𝑋, (51) for Δ𝑋∈Ψ. That is, 𝑔(Ψ)⊆Ψ. By Brouwer’s fixed point theorem, there exists a Δ𝑋∈Ψ such that 𝑔(Δ𝑋)=Δ𝑋. Hence 𝑋+Δ𝑋 is a solution of (1). Moreover, by Lemma 3, we know that the solution 𝑋 of (1) is unique. Then 𝑋−𝑋=‖Δ𝑋‖≤𝑋𝑋^−1/2Δ𝑋𝑋^−1/2=𝜃₂𝑋𝑅𝑋=𝜇𝑅𝑋. (52)

5. Condition Number

In this section, we apply the theory of condition number developed by Rice [36] to study condition numbers of the unique solution to (1).

5.1. The Complex Case

Suppose that 𝑋 and 𝑋 are the solutions of the matrix equations (1) and (7), respectively. Let Δ𝐴_𝑖=𝐴_𝑖−𝐴_𝑖, Δ𝑄=𝑄−𝑄 and Δ𝑋=𝑋−𝑋. Using Theorem 8 and Remark 9, we have Δ𝑋=𝑋−𝑋=𝐋⁻¹Δ𝑄+𝐋^𝑚−1_𝑖=1𝐵^∗_𝑖Δ𝐴_𝑖+Δ𝐴^∗_𝑖𝐵_𝑖+𝑂Δ𝐴₁,Δ𝐴₂,…,Δ𝐴_𝑚,Δ𝑄²_𝐹, (53) as (Δ𝐴₁,Δ𝐴₂,…,Δ𝐴_𝑚,Δ𝑄)→0.

By the theory of condition number developed by Rice [36], we define the condition number of the Hermitian positive definite solution 𝑋 to (1) by 𝑐(𝑋)=lim_𝛿→0sup_{Δ𝐴1/𝜂1,Δ𝐴2/𝜂2,…,Δ𝐴𝑚/𝜂𝑚,Δ𝑄/𝜌𝐹≤𝛿}‖Δ𝑋‖_𝐹,𝜉𝛿 (54) where 𝜉, 𝜌, and 𝜂_𝑖,  𝑖=1,2,…,𝑚, are positive parameters. Taking 𝜉=𝜂_𝑖=𝜌=1 in (54) gives the absolute condition number 𝑐_abs(𝑋), and taking 𝜉=‖𝑋‖_𝐹, 𝜂_𝑖=‖𝐴_𝑖‖_𝐹, and 𝜌=‖𝑄‖_𝐹 in (54) gives the relative condition number 𝑐_rel(𝑋).

Substituting (53) into (54), we get 𝑐(𝑋)=1𝜉max_{(Δ𝐴1/𝜂1,Δ𝐴2/𝜂2,…,Δ𝐴𝑚/𝜂𝑚,Δ𝑄/𝜌)≠0Δ𝐴𝑖∈𝒞𝑛×𝑛,Δ𝑄∈ℋ𝑛×𝑛}𝐋⁻¹Δ𝑄+∑^𝑚_𝑖=1𝐵^∗_𝑖Δ𝐴_𝑖+Δ𝐴^∗_𝑖𝐵_𝑖𝐹Δ𝐴₁/𝜂₁,Δ𝐴₂/𝜂₂,…,Δ𝐴_𝑚/𝜂_𝑚,Δ𝑄/𝜌_𝐹=1𝜉max_{𝐸1,𝐸2,…,𝐸𝑚𝐸,𝐻≠0𝑖∈𝒞𝑛×𝑛,𝐻∈ℋ𝑛×𝑛}𝐋⁻¹𝜌𝐻+∑^𝑚_𝑖=1𝜂_𝑖𝐵^∗_𝑖𝐸_𝑖+𝐸^∗_𝑖𝐵_𝑖𝐹𝐸₁,𝐸₂,…,𝐸_𝑚,𝐻_𝐹. (55)

Let 𝐿 be the matrix representation of the linear operator 𝐋. Then it is easy to see that 𝐿=𝐼+^𝑚_𝑖=1sin𝑝_𝑖𝜋𝜋^∞₀(𝜆𝐼+𝑋)⁻¹𝑋^1/2𝐴_𝑖^𝑇⊗(𝜆𝐼+𝑋)⁻¹𝑋^1/2𝐴_𝑖^∗×𝜆^𝑝𝑖−1𝑑𝜆. (56) Let 𝐿⁻¹𝐿=𝑆+𝐢Σ,₋₁𝐼⊗𝐵^∗_𝑖=𝐿⁻¹𝐼⊗𝑋^𝑝𝑖𝐴_𝑖^∗=𝑈_𝑖1+𝐢Ω_𝑖1,𝐿⁻¹𝐵^𝑇_𝑖⊗𝐼Π=𝐿⁻¹𝑋^𝑝𝑖𝐴_𝑖^𝑇⊗𝐼Π=𝑈_𝑖2+𝐢Ω_𝑖2,𝑆_𝑐=,𝑈𝑆−ΣΣ𝑆_𝑖=𝑈_𝑖1+U_𝑖2Ω_𝑖2−Ω_𝑖1Ω_𝑖1+Ω_𝑖2𝑈_𝑖1−𝑈_𝑖2,𝑖=1,2,…,𝑚, (57)vec𝐻=𝑥+𝐢𝑦,  vec𝐸_𝑖=𝑎_𝑖+𝐢𝑏_𝑖,  𝑔=(𝑥^𝑇,𝑦^𝑇,𝑎^𝑇₁,𝑏^𝑇₁,…,𝑎^𝑇_𝑚,𝑏^𝑇_𝑚)^𝑇,  𝑀=(𝐸₁,𝐸₂,…,𝐸_𝑚,𝐻), where 𝑥,𝑦,𝑎_𝑖,𝑏_𝑖∈ℛ^𝑛2,𝑆,Σ,𝑈_𝑖1,𝑈_𝑖2,Ω_𝑖1,Ω_𝑖2∈ℛ^𝑛2×𝑛2,  𝑖=1,2,…,𝑚,  and  Π is the vec-permutation matrix, such that vec𝐴^𝑇=Πvec𝐴. (58) Then we obtain that 𝑐(𝑋)=1𝜉max_𝑀≠0𝐋⁻¹𝜌𝐻+∑^𝑚_𝑖=1𝜂_𝑖𝐵^∗_𝑖𝐸_𝑖+𝐸^∗_𝑖𝐵_𝑖𝐹𝐸₁,𝐸₂,…,𝐸_𝑚,𝐻_𝐹=1𝜉max_𝑀≠0𝜌𝐿⁻¹+vec𝐻_𝑚𝑖=1𝜂_𝑖𝐿⁻¹𝐼⊗𝐵^∗_𝑖vec𝐸_𝑖+𝐵^𝑇_𝑖⊗𝐼vec𝐸^∗_𝑖×vec𝐸₁,𝐸₂,…,𝐸_𝑚,𝐻⁻¹=1𝜉max_𝑀≠0𝜌()+𝑆+𝐢Σ𝑥+𝐢𝑦_𝑚𝑖=1𝜂_𝑖𝑈_𝑖1+𝐢Ω_𝑖1𝑎_𝑖+𝐢𝑏_𝑖+𝑈_𝑖2+𝐢Ω_𝑖2𝑎_𝑖−𝐢𝑏_𝑖×vec𝐸₁,𝐸₂,…,𝐸_𝑚,𝐻⁻¹=1𝜉max_𝑔≠0𝜌𝑆_𝑐,𝜂₁𝑈₁,𝜂₂𝑈₂,…,𝜂_𝑚𝑈_𝑚𝑔𝑔=1𝜉𝜌𝑆_𝑐,𝜂₁𝑈₁,𝜂₂𝑈₂,…,𝜂_𝑚𝑈_𝑚,𝐸_𝑖∈𝒞^𝑛×𝑛,𝐻∈ℋ^𝑛×𝑛. (59)

Then we have the following theorem.

Theorem 11. If ∑^𝑚_𝑖=1‖𝐴_𝑖‖²/𝛽^1−𝑝𝑖<1, then the condition number 𝑐(𝑋) defined by (54) has the explicit expression 𝑐(𝑋)=1𝜉𝜌𝑆_𝑐,𝜂₁𝑈₁,𝜂₂𝑈₂,…,𝜂_𝑚𝑈_𝑚, (60) where the matrices 𝑆_𝑐 and 𝑈_𝑖 are defined as in (57).

Remark 12. From (60) we have the relative condition number 𝑐_rel(𝑋)=‖𝑄‖_𝐹𝑆_𝑐,𝐴_1𝐹𝑈₁,𝐴_2𝐹𝑈₂,…,𝐴_𝑚F𝑈_𝑚‖𝑋‖_𝐹. (61)

5.2. The Real Case

In this subsection we consider the real case. That is, all the coefficient matrices 𝐴_𝑖, 𝑄 of (1) are real. In such a case the corresponding solution 𝑋 is also real. Completely similar arguments as Theorem 11 give the following theorem.

Theorem 13. Let 𝐴_𝑖, 𝑄 be real and let 𝑐(𝑋) be the condition number defined by (54). If ∑^𝑚_𝑖=1‖𝐴_𝑖‖²/𝛽^1−𝑝𝑖<1, then 𝑐(𝑋) has the explicit expression 𝑐(𝑋)=1𝜉𝜌𝑆_𝑟,𝜂₁𝑈₁,𝜂₂𝑈₂,…,𝜂_𝑚𝑈_𝑚, (62) where 𝑆_𝑟=𝐼+^𝑚_𝑖=1sin𝑝_𝑖𝜋𝜋^∞₀(𝜆𝐼+𝑋)⁻¹𝑋^1/2𝐴_𝑖^𝑇⊗(𝜆𝐼+𝑋)⁻¹𝑋^1/2𝐴_𝑖^𝑇𝜆^𝑝𝑖−1𝑑𝜆⁻¹𝑈_𝑖=𝑆_𝑟𝐼⊗𝐴^𝑇_𝑖𝑋^𝑝𝑖+𝐴^𝑇_𝑖𝑋^𝑝𝑖⊗𝐼Π,𝑖=1,2,…,𝑚. (63)

Remark 14. In the real case the relative condition number is given by 𝑐_rel(𝑋)=‖𝑄‖_𝐹𝑆_𝑟,𝐴_1𝐹𝑈₁,𝐴_2𝐹𝑈₂,…,𝐴_𝑚𝐹𝑈_𝑚‖𝑋‖_𝐹. (64)

6. Numerical Examples

To illustrate the results of the previous sections, in this section three simple examples are given, which were carried out using MATLAB 7.1. For the stopping criterion we take 𝜀_𝑘+1(𝑋)=‖𝑋_𝑘−∑^𝑚_𝑖=1𝐴^∗_𝑖𝑋^𝑝𝑖_𝑘𝐴_𝑖−𝑄‖<1.0𝑒−10.

Example 15. We consider the matrix equation 𝑋−𝐴^∗₁𝑋^1/2𝐴₁−𝐴^∗₂𝑋^1/3𝐴₂=𝐼, (65) with 𝐴₁=(1/3)+2×10⁻²‖𝐴‖𝐴𝐴,₂=(1/6)+3×10⁻²‖𝐴‖.𝐴,𝐴=00.9501 (66)

Suppose that the coefficient matrices 𝐴₁ and 𝐴₂ are perturbed to 𝐴_𝑖=𝐴_𝑖+Δ𝐴_𝑖,𝑖=1,2, where Δ𝐴₁=10^−𝑗𝐶^𝑇𝐶+𝐶^𝑇+𝐶,Δ𝐴₂=3×10^−𝑗−1𝐶^𝑇𝐶+𝐶^𝑇+𝐶, (67) and 𝐶 is a random matrix generated by MATLAB function randn.

We now consider the corresponding perturbation bounds for the solution 𝑋 in Theorems 5 and 8.

The assumptions in Theorem 5 are ass₁=2𝛽−𝑠−𝑏>0,ass₂=𝛽+‖Δ𝑄‖−^𝑚_𝑖=1𝑝_𝑖𝛽^𝑝𝑖𝐴_𝑖²>0,ass₃=𝑏²−4𝛽−𝑠(𝑠+‖Δ𝑄‖)≥0. (68)

The assumptions in Theorem 8 are ass₄=1−𝜎>0,ass₅=(1−𝜎)²𝜁+𝜎𝜁+2𝜃+2((𝜁+𝜃)(𝜎𝜁+𝜃))−𝜖>0. (69) By computation, we list them in Table 1.

Table 1 

Assumptions check for Example 15 with different values of 𝑗.

The results listed in Table 1 show that the assumptions of Theorems 5 and 8 are satisfied.

By Theorems 5 and 8, we can compute the relative perturbation bounds 𝜉₁ and 𝜉₂=𝜈/‖𝑋‖, respectively. These results averaged as the geometric mean of 10 randomly perturbed runs. Some results are listed in Table 2.

Table 2 

Perturbation bounds for Example 15 with different values of 𝑗.

The results listed in Table 2 show that the perturbation bound 𝜉₂ given by Theorem 8 is fairly sharp, while the bound 𝜉₁ given by Theorem 5 which does not depend on the exact solution is conservative.

Example 16. We consider the matrix equation 𝑋−𝐴^∗₁𝑋^0.5𝐴₁−𝐴^∗₂𝑋^0.25𝐴₂=𝐼, (70) with 𝐴₁=(1/3)+2×10⁻²‖𝐴‖𝐴𝐴,₂=(1/6)+3×10⁻²‖𝐴‖.𝐴,𝐴=2100012100012100012100012 (71) Choose 𝑋₁=𝐴, 𝑋₂=2𝐴. Let the approximate solution 𝑋_𝑘 of 𝑋 be given with the iterative method (6), where 𝑘 is the iteration number.

The residual 𝑅(𝑋_𝑘)≡𝐼+𝐴^∗₁𝑋_𝑘^0.5𝐴₁+𝐴^∗₂𝑋_𝑘^0.25𝐴₂−𝑋_𝑘 satisfies the conditions in Theorem 10.

By Theorem 10, we can compute the backward error bound for 𝑋_𝑘 as follows: 𝑋_𝑘−𝑋≤𝜇𝑅𝑋_𝑘, (72) where 2𝑋𝜇=^𝑘𝑋_𝑘⁻¹𝜃₁+𝜃²₁−4𝑋_𝑘⁻¹𝑅𝑋_𝑘,𝜃₁≡1+𝑋_𝑘⁻¹𝑅𝑋_𝑘−0.5𝑋_𝑘^1/4𝐴₁𝑋_𝑘^−1/22+0.25𝑋_𝑘^1/8𝐴₂𝑋_𝑘^−1/22. (73) Some results are listed in Table 3.

Table 3 

Backward error bound for Example 16 with different values of 𝑘.

The results listed in Table 3 show that the error bound given by Theorem 10 is fairly sharp.

Example 17. We study the matrix equation 𝑋−𝐴^∗₁𝑋^1/2𝐴₁−𝐴^∗₂𝑋^1/3𝐴₂=𝑄, (74) with 𝐴₁=00.55+10^−𝑘,𝐴00₂=12𝐴₁,𝑄=.1101 (75) By Remark 14, we can compute the relative condition number 𝑐_rel(𝑋). Some results are listed in Table 4.
Table 4 shows that the unique positive definite solution 𝑋 is well conditioned.

Table 4 

Relative condition number for Example 17 with different values of 𝑘.

Acknowledgments

The author wishes to express her gratitude to the referees for their fruitful comments. The work was supported in part by National Nature Science Foundation of China (11201263), Natural Science Foundation of Shandong Province (ZR2012AQ004), and Independent Innovation Foundation of Shandong University (IIFSDU), China.