A structure-preserving doubling algorithm for the square root of regular M-matrix

1.   Introduction

Let $A\in\mathbb{C}^{n\times n}$, a matrix $X\in\mathbb{C}^{n\times n}$ is called a square root of $A$ if it satisfies the following matrix equation:

The matrix square root is widely encountered in many fields of mathematics, such as nonlinear matrix equations, computation of the matrix logarithm, boundary value problems, Markov chains, and so on (see ^[1,2,3,4]).

In recent years, many scholars have conducted in-depth research on the theories and numerical algorithms of matrix square root, and obtained a large number of results. In terms of theoretical research, people have discussed the existence and the number of square roots for general matrices and obtained many results (see ^[2,3]). In addition, for matrices with special structure and special properties, such as M-matrix (see ^[5]), H-matrix with positive diagonal entries (see ^[6]), central symmetric matrix (see ^[7,8]), P-orthogonal matrix (see ^[9]), Toeplitz matrix (see ^[10]), cyclic matrix (see ^[11]), Boolean matrix (see ^[12]), and so on, the existence of square root is also studied.

In terms of numerical algorithms, the Schur method was first used to calculate the square root of a general matrix (see ^[13,14]). In fact, the command $\tt sqrtm$ in Matlab for calculating the square root of a general matrix was written based on the Schur method. However, when the size of the matrix increases, the complexity of the Schur method increases sharply. Moreover, for some matrices with special structures, iterative methods are more effective, so people turned to iterative methods. Higham first proposed the Newton method to calculate the square root of a general matrix in ^[15]. Although the Newton method has a fast convergence rate, the operation count of each iteration is too large. In order to overcome the defects of the Newton method, the simple Newton method was proposed. However, the simple Newton method has poor stability in calculation (see ^[2]). Later, the Newton method was improved, and some fast and effective algorithms were proposed, such as the DB method, the CR method, and so on (see ^{[16,17,18,19,20]}). Recently, many effective algorithms for calculating matrix square root have also been proposed from different perspectives, such as the power method (see ^[21,22]), gradient neural network method (see ^[23]), Zolotarev iterative method (see ^[24]), Chebyshev iterative method (see ^[25]), and so on (see ^[26,27,28]).

Generally speaking, the theory of matrix square root is very complex, and the existence of a square root of a general matrix is not as evident as it seems. For a general matrix, it may have no square root, or may have many square roots. A sufficient condition for one to exist is that $A$ has no real negative eigenvalue. More generally, any matrix $A$ having no non-positive real eigenvalues has a unique square root for which every eigenvalue has positive real part, and this square root, denoted by $A^{1/2}$, is called the principal square (see ^[1,2]). As for M-matrix, we have the following result.

Lemma 1.1. ^[2] Let $A$ be a nonsingular M-matrix. Then $A$ has exactly one nonsingular M-matrix as its square root.

But for general M-matrix, the above conclusion is not necessarily valid. For example, consider

It is easy to verify that $A$ and $B$ are M-matrices, but direct calculation shows that $A$ has no square root while $B$ has an M-matrix square root. Therefore, it is necessary to extend the above theorem to more general subclasses of M-matrix. The regular M-matrix is an extension of the nonsingular M-matrix and the irreducible singular M-matrix. The regular M-matrix has many beautiful properties similar to the nonsingular M-matrix and the irreducible singular M-matrix, and it also plays an important role in the theories of nonsymmetric algebraic Riccati equations. In this paper, we will prove that for a regular M-matrix, there always exists a regular M-matrix as its square root. In addition, a structure-preserving doubling algorithm is proposed to compute the square root. Theoretical analysis and numerical experiments are given to show that our method is feasible and is effective under certain conditions.

The rest of the paper is organized as follows: In Section 2, we give some preliminary results of M-matrix. In Section 3, based on the properties of M-matrix and quadratic matrix equations, we prove the existence of square root of a regular M-matrix. In Section 4, we propose a structure-preserving doubling algorithm to compute the square root and then give a convergence analysis of it. In Section 5, we use some numerical examples to show the feasibility and effectiveness of our method. Conclusions are given in Section 6.

2.   Preliminaries

In this section, we review the definitions and some properties of M-matrix.

For any matrix $A\in \mathbb{R}^{m\times n}$, if the elements of $A$ satisfy $a_{ij}\geq 0$ for $1\leq i\leq m$, $1\leq j\leq n$, then $A$ is called a nonnegative matrix, denoted by $A\geq0$. For any matrices $A\in\mathbb{R}^{m\times n}$, $B\in\mathbb{R}^{m\times n}$, we write $A\geq B$ if $a_{ij}\geq b_{ij}$ hold for all $1\leq i\leq m$, $1\leq j\leq n$. Let $A = (a_{ij})\in \mathbb{R}^{n\times n}$, then $A$ is called a Z-matrix if $a_{ij}\leq 0$ for all $i\neq j$. A Z-matrix $A$ is called an M-matrix if there exists a nonnegative matrix $B$ such that $A = sI-B$ and $s\geq\rho(B)$, where $\rho(B)$ is the spectral radius of $B$. In particular, $A$ is called a nonsingular M-matrix if $s > \rho(B)$ and is called a singular M-matrix if $s = \rho(B)$.

The following properties of M-matrix are well-known and can be found in ^[29].

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

(i) $A$ is a nonsingular M-matrix;

(ii) $A^{-1}\geq 0$;

(iii) $Au > 0$ for some positive vector $u > 0$.

Lemma 2.2. Let $A$ be a nonsingular M-matrix and $B$ be a Z-matrix. If $A\leq B$, then $B$ is also a nonsingular M-matrix.

Lemma 2.3. Let $A$ be an irreducible singular M-matrix and $B$ be a Z-matrix. If $A\leq B$ and $A\neq B$, then $B$ is a nonsingular M-matrix.

The regular M-matrix is an extension of the nonsingular M-matrix and the irreducible singular M-matrix. The definition of regular M-matrix is introduced in the following:

Definition 2.1. ^[30] An M-matrix $A$ is said to be regular if $Au\geq 0$ for some $u > 0$.

It is easy to verify that nonsingular M-matrices and irreducible singular M-matrices are always regular M-matrices, and any Z-matrix $A$ such that $Au\geq 0$ for some $u > 0$ is a regular M-matrix.

3.   The square root of regular M-matrix

In this section, we consider the square root of a regular M-matrix, and prove that for a regular M-matrix there exists a regular M-matrix as its square root.

Let $A = (a_{ij})\in \mathbb{R}^{n\times n}$ be a regular M-matrix, and let

where $\alpha > 0$ is a parameter to be determined. Then Eq (1.1) can be rewritten as

Here we choose the parameter $\alpha$ such that $\alpha^2I-A$ is a nonnegative matrix. Since $A$ is a regular M-matrix, we can easily verify from definition that $a_{ii}\geq0$ for $i = 1, 2, \cdots, n$, where $a_{ii}$ are the diagonal elements of $A$. In addition, it is evident that $\max\limits_{1\leq i\leq n}\sqrt{a_{ii}} > 0$ unless $A = 0$. Thus, if we choose

then $\alpha^2I-A$ is a nonnegative matrix.

In the following, we discuss the existence of a minimal non-negative solution for Eq (3.2). In order to achieve this goal, we write Eq (3.2) in the following form:

and then consider the iteration

We have the following result for iteration (3.3).

Lemma 3.1. Let $A = (a_{ij})\in\mathbb{R}^{n\times n}$ be a regular M-matrix, and the parameter $\alpha$ satisfy $\alpha\geq\max\limits_{1\leq i\leq n}\sqrt{a_{ii}}$. Then the sequence $\{Y_k\}$ generated by iteration (3.3) is well defined, monotonically increasing, and satisfies $Y_ku\leq \alpha u$, where $u > 0$ is a positive vector such that $Au\geq 0$.

Proof. It is evident that the sequence $\{Y_k\}$ is well defined, since only matrix additions and multiplications are used in the iteration (3.3).

We prove the conclusion by mathematical induction.

(ⅰ) The sequence $\{Y_k\}$ is monotonically increasing, i.e.,

When $k = 0$, it is evident that

Suppose the conclusion is true for $k-1$. From

and

we know the conclusion is true for $k$. So for any $k\geq0$, the conclusion is true.

(ⅱ) Since $A$ is a regular M-matrix, there is a positive vector $u > 0$ such that $Au\geq 0$. In the following, we prove that for any $k\geq0$, $Y_ku\leq \alpha u$ holds true.

When $k = 0$, the conclusion is obvious. Suppose that the conclusion is true for $k-1$. Then

thus, the same is true for $k$. So for any $k\geq0$, the conclusion is true. □

Theorem 3.1. Let $A = (a_{ij})\in\mathbb{R}^{n\times n}$ be a regular M-matrix, and the parameter $\alpha$ satisfy $\alpha\geq\max\limits_{1\leq i\leq n}\sqrt{a_{ii}}$. Then Eq (3.2) has a minimal non-negative solution $Y$, and $\alpha I-Y$ is a regular M-matrix.

Proof. According to Lemma 3.1, the sequence $\{Y_k\}$ obtained from iteration (3.3) is monotonically increasing and upper bounded. So there is a limit $\lim\limits_{k\rightarrow \infty}Y_k = Y$. Taking limit on both sides of (3.3), we know $Y$ is a solution of Eq (3.2), and $Y\geq 0$.

If $Z\geq 0$ is another non-negative solution of Eq (3.2), we can easily verify as in Lemma 3.1 that the sequence $\{Y_k\}$ obtained from iteration (3.3) satisfies $Y_k\leq Z$. Taking limit yields $Y\leq Z$, so $Y$ is the minimal non-negative solution.

In addition, we have obtained $Y_ku\leq \alpha u$ for $k\geq0$ in the proof of Lemma 3.1. Taking limit, we have $Yu\leq \alpha u$, so

Thus, by definition, $\alpha I-Y$ is a regular M-matrix. □

According to the above theorem, we can achieve the following conclusion.

Theorem 3.2. Let $A\in\mathbb{R}^{n\times n}$ be a regular M-matrix. Then there exists a square root of $A$, and the square root is a regular M-matrix.

Proof. By Theorem 3.1, $\alpha I-Y$ is a regular M-matrix, and it is a square root of $A$ since it satisfies Eq (1.1). □

Corollary 3.1. Let $A\in\mathbb{R}^{n\times n}$ be a regular M-matrix, and

be two parameters. Then $Y_1\geq Y_2$, where $Y_1$ and $Y_2$ are minimal non-negative solutions of Eq (3.2) associated with $\alpha_1$ and $\alpha_2$, respectively.

Proof. Let $A^{1/2}$ be the square root of $A$, as in Theorem 3.2. Since

we have $Y_1\geq Y_2$ immediately. □

Corollary 3.2. Let $A\in\mathbb{R}^{n\times n}$ be a regular M-matrix and $Y$ be the minimal non-negative solution of Eq (3.2). If $A$ is a nonsingular M-matrix, then $\rho(Y) < \alpha$; if $A$ is singular, then $\rho(Y) = \alpha$.

Proof. Note that $\alpha I-Y$ is an M-matrix and satisfies Eq (1.1). If $A$ is nonsingular, so is $\alpha I-Y$, and by definition of M-matrix we know $\rho(Y) < \alpha$. If $A$ is singular, so is $\alpha I-Y$, and thus $\rho(Y) = \alpha$. □

Remark 1. We have proved that, for a regular M-matrix $A$, there exists a square root of $A$, and the square root is a regular M-matrix. However, $A$ may have more than one M-matrix as its square root. For example, consider

It is evident that $A$ is a regular M-matrix, and $A$ is a square root of itself since $A^2 = A$. In addition,

is an M-matrix (not regular), and $B$ is also a square root of $A$.

4.   A structure-preserving doubling algorithm

The doubling algorithms are very efficient iterative methods for solving nonlinear matrix equations. The idea of doubling algorithms can be traced back to the 1970s, and recently the doubling algorithms have been successfully applied to various kinds of nonlinear matrix equations. We refer to ^[31] for a complete description of doubling algorithms. In this section, we propose a structure-preserving doubling algorithm to compute the square root of a regular M-matrix and then give a theoretical analysis of it.

The analysis in the previous section shows that in order to calculate the square root of a regular M-matrix, it is only necessary to compute the minimal non-negative solution of Eq (3.2). For effectively solving (3.2), we divide both sides of Eq (3.2) by $\alpha^2$ and let $Z = Y/\alpha$ to obtain

By Theorem 3.1 and Corrolary 3.2, Eq (4.1) has a minimal non-negative solution $Z$ such that $\rho(Z)\leq 1$. In particular, $\rho(Z) < 1$ when $A$ is a nonsingular M-matrix.

It can be verified that Eq (4.1) is equivalent to

Pre-multiply Eq (4.2) by

then we have

where

It is clear that the matrix pencil $M_0-\lambda N_0$ is in the first standard form (SF1) as defined in ^[31], and it is natural to apply the doubling algorithm for SF1 to solve (4.3).

The basic idea of the structure-preserving doubling algorithm for SF1 to solve (4.3) is to recursively construct a sequence of pencils $M_k-\lambda N_k$ for $k\geq1$ that satisfy

and at the same time have the same forms as $M_0-\lambda N_0$:

Moreover, the blocks in $M_k$ and $N_k$ can be produced by the following iteration of doubling algorithms:

For Eq (4.1), another matrix equation can be constructed as follows:

which is called the dual equation of (4.1). It can be verified similarly that Eq (4.6) has a minimal non-negative solution $W$ such that $\rho(W)\leq 1$ as for (4.1). In addition, we have

for matrix pencils $M_k-\lambda N_k$ in (4.5).

Thus, as long as $E_k$, $F_k$, $G_k$, and $H_k$ are well-defined, we will have

or, equivalently,

According to the above derivation, the structure-preserving doubling algorithm to compute the square root of regular M-matrix can be stated as follows:

In the following, we give a convergence analysis of the structure-preserving doubling algorithm.

Lemma 4.1. Let $A\in\mathbb{R}^{n\times n}$ be a regular M-matrix, and $Z$ and $W$ be the minimal non-negative solutions for Eqs (4.1) and (4.6), respectively. Then $I-ZW$ and $I-WZ$ are both regular M-matrices. In particular, $I-ZW$ and $I-WZ$ are nonsingular M-matrices if $A$ is nonsingular, and they are irreducible M-matrices if $A$ is irreducible.

Proof. (ⅰ) Since $A$ is a regular M-matrix, there exists a positive vector $u > 0$ such that $Au\geq 0$. Firstly, we can easily verify as in Lemma 3.1 that the sequences $\{Z_k\}$ and $\{W_k\}$ generated by the following iterations

are well defined, monotonically increasing, and satisfy

for all $k\geq0$. Taking limit, we have $Zu\leq u$ and $Wu\leq u$.

In addition, it is clear $I-ZW$ and $I-WZ$ are Z-matrices, and noting that

we know $I-ZW$ and $I-WZ$ are both regular M-matrices.

(ⅱ) When $A$ is a nonsingular M-matrix, it follows that $Au = v > 0$ for $u > 0$. We can verify further that the sequence $\{Z_k\}$ generated by (4.10) satisfies

Taking limit, we obtain

Thus

which mean $I-ZW$ and $I-WZ$ are nonsingular M-matrices.

When $A$ is irreducible, from (4.10) and (4.11) we can observe

are irreducible. Since $Z\geq Z_1\geq 0$ and $W\geq W_2\geq0$, we can conclude $Z$ and $W$ are irreducible M-matrices. Hence $I-ZW$ and $I-WZ$ are irreducible M-matrices. □

By the general theory of doubling algorithms, we can conclude the following convergence result when $A$ is a nonsingular M-matrix or an irreducible singular M-matrix.

Theorem 4.1. Let $A\in\mathbb{R}^{n\times n}$ be a nonsingular M-matrix or an irreducible singular M-matrix, and $Z$ and $W$ be the minimal non-negative solutions for Eqs (4.1) and (4.6), respectively. Then the sequences $\{E_k\}$, $\{F_k\}$, $\{G_k\}$, $\{H_k\}$ generated by Algorithm 4.1 are well-defined, and for $k\geq 1$

(a) $E_k = (I-G_kZ)Z^{2^k}\geq0$;

(b) $F_k = (I-H_kW)W^{2^k}\geq0$;

(c) $I-G_kH_k$ and $I-H_kG_k$ are nonsingular M-matrices;

(d) $0\leq H_k\leq H_{k+1}\leq Z$, $0\leq G_k\leq G_{k+1}\leq W$, and

In addition, we have

Proof. We prove the conclusions by mathematical induction.

(ⅰ) When $k = 1$, we observe firstly that $E_{0}$, $F_{0}$, $G_{0}$, and $H_{0}$ are well-defined and are all non-negative. In addition, since

both $I-H_0G_0$ and $I-G_0H_0$ are nonsingular M-matrices by Lemmas 2.2 and 2.3. Thus $E_{1}$, $F_{1}$, $G_{1}$, and $H_{1}$ are well-defined. Moreover, from the iteration of doubling algorithms we have

Let $k = 1$ in (4.8) and (4.8) to get

thus we have

From Lemma 4.1 and noting that

we know $I-H_1G_1$ and $I-G_1H_1$ are nonsingular M-matrices. Furthermore, we have

This completes the proof for $k = 1$.

(ⅱ) Suppose now that the conclusions are true for $k = l$. From the iteration of doubling algorithms, $E_{l+1}$, $F_{l+1}$, $G_{l+1}$, and $H_{l+1}$ are well-defined, and satisfy

On the other hand, let $k = l+1$ in (4.8) and (4.8) to obtain

Thus, we have

which mean that $I-H_{l+1}G_{l+1}$ and $I-G_{l+1}H_{l+1}$ are nonsingular M-matrices from Lemma 4.1. Furthermore, we have

Thus, the conclusions are true for $k = l+1$.

By the induction, the conclusions are true for all $k\geq1$.

In addition, from (4.12) we can obtain

Letting $k\rightarrow \infty$, we have

□

Remark 2. We only prove the convergence result of the structure-preserving doubling algorithm for the nonsingular M-matrix and the irreducible singular M-matrix. However, when $A$ is a general regular M-matrix, the convergence analysis of the structure-preserving doubling algorithm is very complicated. We refer to ^[31] for a complete convergence analysis of the structure-preserving doubling algorithm. In addition, from Theorem 4.1, we can observe that when $A$ is a nonsingular M-matrix, the convergence rate is quadratic since $\rho(Z)\cdot\rho(W) < 1$. When $A$ is an irreducible singular M-matrix, the convergence rate is linear.

5.   Numerical examples

In this section, we use some numerical examples to show the feasibility and effectiveness of the structure-preserving doubling algorithm (Algorithm 4.1). We will compare Algorithm 4.1 (SDA) with the basic Newton method (Newton) in ^[2]. The numerical results are presented in terms of iteration numbers (IT), CPU time (CPU) in seconds, and residue (Res), where the residue is defined to be

In our implementations, all iterations are performed in Matlab (R2018a) on a personal computer with a 2G CPU and 8G memory and are terminated when the current iterate satisfies $Res < 10^{-12}$.

Example 5.1. Consider nonsingular M-matrix as follows:

where $n = m^2$. For different orders of $m$, we apply both methods to calculate the square root of $A$. The numerical results are shown in Table 1.

It can be seen from the table that both methods can obtain results satisfying the accuracy for different orders of $m$. In particular, the Newton method is a little faster than SDA, while SDA is much cheaper than the Newton method in terms of CPU time.

Example 5.2. Consider the nonsingular M-matrix generated by the MATLAB command as in the following

For different sizes of $n$, we apply both methods to calculate the square root of $A$. The numerical results are reported in Table 2.

From Table 2, it can be seen that both methods can obtain results satisfying the accuracy. In particular, SDA is much cheaper than the Newton method in terms of CPU time.

Example 5.3. Consider irreducible singular M-matrix as follows:

For different sizes of $n$, the numerical results are reported in Table 3.

In this example, due to the singularity of $A$, the convergence rates of the two iteration methods are linear. In particular, SDA is a little faster than the Newton method and is much cheaper than the Newton method.

Example 5.4. Consider the following regular singular M-matrix

The square root (reserved to four decimal places) is

The numerical results are shown in Table 4. For this example, the Newton method fails to converge, but SDA can still get a result that meets the accuracy.

From the above four examples, it can be seen that the SDA proposed in this paper is feasible. In particular, it requires fewer CPU time than the Newton method.

6.   Conclusions

We studied the square root of M-matrix in this paper and proved that the square root of a regular M-matrix exists and is still a regular M-matrix. In addition, we proposed a structure-preserving doubling algorithm (SDA) to compute the square root. Theoretical analysis and numerical experiments have shown that SDA is feasible and is effective under certain conditions.

Use of AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

Acknowledgments

The authors would like to express their gratitude to the editor and the two anonymous reviewers for their valuable comments and suggestions, which have greatly improved the quality of the paper. This work was supported by the National Natural Science Foundation of China (12001395).

Conflict of interest

The authors declare there is no conflicts of interest.