# A Fast Newton-Shamanskii Iteration for a Matrix Equation Arising from M/G/1-Type Markov Chains

Authors: Pei-Chang Guo
Journal: Mathematical Problems in Engineering
 (2017)

---

## Abstract

For the nonlinear matrix equations arising in the analysis of M/G/1-type and GI/M/1-type Markov chains, the minimal nonnegative solutionG or R can be found by Newton-like methods. We prove monotone convergence results for the Newton-Shamanskii iteration for this class of equations. Starting with zero initial guess or some other suitable initial guess, the Newton-Shamanskii iteration provides a monotonically increasing sequence of nonnegative matrices converging to the minimal nonnegative solution. A Schur decomposition method is used to accelerate the Newton-Shamanskii iteration. Numerical examples illustrate the effectiveness of the Newton-Shamanskii iteration.


---

## Body

## 1. Introduction

Some necessary notation for this article is as follows. For any matrixB=[bij]∈Rn×n, B≥0(B>0) if bij≥0(bij>0) for all i,j; for any matrices A,B∈Rn×n, A≥B(A>B) if aij≥bij(aij>bij) for all i,j; the vector with all entries equal to one is denoted by e; that is, e=(1,1,…,1)T; and the identity matrix is denoted by I. An M/G/1-type Markov Chain (MC) is defined by a transition probability matrix of the form(1)P=B0B1B2B3⋯CA1A2A3⋯A0A1A2⋱A0A1⋱0⋱⋱,while the transition probability matrix of a GI/M/1-type MC is as follows:(2)P=B0C0B1A1A0B2A2A1A0B3A3A2A1⋱⋮⋮⋱⋱⋱,where B0∈Rm0×m0 and A1∈Rm×m, respectively. N is the smallest index i such that Ai, for i>N, is (numerically) zero. The steady-state probability vector of an M/G/1-type MC, if it exists, can be expressed in terms of a matrix G that is the element-wise minimal nonnegative solution to the nonlinear matrix equation [1](3)G=∑i=0NAiGi.Similarly, for the GI/M/1-type MC, a matrix R is of practical interest, which is the element-wise minimal nonnegative solution to the nonlinear matrix equation [2](4)R=∑i=0NRiAi.For a large number of stochastic models, most steady-state and certain transient characteristics of the process can be expressed in terms of the minimal nonnegative solution described above. This problem arises in the analysis of queues with phase type service times, as well as in queues that can be represented as quasi-birth-death processes [3]. It is known that any M/G/1-type MC can be transformed into a GI/M/1-type MC and vice versa through either the Ramaswami [4] or Bright [5] dual, and the G(R) matrix can be obtained directly in terms of the R(G) matrix of the dual chain. The drift of the chain is defined by(5)ρ=pTβ,where p is the stationary probability vector of the irreducible stochastic matrix A=∑i=0NAi and β=∑i=1NiAie. The MC is positive and recurrent if ρ<1, null and recurrent if ρ=1, and transient if ρ>1; and throughout this article it is assumed that ρ≠1.Available algorithms for finding the minimal nonnegative solution to (3) include functional iterations [1], pointwise cyclic reduction (CR) [6], the invariant subspace (IS) approach [7], the Ramaswami reduction (RR) [8], and the Newton iteration (NI) [3, 9–11]. For a detailed comparison of these algorithms, we refer the readers to [11] and the references therein. In [11], Newton’s iteration is revisited and accelerated. From numerical experience, the fast Newton’s iteration in [11] is a very competitive algorithm.In this paper, we consider the Newton-Shamanskii iteration for (3). It is shown that, starting with a suitable initial guess, the sequence generated by the Newton-Shamanskii iteration is monotonically increasing and converges to the minimal nonnegative solution of (3). Similar to Newton’s iteration, equation involved in the Newton-Shamanskii iteration step is a linear equation of the form ∑j=0N-1BjXCj=E, which can be solved fast by a Schur decomposition method. The Newton-Shamanskii iteration differs from Newton’s iteration as the Fréchet derivative is not updated at each iteration; therefore the special coefficient matrix structure form can be reused.The paper is organized as follows. The Newton-Shamanskii iteration and its accelerated iterative procedure using a Schur decomposition method are given in Section2. Then M/G/1-type MCs with low-rank downward transitions and low-rank local and upward transitions are considered in Sections 3 and 4, respectively. Numerical results in Section 6 show that the fast Newton-Shamanskii iteration can be more efficient than the fast Newton’s iteration proposed in [11]. Final conclusions are presented in Section 6.


## 2. Newton-Shamanskii Iteration

In this section, we present the Newton-Shamanskii iteration for (3). First we rewrite (3) as(6)GX=∑v=0NAvXv-X=0.The function G is a mapping from Rm×m into itself and the Fréchet derivative of G at X is a linear map GX′:Rm×m→Rm×m given by(7)GX′Z=∑v=1N∑j=0v-1AvXjZXv-1-j-Z.The second derivative at X, GX′′:Rm×m→Rm×m, is given by(8)GX′′Z1,Z2=∑v=2N∑j=0v-1Av∑i=0j-1XiZ2Xj-1-iZ1Xv-1-j+∑v=2N∑j=0v-2AvXjZ1∑i=0v-2-jXiZ2Xv-2-j-i.For a given initial guessG0,0, the Newton-Shamanskii iteration for the solution of G(x)=0 is as follows.Fork=0,1,…,(9)GGk,0′Xk,s-1=-GGk,s-1,Gk,s=Gk,s-1+Xk,s-1,s=1,2,…,nk,(10)Gk+1=Gk+1,0=Gk,nk.Xk,s-1 is the solution to(11)Xk,s-1-∑v=1N∑j=0v-1AvGkjXk,s-1Gkv-1-j=∑v=0NAvGk,s-1v-Gk,s-1,which, after rearranging the terms, can be rewritten as(12)Xk,s-1-∑j=0N-1∑v=j+1NAvGkv-1-jXk,s-1Gkj=∑v=0NAvGk,s-1v-Gk,s-1.Following the notation of [11], we define Sk,i=∑j=iNAjGkj-i; then the above equation is(13)Sk,1-IXk,s-1+∑j=1N-1Sk,j+1Xk,s-1Gkj=Gk,s-1-∑v=0NAvGk,s-1v,which is a linear equation of the same form ∑j=0N-1BjXCj=E as the Newton’s iteration step. It can be solved fast by applying a Schur decomposition on the matrix C, which is the m×m matrix Gk here, and then solving m linear systems with m unknowns and equations. For the detailed description for solving ∑j=0N-1BjXCj=E, we refer the reader to [11, 12]. We stress that, for Newton-Shamanskii iteration, the coefficient matrices are updated once after every nk iteration step and the special coefficient structure can be reused, so the cost per iteration step is reduced significantly.


## 3. The Case of Low-Rank Downward Transitions

When the matrixA0 is of rank r, meaning that it can be decomposed as A0=A^0Γ with A^0∈Rm×r and Γ∈Rr×m, we refer to the MC as having low-rank downward transitions. If Newton-Shamanskii iteration is applied to this case, all the matrices Xk,s-1 can be written as X^k,s-1Γ. This can be shown by making induction on the index s. X0,0 can be written as X^0,0Γ and we assume that it is true for all Xl,j-1 for l=0,…,k and j=1,…,s-1. Hence, Gk,s-1 can be written as G^k,s-1Γ, since Gk,s-1=∑l=0k-1∑j=1nlXl,j-1+∑j=1s-1Xk,j-1=(∑l=0k-1∑j=1nlX^l,j-1+∑j=1s-1X^k,j-1)Γ. Then (12) can be rewritten as(14)Xk,s-1=A^0Γ+∑j=1NAjGk,s-1j-1G^k,s-1Γ-G^k,s-1Γ+∑v=1NAvGkv-1Xk,s-1+∑j=1N-1∑v=j+1NAvGkv-1-jXk,s-1Gkj-1G^kΓ=I-∑v=1NAvGkv-1-1×A^0+∑j=1NAjGk,s-1j-1G^k,s-1-G^k,s-1∑j=1N-1∑v=j+1NAvGkv-1-jXk,s-1Gkj-1G^kΓ;therefore Xk,s-1 can be decomposed as the product of an m×r matrix X^k,s-1 and an r×m matrix Γ. The inverse on the right-hand side exists since 0≤∑v=1NAvGkv-1≤∑v=1NAvGv-1 and the spectral radius of ∑v=1NAvGv-1   is strictly less than one [13]. Therefore we will concentrate on finding X^k,s-1 as the solution to(15)X^k,s-1=A^0+∑j=1NAjGk,s-1j-1-IG^k,s-1+∑v=1NAvGkv-1X^k,s-1+∑j=1N-1∑v=j+1NAvGkv-1-jX^k,s-1ΓGkj-1G^k=A^0+∑j=1NAjGk,s-1j-1-IG^k,s-1+∑j=0N-1Sk,j+1X^k,s-1ΓG^kj,which can be rewritten as(16)Sk,1-IX^k,s-1+∑j=1N-1Sk,j+1X^k,s-1ΓG^kj=I-∑j=1NAjGk,s-1j-1G^k,s-1-A^0.We can use the Schur decomposition method in [11, 12] to solve the above equation. Different from Newton’s iteration in [11], the special coefficient structure can be reused here, thus saving the overall computational cost. We will report the numerical performance of the Newton-Shamanskii iteration in Section 6.


## 4. The Case of Low-Rank Local and Upward Transitions

In this section, the case of low-rank local and upward transitions is considered, where them×m matrices {Ai,1≤i≤N} can be decomposed as Ai=ΓA^i with Γ∈Rm×r and A^i∈Rr×m. To exploit low-rank local and upward transitions, we introduce the matrix U, which is the generator of the censored Markov chain on level i, starting from level i, before the first transition on level i-1. The following equality holds based on a level crossing argument:(17)U=∑i=1NAiGi-1=∑i=1NAiI-U-1A0i-1.For the case of low-rank local and upward transitions, we can rewrite U as(18)U=∑i=1NAiI-U-1A0i-1=Γ∑i=1NA^iI-U-1A0i-1=ΓU^,which means that U is of rank r, while G=(I-U)-1A0 is generally of rank m.Therefore, we findU as the solution to(19)FX=X-∑i=1NAiI-X-1A0i-1=0and get G from G=(I-U)-1A0 [11, 14]. The Newton-Shamanskii iteration step for (19) is as follows.Fork=0,1,…,(20)FUk′Yk,s-1=-FUk,s-1,Uk,s=Uk,s-1+Yk,s-1,s=1,2,…,nk,Uk+1=Uk+1,0=Uk,nk.Yk,s-1 is the solution to(21)Yk,s-1-∑i=2NAi∑j=1i-1I-Uk-1A0j-1I-Uk-1Yk,s-1I-Uk-1A0i-j=∑i=1NAiI-Uk,s-1-1A0i-1-Uk,s-1.If we define Rk,j=∑i=j+1NAi((I-Uk)-1A0)i-1-j(I-Uk)-1 and rearrange the terms, (21) can be rewritten as(22)Yk,s-1-∑j=1N-1Rk,jYk,s-1I-Uk-1A0j=∑i=1NAiI-Uk,s-1-1A0i-1-Uk,s-1,which is of the form ∑j=0N-1BjXCj=E. This iteration enables us to exploit low-rank local and upward transitions. The iterates Uk,s=Uk,s-1+Yk,s-1, where Yk,s-1 solves (21), can be rewritten as Uk,s=ΓU^k,s. This can be shown by making induction on the index s. It obviously holds for U0,0. Assuming that Uk,s-1=ΓU^k,s-1, from (21), we get(23)Yk,s-1=Γ∑i=2NA^i∑j=1i-1I-Uk-1A0j-1I-Uk-1Yk,s-1I-Uk-1A0i-j+∑i=1NA^iI-Uk,s-1-1A0i-1-U^k,s-1,which tells us that Yk,s-1 can be decomposed as ΓY^k,s-1, and the same holds for Uk,s=Uk,s-1+Yk,s-1. Therefore, from (21), we will focus on finding Y^k,s-1 as the solution to(24)Y^k,s-1-∑i=2NA^i∑j=1i-1I-Uk-1A0j-1I-Uk-1Yk,s-1I-Uk-1A0i-j=∑i=1NA^iI-Uk,s-1-1A0i-1-U^k,s-1.DefiningR^k,j=∑i=j+1NA^i((I-Uk)-1A0)i-1-j(I-Uk)-1Γ, we can rewrite the above equation as(25)Y^k,s-1-∑j=1N-1R^k,jY^k,s-1I-Uk-1A0j=∑i=1NA^iI-Uk,s-1-1A0i-1-U^k,s-1,which is of the form ∑j=0N-1BjXCj=E.


## 5. Convergence Analysis

There is monotone convergence when the Newton-Shamanskii method is applied to (3).


### 5.1. Preliminary

Let us first recall that a real square matrixA is a Z-matrix if all its off-diagonal elements are nonpositive and can be written as sI-B with B≥0. Moreover, a Z-matrix A is called an M-matrix if s≥ρ(B), where ρ(·) is the spectral radius; it is a singular M-matrix if s=ρ(B) and a nonsingular M-matrix if s>ρ(B). The following result from [15] is to be exploited.Lemma 1.
For aZ-matrix A, the following statements are equivalent:(a)
A is a nonsingular M-matrix.(b)
A-1≥0.(c)
Av>0 for some vector v>0.(d)
All eigenvalues ofA have positive real parts.The following result is also well known [15].Lemma 2.
LetA be a nonsingular M-matrix. If B≥A is a Z-matrix, then B is a nonsingular M-matrix. Moreover, B-1≤A-1.The minimal nonnegative solutionS for (3) may also be recalled—see [9] for details.Theorem 3.
If the rateρ defined by (5) satisfies ρ≠1, then the matrix(26)I-∑v=1N∑j=0v-1Gv-1-jT⨂AvGjis a nonsingular M-matrix.


### 5.2. Monotone Convergence

The following lemma displays the monotone convergence properties of the Newton iteration for (3).Lemma 4.
Consider a matrixX such that(i)
G(X)≥0,(ii)
0≤X≤G,(iii)
I-∑v=1N∑j=0v-1(Xv-1-j)T⊗AvXj is a nonsingular M-matrix.Then the matrix(27)Y=X-GX′-1GXis well defined, and(a)
G(Y)≥0,(b)
0≤X≤Y≤G,(c)
I-∑v=1N∑j=0v-1(Yv-1-j)T⊗AvYj is a nonsingular M-matrix.Proof.
G
X
′ is invertible and the matrix Y is well defined from (iii) and Lemma 1. Since(28)I-∑v=1N∑j=0v-1Xv-1-jT⊗AvXj-1≥0from (iii) and Lemma 1 and G(X)≥0, we get that vec(Y)≥vec(X) and thus Y≥X. From (27) and the Taylor formula, there exists a number θ1, 0<θ1<1, such that(29)GY=GX+GX′Y-X+12GX′′θ1Y-X,θ1Y-X=12GX′′θ1Y-X,θ1Y-X≥0,so (a) is proven, where the last inequality is obtained from Y-X≥0 and (8). (b) may be proven as follows. From(30)0=GG=GX+GX′G-X+12GX′′θ2G-X,θ2G-X,where 0<θ2<1, we have(31)-GX′G-Y=GX′Y-X-GX′G-X=-GX-GX′G-X=12GX′′θ2G-X,θ2G-X≥0,where the last inequality is from G-X≥0 by (ii). It is notable that(32)I-∑v=1N∑j=0v-1Xv-1-jT⨂AvXjis a nonsingular M-matrix, so vec(G-Y)≥0 from Lemma 1; that is, G-Y≥0. Now Y≥X, so (b) follows. Next we prove (c). From 0≤Y≤G, we have(33)I-∑v=1N∑j=0v-1Yv-1-jT⨂AvYj≥I-∑v=1N∑j=0v-1Gv-1-jT⨂AvGj,and we know that I-∑v=1N∑j=0v-1(Gv-1-j)T⊗AvGj is a nonsingular M-matrix. Consequently, from Lemma 2, I-∑v=1N∑j=0v-1(Yv-1-j)T⊗AvYj is a nonsingular M-matrix.A generalization of Lemma4 provides the theoretical basis for the monotone convergence of the Newton-Shamanskii method for (3).Lemma 5.
Consider a matrixX such that (i)
G(X)≥0,(ii)
0≤X≤G,(iii)
I-∑v=1N∑j=0v-1(Xv-1-j)T⊗AvXj is a nonsingular M-matrix.Then, for any matrix Z, where 0≤Z≤X, the matrix(34)Y=X-GZ′-1GXexists, such that(a)
G(Y)≥0,(b)
0≤X≤Y≤G,(c)
I-∑v=1N∑j=0v-1(Yv-1-j)T⊗AvYj is a nonsingular M-matrix.Proof.
Since0≤Z≤X, we have(35)I-∑v=1N∑j=0v-1Zv-1-jT⨂AvZj≥I-∑v=1N∑j=0v-1Xv-1-jT⨂AvXj.From (iii) and Lemma 2, GZ′ is invertible and the matrix Y is well defined such that 0≤X≤Y.Let(36)Y^=X-GX′-1GX;we know that Y^≥Y from Lemma 2. Because Y^≤G from Lemma 4, (b) follows. Now(37)I-∑v=1N∑j=0v-1Y^v-1-jT⨂AvY^jis a nonsingular M-matrix from Lemma 4 and Y^≥Y; therefore I-∑v=1N∑j=0v-1(Yv-1-j)T⊗AvYj is a nonsingular M-matrix from Lemma 2. Next we show that (a) is true. From the Taylor formula, there exists two numbers θ3 and θ4, where 0<θ3,θ4<1, such that(38)GY=GX+GX′Y-X+12GX′′θ3Y-X,θ3Y-X=GX+GZ′Y-X+GX′-GZ′Y-X+12GX′′θ3Y-X,θ3Y-X=GZ′′Y-X,θ4X-Z+12GX′′θ3Y-X,θ3Y-X≥0,where the latter inequality holds, since X-Z≥0 and Y-X≥0.The monotone convergence result for the Newton-Shamanskii method applied to (3) follows.Theorem 6.
Suppose that a matrixG0 is such that(i)
G(G0)≥0,(ii)
0≤G0≤G,(iii)
I-∑v=1N∑j=0v-1(G0v-1-j)T⊗AvG0j is a nonsingular M-matrix.Then the Newton-Shamanskii algorithm (9)-(10) generates a sequence Gk such that Gk≤Gk+1≤G for all k≥0, and limk→∞Gk=G.Proof.
The proof is by mathematical induction. From Lemma5,(39)G0=G0,0≤⋯≤G0,n0=G1≤G,GG1≥0,I-∑v=1N∑j=0v-1G1v-1-jT⨂AvG1jis a nonsingular M-matrix. Assuming that(40)GGi≥0,G0=G0,0≤⋯≤G0,n0=G1≤⋯≤Gi-1,ni-1=Gi≤Gand that I-∑v=1N∑j=0v-1(Giv-1-j)T⊗AvXij is a nonsingular M-matrix, from Lemma 5,(41)GGi+1≥0,Gi=Gi,0≤⋯≤Gi,ni=Gi+1≤G,and I-∑v=1N∑j=0v-1(Gi+1v-1-j)T⊗AvGi+1j is a nonsingular M-matrix. By induction, the sequence Gk is therefore monotonically increasing and bounded above by G and so has a limit G∗ such that G∗≤G. Letting i→∞ in Gi+1≥Gi,1=Gi-(GGi′)-1G(Gi)≥0, it follows that G(G∗)=0. Consequently, G∗=G since G∗≤G and G is the minimal nonnegative solution of (3).


## 5.1. Preliminary

Let us first recall that a real square matrixA is a Z-matrix if all its off-diagonal elements are nonpositive and can be written as sI-B with B≥0. Moreover, a Z-matrix A is called an M-matrix if s≥ρ(B), where ρ(·) is the spectral radius; it is a singular M-matrix if s=ρ(B) and a nonsingular M-matrix if s>ρ(B). The following result from [15] is to be exploited.Lemma 1.
For aZ-matrix A, the following statements are equivalent:(a)
A is a nonsingular M-matrix.(b)
A-1≥0.(c)
Av>0 for some vector v>0.(d)
All eigenvalues ofA have positive real parts.The following result is also well known [15].Lemma 2.
LetA be a nonsingular M-matrix. If B≥A is a Z-matrix, then B is a nonsingular M-matrix. Moreover, B-1≤A-1.The minimal nonnegative solutionS for (3) may also be recalled—see [9] for details.Theorem 3.
If the rateρ defined by (5) satisfies ρ≠1, then the matrix(26)I-∑v=1N∑j=0v-1Gv-1-jT⨂AvGjis a nonsingular M-matrix.


## 5.2. Monotone Convergence

The following lemma displays the monotone convergence properties of the Newton iteration for (3).Lemma 4.
Consider a matrixX such that(i)
G(X)≥0,(ii)
0≤X≤G,(iii)
I-∑v=1N∑j=0v-1(Xv-1-j)T⊗AvXj is a nonsingular M-matrix.Then the matrix(27)Y=X-GX′-1GXis well defined, and(a)
G(Y)≥0,(b)
0≤X≤Y≤G,(c)
I-∑v=1N∑j=0v-1(Yv-1-j)T⊗AvYj is a nonsingular M-matrix.Proof.
G
X
′ is invertible and the matrix Y is well defined from (iii) and Lemma 1. Since(28)I-∑v=1N∑j=0v-1Xv-1-jT⊗AvXj-1≥0from (iii) and Lemma 1 and G(X)≥0, we get that vec(Y)≥vec(X) and thus Y≥X. From (27) and the Taylor formula, there exists a number θ1, 0<θ1<1, such that(29)GY=GX+GX′Y-X+12GX′′θ1Y-X,θ1Y-X=12GX′′θ1Y-X,θ1Y-X≥0,so (a) is proven, where the last inequality is obtained from Y-X≥0 and (8). (b) may be proven as follows. From(30)0=GG=GX+GX′G-X+12GX′′θ2G-X,θ2G-X,where 0<θ2<1, we have(31)-GX′G-Y=GX′Y-X-GX′G-X=-GX-GX′G-X=12GX′′θ2G-X,θ2G-X≥0,where the last inequality is from G-X≥0 by (ii). It is notable that(32)I-∑v=1N∑j=0v-1Xv-1-jT⨂AvXjis a nonsingular M-matrix, so vec(G-Y)≥0 from Lemma 1; that is, G-Y≥0. Now Y≥X, so (b) follows. Next we prove (c). From 0≤Y≤G, we have(33)I-∑v=1N∑j=0v-1Yv-1-jT⨂AvYj≥I-∑v=1N∑j=0v-1Gv-1-jT⨂AvGj,and we know that I-∑v=1N∑j=0v-1(Gv-1-j)T⊗AvGj is a nonsingular M-matrix. Consequently, from Lemma 2, I-∑v=1N∑j=0v-1(Yv-1-j)T⊗AvYj is a nonsingular M-matrix.A generalization of Lemma4 provides the theoretical basis for the monotone convergence of the Newton-Shamanskii method for (3).Lemma 5.
Consider a matrixX such that (i)
G(X)≥0,(ii)
0≤X≤G,(iii)
I-∑v=1N∑j=0v-1(Xv-1-j)T⊗AvXj is a nonsingular M-matrix.Then, for any matrix Z, where 0≤Z≤X, the matrix(34)Y=X-GZ′-1GXexists, such that(a)
G(Y)≥0,(b)
0≤X≤Y≤G,(c)
I-∑v=1N∑j=0v-1(Yv-1-j)T⊗AvYj is a nonsingular M-matrix.Proof.
Since0≤Z≤X, we have(35)I-∑v=1N∑j=0v-1Zv-1-jT⨂AvZj≥I-∑v=1N∑j=0v-1Xv-1-jT⨂AvXj.From (iii) and Lemma 2, GZ′ is invertible and the matrix Y is well defined such that 0≤X≤Y.Let(36)Y^=X-GX′-1GX;we know that Y^≥Y from Lemma 2. Because Y^≤G from Lemma 4, (b) follows. Now(37)I-∑v=1N∑j=0v-1Y^v-1-jT⨂AvY^jis a nonsingular M-matrix from Lemma 4 and Y^≥Y; therefore I-∑v=1N∑j=0v-1(Yv-1-j)T⊗AvYj is a nonsingular M-matrix from Lemma 2. Next we show that (a) is true. From the Taylor formula, there exists two numbers θ3 and θ4, where 0<θ3,θ4<1, such that(38)GY=GX+GX′Y-X+12GX′′θ3Y-X,θ3Y-X=GX+GZ′Y-X+GX′-GZ′Y-X+12GX′′θ3Y-X,θ3Y-X=GZ′′Y-X,θ4X-Z+12GX′′θ3Y-X,θ3Y-X≥0,where the latter inequality holds, since X-Z≥0 and Y-X≥0.The monotone convergence result for the Newton-Shamanskii method applied to (3) follows.Theorem 6.
Suppose that a matrixG0 is such that(i)
G(G0)≥0,(ii)
0≤G0≤G,(iii)
I-∑v=1N∑j=0v-1(G0v-1-j)T⊗AvG0j is a nonsingular M-matrix.Then the Newton-Shamanskii algorithm (9)-(10) generates a sequence Gk such that Gk≤Gk+1≤G for all k≥0, and limk→∞Gk=G.Proof.
The proof is by mathematical induction. From Lemma5,(39)G0=G0,0≤⋯≤G0,n0=G1≤G,GG1≥0,I-∑v=1N∑j=0v-1G1v-1-jT⨂AvG1jis a nonsingular M-matrix. Assuming that(40)GGi≥0,G0=G0,0≤⋯≤G0,n0=G1≤⋯≤Gi-1,ni-1=Gi≤Gand that I-∑v=1N∑j=0v-1(Giv-1-j)T⊗AvXij is a nonsingular M-matrix, from Lemma 5,(41)GGi+1≥0,Gi=Gi,0≤⋯≤Gi,ni=Gi+1≤G,and I-∑v=1N∑j=0v-1(Gi+1v-1-j)T⊗AvGi+1j is a nonsingular M-matrix. By induction, the sequence Gk is therefore monotonically increasing and bounded above by G and so has a limit G∗ such that G∗≤G. Letting i→∞ in Gi+1≥Gi,1=Gi-(GGi′)-1G(Gi)≥0, it follows that G(G∗)=0. Consequently, G∗=G since G∗≤G and G is the minimal nonnegative solution of (3).


## 6. Numerical Experiments

The Newton-Shamanskii iteration differs from Newton’s method in that the evaluation of the Fréchet derivative is not done at every iteration step. So, while more iterations will be needed than for Newton’s method, the overall cost of the modified Newton method could be less. Our numerical experiments confirm the feasibility of the Newton-Shamanskii iteration for (6).About how to choose the optimal scalarsni in the Newton-like algorithm (see (9) and (10)), now we have no theoretical results. This is a goal for our future research. In our extensive numerical experiments, we update the Fréchet derivative every two iteration steps. That is, for i=0,1,…, we choose ni=2 in the Newton-like algorithm (9).The elapsed CPU time in seconds (denoted as “time") is used to measure the feasibility of our new method. In our numerical experiments, we use zero matrix as the initial iteration value and choose the following stopping criterion:(42)Gnew-Gold<1e-14,where Gnew and Gold denote the iteration values of solution after and before one iteration step, respectively, and · denotes the infinity-norm of a matrix. The numerical tests were performed on a laptop (2.4 GHz and 2 G Memory) with MATLAB R2013b. Numerical experiments show that the modified Newton method could be more efficient than Newton iteration proposed in [11]. We present the numerical results for a random-generated problem in Figure 1. In this figure, we fix the number of coefficient matrices (N+1) and vary the problem size (n) and plot the CPU time of the two algorithms in seconds for different parameters n.Figure 1
Comparisons about CPU time.The MATLAB code used for the problem construction is reported as follows, which generates theN+1 coefficient matrices for (6), that is, Ai∈Rn×n,  i=0,1,2,…,N:function [A,er,rho]=problem(n,N)rand(’state’,0);A=rand(n,n∗(N+1));e=ones(n∗(N+1),1);s=A∗e;e=ones(n,1);for i=1:nA(i,:)=A(i,:)/s(i);endB=zeros(n);beita=0;for i=1:(N+1)B=B+A(:,((i-1)∗n+1):(i∗n));beita=beita+(i-1)∗A(:,((i-1)∗n+1):(i∗n))∗e;end[V,D]=eig(B’);pai=V(:,1);s=sum(pai);pai=pai/ser=max(B’∗pai-pai)rho=pai’∗beitaendUp to now, we do not have theoretical results about when this modified Newton method can outperform the Newton iteration. This is not an easy question. We leave it as a goal for future research.



---

*Source: 4018239-2017-10-19.xml*