GUAN Jinrui(關(guān)晉瑞), SHAO Rongxia(邵榮俠), ZUBAIR Ahmed

(1.Department of Mathematics, Taiyuan Normal University, Jinzhong 030619, China;2.School of Statistics and Data Science, Xinjiang University of Finance and Economics,Urumqi 830012, China; 3.Institute of Mathematics and Computer Science,University of Sindh, Jamshoro 71000, Pakistan)

Abstract: Matrix square root plays an important role in many applications of mathematics.In this paper, we consider the problem of computing the square root of M-matrix and propose an iteration method for computing the square root of regular M-matrices.We first transform this problem into an M-matrix algebraic Riccati equation (MARE) and then propose an efficient method to solve this special MARE.Theoretical analysis shows that our method is convergent under certain conditions and numerical experiments are given to show that our method is feasible and outperforms the Binomial iteration.

Key words: Matrix square root;Regular M-matrix;Algebraic Riccati equation;Binomial iteration

1.Introduction

LetA=(aij)∈Rn×n, a matrixX ∈Rn×nsatisfies the equation

is called a square root ofA.The square root of a matrix play an important role in many applications of mathematics, such as nonlinear matrix equations, computation of the matrix logarithm, the boundary value problems, Markov chains, and so on.[6,13-15,20]Research on the theory and efficient numerical methods for the matrix square root is a hot topic in recent years.[1,4,6,9-14,16-20]

The existence of square root of a general matrix is not evident as it seems.For a general matrixA, it may or may not have a square root.A sufficient condition for one to exist is thatAhave no real negative eigenvalue.More generally, any matrixAhaving no nonpositive real eigenvalues has a unique square root for which every eigenvalue has positive real part, and this square root, denoted byA1/2, is called the principal square.[13,15]

WhenAis an M-matrix, Alefeld and Schneider proved the following results.

Theorem 1.1[1]LetAbe an M-matrix.ThenAhas an M-matrix as a square root if and only ifAhas “property c”.In particular,Ahas exactly one M-matrix as a square root if 0 is at most a simple eigenvalue ofA.

The computation of square root of a matrix has been studied for many years by several authors.[13-14]Efficient numerical methods for computing the matrix square root include Schur method[4,11], matrix sign function method[13], Newton method[10], DB method[5], Newton-Schulz method[13], CR method[20]and so on.These methods can be divided into two classes,direct methods and iterative methods.Schur method and matrix sign function method are direct methods, while others are iterative methods.In this paper we only consider iterative methods.

For computing the square root of an M-matrixA=s(I -C), a simple iteration called Binomial iteration was proposed in [1, 13] as follows:

Convergence analysis in [13] shows that Binomial iteration enjoys monotonic convergence, is structure-preserving, and has linear convergence rate for nonsingular M-matrix withA1/2=and limk→∞Pk=P.However, the Binomial iteration is too slowly to converge in some cases.

Other iterative method for computing the square root of M-matrix include the Newton method and Newton-Schulz method, which are both quadratically convergent.However, the Newton method is not stable in general and the Newton-Schulz method is not structurepreserving for M-matrix.

In this paper, we will propose a new iteration method for computing the square root of M-matrix, which is structure-preserving and has linear convergence rate for nonsingular M-matrix but is faster than the Binomial iteration.

The rest of the paper is organized as follows.In Section 2, we review some basic results of M-matrix and M-matrix algebraic Riccati equation.In Section 3, we propose the new iteration method for computing the square root of regular M-matrix and give its convergence analysis.In Section 4, we use some numerical examples to show the effectiveness of the new method.Conclusions is given in Section 5.

2.Preliminaries

In this section, we review some basic results of M-matrix and M-matrix algebraic Riccati equation.

LetA= (aij)∈Rn×n.Ifaij ≤0 for all, thenAis called a Z-matrix.A Z-matrixAis called an M-matrix if there exists a nonnegative matrixBsuch thatA=sI -Bands ≥ρ(B) whereρ(B) is the spectral radius ofB.In particular,Ais called a nonsingular M-matrix ifs ＞ρ(B) and is called a singular M-matrix ifs=ρ(B).

Definition 2.1[8]An M-matrixAis said to be regular ifAv ≥0 for somev ＞0.

The following lemmas can be found in [2].

Lemma 2.1LetA ∈Rn×nbe a Z-matrix.Then the following statements are equivalent:

1)Ais a nonsingular M-matrix;

2)A-1≥0;

3)Av ＞0 for some vectorsv ＞0;

4) All eigenvalues ofAhave positive real part.

Lemma 2.2LetA,Bbe Z-matrices.IfAis a nonsingular M-matrix andA ≤B,thenBis also a nonsingular M-matrix.In particular, for any nonnegative real numberα,B=αI+Ais a nonsingular M-matrix.

Lemma 2.3LetAbe an M-matrix,B ≥Abe a Z-matrix.IfAis nonsingular or irreducible singular withAB, thenBis also a nonsingular M-matrix.

Lemma 2.4LetA,Bbe nonsingular M-matrices andA ≤B, thenA-1≥B-1.

The M-matrix algebraic Riccati equation (MARE) is of the form

whereA,B,CandDare real matrices of sizesm×m,m×n,n×mandn×n, respectively and

is an M-matrix.MARE appears in many branches of applied mathematics, such as transport theory, Markov chains, stochastic process, and so on.See [3, 7] and the references therein for details.For the MARE (2.1), the solution of practical interest is its minimal nonnegative solution.The following basic result is obtained in [8].

Lemma 2.5[8]For the MARE (2.1), ifKin (2.2) is a regular M-matrix, then the equation (2.1) has a unique minimal nonnegative solutionSandD -CSis a regular Mmatrix.

Efficient numerical methods for solving the MARE(2.1)include the Schur method,fixedpoint iteration, Newton iteration, doubling algorithms and etc.For details, see [3] and the references therein.

3.A New Iteration Method

In this section, we will propose a new iteration method to compute the square root of a regular M-matrix.

LetX=αI-Y, whereα ＞0 is a parameter to be determined.Then the equation (1.1)can be transformed into

Ifα2I-Ais a nonnegative matrix, then equation (3.1) is a special type of MARE.In fact, if we chooseαsuch that

then we can easily verify thatα2I-Ais a nonnegative matrix.

The matrixKin (2.2) associated with the MARE (3.1) is

Lemma 3.1LetA=(aij)∈Rn×nbe a regular M-matrix and the parameterαsatisfy(3.2).ThenKdefined in (3.3) is a regular M-matrix.

ProofWe can easily verify thatKin(3.3)is a Z-matrix.SinceAis a regular M-matrix,there is a vectorv ＞0 such thatAv ≥0.Thus we have

By definition,Kis a regular M-matrix.

Theorem 3.1LetA ∈Rn×nbe a regular M-matrix and the parameterαsatisfy (3.2).Then

(i) The MARE (3.1) has a unique minimal nonnegative solutionSαandαI -Sαis a regular M-matrix;

(ii)Ahas an M-matrix square rootX=αI -Sα.

In the following,we propose a new iteration method to compute the minimal nonnegative solution of equation (3.1).

Write the equation (3.1) as a fixed-point form

(2αI -Y)Y=α2I -A,

then we have the following iterations:

Compared with the Binomial iteration,the new iteration method is a little expensive,since at each iteration the cost is 8/3n3operations,while it is 2n3in the Binomial iteration.However,the new iteration method will need less iteration numbers than the Binomial iteration.So it is feasible.

In the following, we give convergence analysis of the new iteration method (3.4).

11. Our Lady washes the Christ-child s little shirts, and wants to dry them: Tatar writes of the saying, The insertion of a reference to the Madonna anchors the tale in a culture where weather was described in religious terms (152). Return to place in story.

Theorem 3.2LetA ∈Rn×nbe a regular M-matrix and the parameterαsatisfy (3.2).Then the sequence{Yk}generated by (3.4) is well defined, converges toS, and satisfy

whereSis the unique minimal nonnegative solution of equation (3.1).

ProofWe first prove (3.5) by induction.

Whenk= 0, we haveY1= (2α)-1(α2I - A).SinceSis the minimal nonnegative solution of (3.1), it holds that (2αI -S)S=α2I -A.By Theorem 3.1, we knowαI -Sis a regular M-matrix, and hence 2αI -Sis a nonsingular M-matrix by Lemma 2.2.ThusS=(2αI -S)-1(α2I -A).It is evident (2αI)-1≤(2αI -S)-1, thusY1≤S.

Suppose that the assertions (3.5) hold fork=l-1, i.e., 0≤Yl-1≤Yl, Yl ≤S.

SinceYl ≤S, we know 2αI -Yl ≥2αI -Sis a nonsingular M-matrix by Lemma 2.3.ThusYl+1is well defined.From 0≤Yl-1≤Yl ≤Sand Lemma 2.4,we know(2αI-Yl-1)-1≤(2αI-Yl)-1≤(2αI -S)-1.We have

and

Hence the assertions (3.5) hold fork=l+1.Thus we have proved by induction that the assertions (3.5) hold for allk ≥0.

Since{Yk}is nonnegative, monotonically increasing and bounded from above, there is a nonnegative matrixYsuch that limk→∞Yk=Y.From (3.5) we knowY ≤S.On the other hand, take the limit in (3.4), we knowYis a solution of (3.1), thusS ≤Y.HenceS=Y.

Theorem 3.3LetA ∈Rn×nbe a regular M-matrix and the parameterαsatisfy (3.2).Then the convergence rate of (3.4) is given by

whereλ=ρ(S).

ProofWe have

Hence

Taking limit on both side and noting thatρ(A)=limk→∞, we have

It is easy to verify thatρ((2αI -S)-1) =, whereλ=ρ(S) is the Perron eigenvalue of the nonnegative matrixS.Thus the conclusion (3.6) holds.

Corollary 3.1IfAis a nonsingular M-matrix, then the convergence rate of (3.4) is linear.IfAis a regular singular M-matrix, then the convergence rate of (3.4) is sublinear.

ProofSince (αI -S)2=A, ifAis nonsingular,αI -Sis a nonsingular M-matrix.Thus we haveα ＞ρ(S) and

The convergence rate of (3.4) is linear in this case.

On the other hand, ifAis singular,αI-Sis a singular M-matrix.By definition we haveα=ρ(S).Thus

the convergence rate of (3.4) is sublinear in this case.

Corollary 3.2IfAis a regular M-matrix with max1≤i≤n{aii} ＞0, then the optimal parameter of (3.4) is given byαopt=max1≤i≤n

ProofLetα1＞0,α2＞0 be two different parameters andα1＞α2.The convergence factor of (3.4) are

respectively, whereλ1=ρ(S1),λ2=ρ(S2) andS1,S2are the minimal nonnegative solutions of (3.1) respectively.Sinceα1I -S1=α2I -S2=A1/2, we haveS1= (α1-α2)I+S2andλ1=(α1-α2)+λ2.Thus

and the conclusion follows.

4.Numerical Experiments

In this section, numerical experiments are given to verify our theoretical analysis.We compare our new method,Linear iteration(LI)in(3.4),with Binomial iteration(BI)in(1.2).We present numerical results of each experiment in terms of iteration numbers (IT), CPU time (CPU) in seconds and residue (Res), where the residue is defined to be Res=In our implementations all iterations are run in MATLAB (R2012a) on a personal computer,and are terminated when the current iterate satisfies‖Xk-1-Xk‖∞/‖Xk-1‖∞＜10-6.

Example 4.1In the first example, we consider the nonsingular M-matrix

which is from discretization of the two-point boundary value problem[21].The square root ofAis unique.For different sizes ofn, the numerical results are concluded in Tab.4.1.

Tab.4.1 Numerical Results of Example 4.1

Example 4.2In this example, we consider the square root of a nonsingular M-matrix generated as follows:

wheree=(1,1,··· ,1)T.For different sizes ofn, the numerical results are concluded in Tab.4.2.

Tab.4.2 Numerical Results of Example 4.2

Example 4.3In this example, we consider an irreducible singular M-matrix

In this case,the square root ofAis also unique and both methods have sublinear convergence rate.Numerical results are reported in Tab.4.3.

Tab.4.3 Numerical Results of Example 4.3

From Tabs.4.1-4.3, we can conclude that the new method is feasible.In addition, it is effective than Binomial iteration in iteration numbers and CPU times.

5.Conclusions

We have proposed a new iteration method for computing the square root of M-matrix.The new method is structure-preserving and has linear convergence rate for nonsingular Mmatrix.Theoretical analysis and numerical experiments showed that it is feasible and efficient than the Binomial iteration.However, since the new method is only linearly convergent, it will be very slowly as compared with quadratic convergent methods.Hence, for the singular case or near-singular case, the quadratic convergent methods will be more preferable.