A structure-preserving doubling algorithm for solving a class of quadratic matrix equation with $M$-matrix

1.   Introduction

In this paper, we consider the problem of finding the maximal nonpositive solvent of the following quadratic matrix equation (QME)

where

Such QME arises from an overdamped vibrating system ^[1,2]. By left multiplying $\widetilde{A}^{-1}$ ^[3], without changing the $M$-matrix structure of it, QME (1.1) can be reduced to the following form

where $B$ is a nonsingular $M$-matrix and $C$ is an $M$-matrix such that $B^{-1}C \ge 0$. It is known that Eq (1.2) has a maximal nonpositive solvent $\varPhi$ under the condition that ^[3]

This solvent $\varPhi$ is the one of interest.

Various iterative methods have been developed to obtain the maximal nonpositive solvent of QME (1.2) with assumption (1.3), including the Newton's method and Bernoulli-like methods (fixed-point iterative methods) ^[3], modified Bernoulli-like methods with diagonal update skill ^[4]. Newton's method is not competitive in terms of CPU time because it is to solve a generalized Sylvester matrix equation in each Newton's iterative step. The fixed-point iterative methods are usually linearly or sublinearly convergent and sometimes can be very slow ^[3].

There are many researches on iterative methods for other QME s. For example, the cyclic reduction algorithm and the invariant subspace method for solving the quadratic matrix equation arising from quasi-birth-death problems ^[5,6]; the methods for solving the quadratic matrix equation from quadratic eigenvalue problems ^{[7,8,9,10,11,12]}; the fixed-point iteration and the Schur method for solving the quadratic matrix equation arising in noisy Wiener-Hopf problems for Markov chains ^[13,14] and others; see ^{[15,16,17,18]} and the references therein. Our work here is mainly inspired by recent study on highly accurate structure-preserving doubling algorithm for quadratic matrix equation from quasi-birth-and-death process ^[19]. Structure-preserving doubling algorithms are very efficient iterative methods for solving nonlinear matrix equations. For instance, some structure-preserving doubling algorithms are presented to solve continuous-time algebraic Riccati equations (ARE) ^[20], periodic discrete-time ARE ^[21], nonsymmetric ARE ^[22,23] and $M$-matrix ARE ^[24,25]. For more applications, please refer to ^[26,27] and the references therein.

Yu et al. ^[3] proved $\rho(\varPhi)\le 1$ under (1.3). In this paper, we will slightly improve their result and prove that $\rho(\varPhi) < 1$ under the same condition. The property $\rho(\varPhi) < 1$ is important since it is desired for the quadratic convergence of structure-preserving doubling algorithms. Based on the new property $\rho(\varPhi) < 1$, furthermore, we extend the structure-preserving doubling algorithm for the first standard form (SF1) ^[27] to solve QME (1.2) and give the quadratically convergent result.

The rest of this paper is organized as follows. In Section 2 we give some notations and state a few basic results on nonnegative and $M$-matrices. The main results of this paper are presented in Section 3. Numerical examples are given in Section 4 to demonstrate the performance of our method. Finally, conclusions are made in Section 5.

2.   Notations and preliminaries

In this section, we first introduce some necessary notations and terminologies for this paper. ${\mathbb{R}}^{m\times m}$ is the set of all $m\times m$ real matrices, ${\mathbb{R}}^n = {\mathbb{R}}^{n\times 1}$, and ${\mathbb{R}} = {\mathbb{R}}^1$. $I_n$ (or simply $I$ if its dimension is clear from the context) is the $n\times n$ identity matrix. For $X\in{\mathbb{R}}^{m\times n}$, $X_{(i, j)}$ refers to its $(i, j)$th entry. Inequality $X\le Y$ means $X_{(i, j)}\le Y_{(i, j)}$ for all $(i, j)$, and similarly for $X < Y$, $X\ge Y$, and $X > Y$. In particular, $X\ge 0$ means that $X$ is entrywise nonnegative and it is called a nonnegative matrix. $X$ is entrywise nonpositive if $-X$ is entrywise nonnegative. A matrix $A\in{\mathbb{R}}^{m\times n}$ is positive, denoted by $A > 0$, if all its entries are positive. The same understanding goes to vectors. For a square matrix $X$, denote by $\rho(X)$ its spectral radius. A matrix $A\in{\mathbb{R}}^{n\times n}$ is called a $Z$-matrix if $A_{(i, j)}\le 0$ for all $i\ne j$. Any $Z$-matrix $A$ can be written as $sI-N$ with $N\ge 0$, and it is called an $M$-matrix if $s\ge \rho(N)$. Specifically, it is a singular $M$-matrix if $s = \rho(N)$, and a nonsingular $M$-matrix if $s > \rho(N)$.

The following results on nonnegative matrices and $M$-matrices can be found in, e.g., ^[28,29].

Theorem 2.1. Let $A\in {\mathbb{R}}^{n\times n}$ be a nonnegative matrix. Then the spectral radius, $\rho(A)$, is an eigenvalue of $A$ and there exists a nonnegative right eigenvector ${\boldsymbol{x}}$ associated with the eigenvalue $\rho(A)$: $A{\boldsymbol{x}} = \rho(A){\boldsymbol{x}}$.

Theorem 2.2. Let $A\in {\mathbb{R}}^{n\times n}$ be a $Z$-matrix. Then the following statements are equivalent:

${\rm (a)}$ $A$ is a nonsingular $M$-matrix;

${\rm (b)}$ $A^{-1}\ge 0$;

${\rm (c)}$ $A{\boldsymbol{u}} > 0$ holds for some positive vector ${\boldsymbol{u}}\in {\mathbb{R}}^n$.

Theorem 2.3. Let $A\in {\mathbb{R}}^{n\times n}$ be an $M$-matrix. Let $B\in {\mathbb{R}}^{n\times n}$ be a $Z$-matrix. If $A$ is nonsingular and $B\ge A$, then $B$ is also a nonsingular $M$-matrix.

3.   The main results

In this section, we give the main results of this paper. Lemma 3.1 below can be found in [3,Theorem 3.1]. The first goal of this paper is to further prove that $\rho(\varPhi) < 1$.

Lemma 3.1. Suppose $(1.3)$, then QME $(1.2)$ has a maximal nonpositive solvent $\varPhi$ with $\rho(\varPhi)\le 1$, also $B+\varPhi$ and $B+\varPhi-C$ are both nonsingular $M$-matrices.

Since $B - C - I$ is a nonsingular $M$-matrix, by Theorems 2.2, there exists a vector $0 < {\boldsymbol{u}} \in {\mathbb{R}}^n$ such that

Unless stated otherwise, throughout the rest of this paper, ${\boldsymbol{u}}$ and ${\boldsymbol{v}}$ are reserved for the ones here. The following lemma is inspired by [19,Lemma 3.2], we still give the proof for completeness.

Lemma 3.2. Suppose $(1.3)$, i.e., $B - C - I$ is a nonsingular $M$-matrix.Then $\rho(X)\neq 1$ for any nonpositive solvent $X$ of Eq $(1.2)$.

Proof. Suppose, to the contrary, that $\rho(X) = 1$ (which is equivalent to $\rho(-X) = 1$), where $X$ is a nonpositive solvent of Eq (1.2). Then according to Theorem 2.1, there exists a nonzero and nonnegative vector ${\boldsymbol{z}}\in {\mathbb{R}}^n$ such that $-X{\boldsymbol{z}} = {\boldsymbol{z}}$ and thus

implies that $(B-C-I){\boldsymbol{z}} = 0$, which contradicts with the fact that $B-C-I$ is nonsingular.

Combining Lemmas 3.1 and 3.2, we immediately finish our first goal of this paper. Moreover, we have the following theorem. The theorem is implied by [19,Theorem 3.1] or [30,Theorem 2.3].

Theorem 3.3. Under the assumption $(1.3)$, QME $(1.2)$ has a unique maximal nonpositive solvent $\varPhi$.Moreover, it holds that $\varPhi\le X_0$ and

where $X_0 = -B^{-1}C$ is as defined in Eq $(3.3\text{a})$.

It can be checked that Theorem 3.3 is applicable to

which is called the dual equation of Eq (1.2). In particular, under assumption (1.3), the dual equation Eq (3.1) also has a unique maximal nonpositive solvent, denoted by $\varPsi$ hereafter. In conclusion, Theorem 3.4 below gives some of the important results, the proof is similar to that of [19,Theorem 3.4] and thus it is omitted here.

Theorem 3.4. Suppose $(1.3)$. The following statements hold.

${\rm (a)}$ We have

${\rm (b)}$ $\rho(\varPhi) < 1$ and $\rho(\varPsi) < 1$.

${\rm (c)}$ $I-\varPhi\varPsi$ and $I-\varPsi\varPhi$ are nonsingular $M$-matrices.

Now we are in position to develop a new structure-preserving doubling algorithm for solving the QME (1.2). Similar to the discussion in the introduction of ^[19], QME (1.2) is connected with the matrix pencil

where

Now that the matrix pencil ${\mathscr{A}}_0-\lambda{\mathscr{B}}_0$ is in (SF1), it is natural for us to apply the doubling algorithm (see Algorithm 1) for (SF1) ^[27] to solve Eq (3.2).

The basic idea of the structure-preserving doubling algorithm for (SF1) ^[23,27] for solving Eq (3.2) is to recursively construct a sequence of matrix pencils ${\mathscr{A}}_i-\lambda{\mathscr{B}}_i$ for $i\ge 1$ that have the same block structure as ${\mathscr{A}}_0-\lambda{\mathscr{B}}_0$:

and at the same time

where $M = X$.

We observe that as long as $E_k, \, F_k, \, X_k, \, Y_k$ are well-defined (so are ${\mathscr{A}}_k$ and ${\mathscr{B}}_k$), we will have

where ${\mathscr{A}}_k$ and ${\mathscr{B}}_k$ are defined as in Eq (3.5). Or, equivalently,

In the following we will analysis the convergence of Algorithm 1 for solving the QME (1.2) under the assumption (1.3). Theorem 3.5 below is essentially [19,Theorem 6.1] or [23,Theorem 4.1]. The only difference lies in the initial matrices $(E_0, F_0, X_0, Y_0)$.

Theorem 3.5. Under $(1.3)$, the matrix sequences $\{E_k\}, \{F_k\}, \{X_k\}$and $\{Y_k\}$ generated by Algorithm 1 are well-defined and, moreover, for $k\ge 1$,

${\rm (a)}$ $E_k = (I-Y_k\varPhi)\varPhi^{2^k}\ge 0$;

${\rm (b)}$ $F_k = (I-X_k\varPsi)\varPsi^{2^k}\ge 0$;

${\rm (c)}$ $I-X_kY_k$ and $I-Y_kX_k$ are nonsingular $M$-matrices;

${\rm (d)}$ $\varPhi\le X_{k}\le X_{k-1}\le 0, \; \varPsi\le Y_k\le Y_{k-1}\le 0$, and

Proof. We prove the theorem by mathematical induction.

Since $B$ is a nonsingular $M$-matrix, we immediately conclude that $X_0$, $Y_0$, $E_0$ and $F_0$ are well-defined as in Eq (3.3) and they are all nonpositive. From Theorem 3.4(a), we obtain that $\varPhi\le X_0\le 0$, $\varPsi\le Y_0\le 0$. Therefore,

By Theorem 3.4(c), both $I-\varPhi\varPsi$ and $I-\varPsi\varPhi$ are nonsingular $M$-matrices; so are $I-X_0Y_0$ and $I-Y_0X_0$ according to Theorem 2.3. Hence, the matrices $E_1, \, X_1, \, F_1, \, Y_1$ generated by the Algorithm 1 are well-defined. Moreover, from Eq (3.4) we have

Let $k = 1$ in Eq (3.6), we have

Noting that $F_1, \, E_1\ge 0$ and $\varPhi, \, \varPsi\le 0$, it follows from Eq (3.8) that $\varPhi\le X_1\le X_0\le 0$, $\varPsi\le Y_1\le Y_0\le 0$. By the same reasoning above, we can conclude that $I-X_1Y_1$ and $I-Y_1X_1$ are nonsingular $M$-matrices. Furthermore, it follows from Eq (3.8), $\varPhi\le X_1\le 0$ and $\varPsi\le Y_1\le 0$ that

This completes the proof of our results for $k = 1$.

Next, suppose that the results hold for all positive integers $k\le\ell$. Hence $E_{\ell+1}, \, X_{\ell+1}, \, F_{\ell+1}, \, Y_{\ell+1}$ are well-defined by Eq (3.4), which, together with the induction hypothesis, guarantee that

On the other hand, Eqs (3.9) and (3.6) for $k = \ell+1$ say

Thus we have

Following the same line as the proof of the $k = 1$ case, we conclude that $I-X_{\ell+1}Y_{\ell+1}$ and $I-Y_{\ell+1}X_{\ell+1}$ are nonsingular $M$-matrices. Similarly, we deduce from Eq (3.10) that

By the induction principle, we have finished the proof.

From Eq (3.7) and Theorem 3.4(b), we can conclude that $X_k$ and $Y_k$ generated by Algorithm 1 converge quadratically to $\varPhi$ and $\varPsi$, respectively, under assumption (1.3).

4.   Numerical examples

In this section, we will present numerical results applying Algorithm 1 to solve QME (1.2). We will compare Algorithm 1 (referred to as DA) with two Bernoulli-like methods presented in ^[3] (referred to, respectively, as BL1 and BL2 as in ^[4]) and three modified Bernoulli-like methods with diagonal update skill ^[4] (referred to as BL1-DU, BL2-DU1 and BL2-DU2, respectively). In reporting numerical results, we will record the numbers of iterations (denoted by "Iter"), the elapsed CPU time in seconds (denoted as "CPU") and plot iterative history curves for normalized residual NRes defined by

All runs terminate if the current iteration satisfies either ${\mbox{NRes}} < 10^{-12}$ or the number of the prescribed iteration $k_{\max} = 1000$ is exceeded. According to DA, we set $X_0 = -B^{-1}C$ for all methods. All experiments are implemented in MATLAB R2018b with a machine precision $2.22\times 10^{-16}$ on a PC Windows 10 operating system with an Intel i7-9700 CPU and 8GB RAM.

Example 1 (^[4]). Consider the Eq $(1.2)$ with

In Table 1, we record the numerical results for Example 1. We find that DA uses the smallest number of iterations and delivers the lowest value of NRes within all the tested methods. For this example, however, DA is not the fastest one in terms of the elapsed CPU time. The reason is that each step of the DA iteration is expensive than that of the other methods and the iteration number of DA can not compensate its additional cost such that DA beats the other methods. Figure 1 plots the convergent history for Example 1. Quadratic monotonic convergence of DA and monotonic linear convergence of Bernoulli-like methods clearly show.

Example 2 (^[4]). Consider the Eq (1.2) with

Table 2 displays the the numerical results for Example 2. We find that DA is the best one for this example in terms of Iter, CPU and NRes. The iteration number of DA compensates its additional cost such that DA beats the other methods in terms of the elapsed CPU time. Figure 2 shows the convergent history for Example 2. Quadratic monotonic convergence of DA and monotonic linear convergence of Bernoulli-like methods again clearly show.

5.   Conclusions

The structure-preserving doubling algorithm for (SF1) ^[27] is extended to compute the maximal nonpositive solvent of a type of QME s. It is shown that the approximations generated by the algorithm are globally monotonically and quadratically convergent. Two numerical examples are presented to demonstrate the feasibility and effectiveness of our method. Our work here can be seen as a new application of the structure-preserving doubling algorithm for (SF1).

Acknowledgments

The author would like to thank Prof. Ren-Cang Li from University of Texas at Arlington for his useful guidance. The author is grateful to both reviewers for their helpful comments and suggestions which have helped to improve the paper. The author is partially supported by the National Natural Science Foundation of China (Grant No. 11901024) and the Natural Science Foundation of Fujian Province (Grant No. 2021J01661).

Conflict of interest

The author declares there is no conflicts of interest.