An inexact Noda iteration for computing the smallest eigenpair of a large irreducible monotone matrix
aa r X i v : . [ m a t h . NA ] M a y NUMERICAL LINEAR ALGEBRA WITH APPLICATIONS
Numer. Linear Algebra Appl. An inexact Noda iteration for computing the smallest eigenpair ofa large irreducible monotone matrix
Ching-Sung Liu
SUMMARYIn this paper, we present an inexact Noda iteration with inner-outer iterations for finding the smallesteigenvalue and the associated eigenvector of an irreducible monotone matrix. The proposed inexact Nodaiteration contains two main relaxation steps for computing the smallest eigenvalue and the associatedeigenvector, respectively. These relaxation steps depend on the relaxation factors, and we analyze how therelaxation factors in the relaxation steps affect the convergence of the outer iterations. By considering twodifferent relaxation factors for solving the inner linear systems involved, we prove that the convergenceis globally linear or superlinear, depending on the relaxation factor, and that the relaxation factor alsoinfluences the convergence rate. The proposed inexact Noda iterations are structure preserving and maintainthe positivity of approximate eigenvectors. Numerical examples are provided to illustrate that the proposedinexact Noda iterations are practical, and they always preserve the positivity of approximate eigenvectors.Copyright c (cid:13)
Received . . .
KEY WORDS: Inexact Noda iteration, Modified inexact Noda iteration, M -matrix, non-negative matrix,monotone matrix, smallest eigenpair, singular value, Perron vector, Perron root. Department of Applied Mathematics, National Chiao Tung University, Taiwan.Ching-Sung Liu, Department of Applied Mathematics, National Chiao Tung University, Taiwan.E-mail: [email protected] 1. INTRODUCTIONMonotone matrices arise in many areas of mathematics, such as stability analysis [19], and boundsfor eigenvalues and singular values [3, 4]. In many applications, one is interested in finding thesmallest eigenvalue l and the associated eigenvector x of an irreducible nonsingular monotone matrix A ∈ R n × n . The smallest eigenvalue l of a monotone matrix A is defined as s min ( A ) = min {| l | | l ∈ s ( A ) } , where s ( A ) denotes the set of eigenvalues of A . In [23, 25], a real matrix A is called monotone if and only if A − is a non-negative matrix. The irreducible nonsingular M -matrices are one of the most important classes of matrices for applications such as discretized PDEs,Markov chains [2] and electric circuits [24], and they have been studied extensively in the literature[5, Chapter 6]. It is well known that there exist some monotone matrices that are not M -matrices,such as matrices that can be written as a product of M -matrices.There are some differences between an M -matrix and a monotone matrix. For example, an M -matrix can be expressed in the form s I − B with a non-negative matrix B and some constant s > r ( B ) , where r ( · ) denotes the spectral radius, see [5]. Thus, the smallest eigenvalue l of anirreducible nonsingular M -matrix A is equal to s − r ( B ) >
0. In contrast, the smallest eigenvalueof a monotone matrix A can only be expressed as s min ( A ) = r ( A − ) − . However, the smallest eigenvalue retains the same properties [12, p. 487], that is, the largest eigenvalue of an irreduciblenon-negative matrix A − is the Perron root, which is simple and equal to the spectral radius of A − with a positive associated eigenvector. Copyright c (cid:13)
Prepared using nlaauth.cls [Version: 2010/05/13 v2.00]
C.-S. LIU
For the computation of the Perron vector of a non-negative matrix B , many methods exist[21, 22, 28, 20, 14, 13, 17, 6, 26, 16] but the power methods are not structure preserving andcannot guarantee the desired positivity of approximations when the Perron vector x has very smallcomponents. Therefore, a central concern is how to preserve strict positivity of approximations to thePerron vector. In 1971, Noda introduced an inverse iteration method with shifted Rayleigh quotient-like approximations [18] for non-negative matrix eigenvalue problems. This iteration method iscalled Noda iteration (NI), and it has also been adapted to the computation of the smallest eigenvalueand the eigenvector of an irreducible nonsingular M -matrix [29, 1]. The major advantages ofNoda iteration are structure preservation and global convergence. More precisely, it generates amonotonically decreasing sequence of approximate eigenvalues that is guaranteed to converge to r ( B ) , and maintains the positivity of approximate eigenvectors. Furthermore, the convergence hasbeen proven to be superlinear [18] and asymptotically quadratic [9]. In [15], the authors introducedtwo inexact strategies for Noda iteration, which are called inexact Noda iteration (INI) to find thePerron vector of a non-negative matrix (or M -matrix). The proposed INI algorithms are practical,and they always preserve the positivity of approximate eigenvectors. Moreover, the convergenceof INI with these two strategies is globally linear and superlinear with convergence order + √ ,respectively. In this paper, we propose an inexact Noda iteration (INI) to find the smallest eigenvalue andthe associated eigenvector of an irreducible monotone matrix A . The major contribution of thispaper is to provide two main relaxation steps for computing the smallest eigenvalue l and theassociated eigenvector x , respectively. The first step is to use O ( g k min ( x k )) as a stopping criterionfor inner iterations, with 0 < g k <
1, where x k is the current positive approximate eigenvector. Thesecond step is to update the approximate eigenvalues by using the recurrence relations l k + = l k − ( − g k ) min (cid:16) x k y k + (cid:17) , where y k + is the next normalized positive approximate eigenvector, soresulting INI algorithms are structure preserving and globally convergent. The above parameter g k is called the “relaxation factor”. We then establish a rigorous convergence theory of INI with twodifferent relaxation factors g k , and prove that the convergence of the resulting INI algorithms isglobally linear, and superlinear with the relaxation factor g k as the convergence rate, respectively.In fact, the inner iterations of INI (or NI) require the solution of ill-conditioned linear systemswhen the sequence of approximate eigenvalues converges to r ( A − ) (or r ( B ) ). In order to reducethe condition number of the inner linear system, we propose a modified Noda iteration (MNI) byusing rank one update for the inner iterations, and we show that MNI and NI are mathematicallyequivalent. For monotone matrix eigenvalue problems, we also develop an integrated algorithm thatcombines INI with MNI, and we call this modified inexact Noda iteration (MINI). This hybriditerative method can significantly improve the condition number of inner linear systems of INI.The paper is organized as follows. In Section 2, we introduce the Noda iteration and somepreliminaries. Section 3 contains the new strategy for inexact Noda iteration, and proves somebasic properties for it. In Section 4, we establish its convergence theory, and derive the asymptotic convergence factor precisely. In Section 5, we present the integrated algorithm that combines INIwith MNI. Finally, in Section 6 we present some numerical examples illustrating the convergencetheory and the effectiveness of INI, and we make some concluding remarks in Section 7.2. PRELIMINARIES AND NOTATIONFor any real matrix B = [ b i j ] ∈ R n × n , we write B ≥ ( > ) if b i j ≥ ( > ) for all 1 ≤ i , j ≤ n .We define | B | = [ | b i j | ] . If B ≥
0, we say B is a non-negative matrix, and if B >
0, we say B is apositive matrix. For real matrices B and C of the same size, if B − C is a non-negative matrix, wewrite B ≥ C . A non-negative (positive) vector is similarly defined. A non-negative matrix B is saidto be reducible if it can be placed into block upper-triangular form by simultaneous row/columnpermutations; otherwise it is irreducible. If m is not an eigenvalue of B , the function sep ( m , B ) isdefined as sep ( m , B ) = k ( m I − B ) − k − . (1) Copyright c (cid:13)
Numer. Linear Algebra Appl. (2010)
Prepared using nlaauth.clsnlaauth.cls
Prepared using nlaauth.clsnlaauth.cls
DOI: 10.1002/nlaNEXACT NODA ITERATION ∠ ( w , z ) denotes the acute angle of any two nonzero vectors w and z . Throughout the paper, we usea 2-norm for vectors and matrices, and the superscript T denotes its transpose.We review some fundamental properties of non-negative matrices, monotone matrices and M -matrices. Definition 1
A matrix A is said to be “monotone” if Ax ≥ x ≥ Theorem 1 ([7]) A is monotone if and only if A is non-singular and A − ≥ . Definition 2
A monotone matrix M is an M -matrix if M = ( m i j ) , m i j ≤ i = j . Lemma 1 ([5])Let M is a nonsingular M -matrix. Then the following statements are equivalent:1. M = ( a i j ) , a i j ≤ i = j , and M − ≥ M = s I − B with some B ≥ s > r ( B ) .For a pair of positive vectors v and w , definemax (cid:16) wv (cid:17) = max i w ( i ) v ( i ) ! , min (cid:16) wv (cid:17) = min i w ( i ) v ( i ) ! , where v = [ v ( ) , v ( ) , . . ., v ( n ) ] T and w = [ w ( ) , w ( ) , . . ., w ( n ) ] T . The following lemma gives boundsfor the spectral radius of a non-negative matrix B . Lemma 2 ([12, p. 493])Let B be an irreducible non-negative matrix. If v > is not an eigenvector of B , thenmin (cid:18) B vv (cid:19) < r ( B ) < max (cid:18) B vv (cid:19) . (2) The Noda iteration [18] is an inverse iteration shifted by a Rayleigh quotient-like approximation ofthe Perron root of an irreducible non-negative matrix B .Given an initial vector x > with k x k =
1, the Noda iteration (NI) consists of three steps: ( b l k I − B ) y k + = x k , (3) x k + = y k + / k y k + k , (4) b l k + = max (cid:18) B x k + x k + (cid:19) . (5)The main task is to compute a new approximation x k + to x by solving the inner linear system(3). From Lemma 2, we know that b l k > r ( B ) if x k is not a scalar multiple of the eigenvector x . Thisresult shows that b l k I − B is an irreducible nonsingular M -matrix, and its inverse is an irreduciblenon-negative matrix. Therefore, we have y k + > and x k + > , i.e., x k + is always a positive vectorif x k is. After variable transformation, we get b l k + from the following relation: b l k + = b l k − min (cid:18) x k y k + (cid:19) , so { b l k } is monotonically decreasing. Copyright c (cid:13)
Numer. Linear Algebra Appl. (2010)
Prepared using nlaauth.clsnlaauth.cls
Prepared using nlaauth.clsnlaauth.cls
DOI: 10.1002/nla
C.-S. LIU
The inexact Noda iteration. Based on the Noda iteration, in [15] the authors propose an inexactNoda iteration (INI) for the computation of the spectral radius of a non-negative irreducible matrix B . In this paper, since A is a monotone matrix, i.e., A − is a non-negative matrix, we replace B by A − in (3), i.e., ( b l k I − A − ) y k + = x k . (6)When A is large and sparse, we see that we must resort to an iterative linear solver to get anapproximate solution. In order to reduce the computational cost of (6), we solve y k + in (6) byinexactly satisfying ( b l k I − A − ) y k + = x k + A − f k , (7)which is equivalent to ( b l k A − I ) y k + = A x k + f k , (8) x k + = y k + / k y k + k , (9) where f k is the residual vector between ( b l k A − I ) y k + and A x k . Here, the residual norm (innertolerance) x k : = k f k k can be changed at each iterative step k . Theorem 2 ([15])Let A be an irreducible monotone matrix and 0 ≤ g < x k > , if x k = x and f k in (8) satisfies (cid:13)(cid:13) A − f k (cid:13)(cid:13) ≤ g min ( x k ) , (10)then { b l k } is monotonically decreasing and lim k → ¥ b l k = r ( A − ) . Moreover, the convergence of INIis at least globally linear.Based on (8)-(10), we describe INI as Algorithm 1. Algorithm 1
Inexact Noda Iteration (INI)1. Given b l , x > with k x k =
1, 0 ≤ g < tol > for k = , , , . . .
3. Solve ( b l k A − I ) y k + = A x k inexactly such that the inner tolerance x k satisfies condition(10)
4. Normalize the vector x k + = y k + / k y k + k .5. Compute b l k + = max (cid:16) A − x k + x k + (cid:17) .6. until convergence: Resi = k A x k + − b l − k x k + k < tol .Using the relation (7), step 5 in Algorithm 1 can be rewritten as b l k + = b l k − min (cid:18) x k + A − f k y k + (cid:19) . Unfortunately, A − is not explicitly available; in other words, we need to compute “ A − f k ”exactly for the required approximate eigenvalue b l k + . Hence, in the next section, we propose a new strategyto estimate the approximate eigenvalues without increasing the computational cost. This strategy ispractical and preserves the strictly decreasing property of the approximate eigenvalue sequence. Copyright c (cid:13)
Numer. Linear Algebra Appl. (2010)
Prepared using nlaauth.clsnlaauth.cls
Prepared using nlaauth.clsnlaauth.cls
DOI: 10.1002/nlaNEXACT NODA ITERATION
53. THE RELAXATION STRATEGY FOR INI AND SOME BASIC PROPERTIESIn order to ensure that INI is correctly implemented, we now propose two main relaxation steps todefine Algorithm 2: • The residual norm satisfies x k = k f k k ≤ g k sep ( , A ) min ( x k ) , (11)where 0 ≤ g k ≤ g < g . • The update of the approximate eigenvalue satisfies l k + = l k − ( − g k ) min (cid:18) x k y k + (cid:19) . (12) Algorithm 2
Inexact Noda Iteration for monotone matrices (INI)
1. Given l , x > with k x k =
1, 0 ≤ g < tol > for k = , , , . . .
3. Solve ( l k A − I ) y k + = A x k inexactly such that the inner tolerance x k satisfies condition(11)4. Normalize the vector x k + = y k + / k y k + k .5. Compute l k + that satisfies condition (12).6. until convergence: k A x k + − l − k x k + k < tol .In step 3 of Algorithm 2, it leads to two equivalent inexact relation satisfying ( l k I − A − ) y k + = x k + A − f k (13) ( l k A − I ) y k + = A x k + f k , (14)We remark that l k + in (12) is no longer equal to max (cid:16) A − x k + x k + (cid:17) , and therefore that l k + cannot beclearly demonstrated to be greater than its lower bound r ( A − ) . The following lemma ensures that r ( A − ) is still the lower bound of l k . Lemma 3
Let A be an irreducible monotone matrix. For the unit length x k = x > and the relaxation factor g k ∈ [ , ) , if l k > r ( A − ) , f k in (
14) satisfies condition (11) and the approximate eigenvalue satisfies(12), then the new approximation x k + > and the sequence n l k o is monotonically decreasing andbounded below by r ( A − ) , i.e., l k > l k + ≥ r ( A − ) . Proof
From (12) and g k ∈ [ , ) , it is easy to know that n l k o is monotonically decreasing, i.e., l k > l k + . From (11), (cid:13)(cid:13) A − f k (cid:13)(cid:13) ≤ (cid:13)(cid:13) A − (cid:13)(cid:13) k f k k ≤ g k min ( x k ) , (15)which implies (cid:12)(cid:12) A − f k (cid:12)(cid:12) ≤ g k x k , then x k + A − f k > . Consequently, l k I − A − is a nonsingular M -matrix, and the vector y k + satisfies y k + = ( l k I − A − ) − (cid:0) x k + A − f k (cid:1) > . This implies x k + = y k + / k y k + k > . Copyright c (cid:13)
Numer. Linear Algebra Appl. (2010)
Prepared using nlaauth.clsnlaauth.cls
Prepared using nlaauth.clsnlaauth.cls
DOI: 10.1002/nla
C.-S. LIU
We now prove l k is bounded below by r ( A − ) . From (15), we have ( − g k ) x k ≤ x k + A − f k ≤ ( + g k ) x k , (16)and ( − g k ) x k y k + ≤ x k + A − f k y k + ≤ ( + g k ) x k y k + . Hence, we obtain ( − g k ) min (cid:18) x k y k + (cid:19) ≤ min (cid:18) x k + A − f k y k + (cid:19) ≤ ( + g k ) min (cid:18) x k y k + (cid:19) , then (cid:12)(cid:12)(cid:12) (cid:12) min (cid:18) x k + A − f k y k + (cid:19) − min (cid:18) x k y k + (cid:19)(cid:12)(cid:12)(cid:12) (cid:12) ≤ g k min (cid:18) x k y k + (cid:19) . (17)From (13), it follows that r ( A − ) ≤ b l k + = max (cid:18) A − x k + x k + (cid:19) = l k − min (cid:18) x k + A − f k y k + (cid:19) . (18)Combine (12), (18) and (17), then l k + = l k − ( − g k ) min (cid:18) x k y k + (cid:19) = b l k + + min (cid:18) x k + A − f k y k + (cid:19) − ( − g k ) min (cid:18) x k y k + (cid:19) (19) ≥ b l k + + ( − g k − ( − g k )) min (cid:18) x k y k + (cid:19) > r ( A − ) . By induction, n l k o is bounded below by r ( A − ) , i.e., l k > r ( A − ) for all k . From Lemma 3, since n l k o is a monotonically decreasing and bounded sequence, we musthave lim k → ¥ l k = a ≥ r ( A − ) , where a = r ( A − ) or a > r ( A − ) . We next investigate the case a > r ( A − ) , and present some basic results; this plays an important role later in proving theconvergence of INI. Lemma 4
For Algorithm 2, if l k is converge to a > r ( A − ) , then (i) k y k k is bounded;(ii) lim k → ¥ min ( x k ) = ∠ ( x k , x ) ≥ m > m >
0, where ∠ ( x k , x ) the acute angle of x k and x . Proof (i). From ( k y k + k = k (cid:16) l k I − A − (cid:17) − (cid:0) x k + A − f k (cid:1) k ≤ ( + g k ) k ( l k I − A − ) − k = ( l k , A − ) − ≤ ( a , A − ) − < ¥ . (20) Copyright c (cid:13)
Numer. Linear Algebra Appl. (2010)
Prepared using nlaauth.clsnlaauth.cls
Prepared using nlaauth.clsnlaauth.cls
DOI: 10.1002/nlaNEXACT NODA ITERATION k → ¥ min (cid:18) x k y k + (cid:19) = lim k → ¥ l k − l k + − g k ! = . (21)On the other hand, from (20) and (21) we havemin (cid:18) x k y k + (cid:19) ≥ min ( x k ) max ( y k + ) ≥ min ( x k ) sep ( a , A − ) > . Thus, it is holds that lim k → ¥ min ( x k ) = . (22)(iii) Suppose there is a subsequence { sin ∠ ( x k j , x ) } which converges to zero, thenlim j → ¥ b l k j = lim j → ¥ max A − x k j x k j ! = max lim j → ¥ A − x k j x k j ! = r ( A − ) . By the definition of b l k , from (19) and (17), we have (cid:12)(cid:12)(cid:12)b l k − l k (cid:12)(cid:12)(cid:12) = (cid:12)(cid:12)(cid:12)(cid:12) min (cid:18) x k − + A − f k − y k (cid:19) − ( − g k ) min (cid:18) x k − y k (cid:19)(cid:12)(cid:12)(cid:12)(cid:12) ≤ (cid:12)(cid:12)(cid:12)(cid:12) ( + g k − ( − g k )) min (cid:18) x k − y k (cid:19)(cid:12)(cid:12)(cid:12)(cid:12) for any k . (23)From (21) and (23), lim j → ¥ l k j = lim j → ¥ (cid:16) l k j − b l k j + b l k j (cid:17) = r ( A − ) . This is a contradiction.
Lemma 5 ([15])Let x > be the unit length eigenvector of A associated with s min ( A ) . For any vector z > with k z k =
1, it holds that cos ∠ ( z , x ) > min ( x ) andinf k z k = , z > cos ∠ ( z , x ) = min ( x ) . (24)4. CONVERGENCE ANALYSIS FOR INI In Sections 4.1–4.2, we will prove the global convergence and the convergence rate of INI.Furthermore, we will derive the explicit linear convergence factor and the superlinear convergenceorder with different g k . For an irreducible non-negative matrix A − , recall that the largest eigenvalue r ( A − ) of A − issimple. Let x be the unit length positive eigenvector corresponding to r ( A − ) . Then for anyorthogonal matrix (cid:2) x V (cid:3) it holds (cf. [10]) that (cid:20) x T V T (cid:21) A − (cid:2) x V (cid:3) = (cid:20) r ( A − ) c T L (cid:21) (25) with L = V T A − V whose eigenvalues constitute the other eigenvalues of A − . Therefore, we nowdefine e k = l k − r ( A − ) , A k = l k I − A − . (26) Copyright c (cid:13)
Numer. Linear Algebra Appl. (2010)
Prepared using nlaauth.clsnlaauth.cls
Prepared using nlaauth.clsnlaauth.cls
DOI: 10.1002/nla
C.-S. LIU
Similar to (25) we also have the spectral decomposition (cid:20) x T V T (cid:21) A k (cid:2) x V (cid:3) = (cid:20) e k L k (cid:21) , (27)where L k = l k I − L . For l k = r ( A − ) , it is easy to verify that (cid:20) x T V T (cid:21) A − k (cid:2) x V (cid:3) = (cid:20) e − k b Tk L − k (cid:21) with b Tk = − c T L − k e k , (28)from which we get A − k V = xb Tk + V L − k = − x c T L − k e k + V L − k , (29) A − k = e − k xx T − e − k xc T L − k V T + V L − k V T , (30)and x T A − k = e − k x T − e − k c T L − k V T . Let { x k } be generated by Algorithm 2. We decompose x k into the orthogonal direct sum by x k = x cos ( j k ) + p k sin ( j k ) , p k ∈ span ( V ) ⊥ x (31)with k p k k = j k = ∠ ( x k , x ) being the acute angle between x k and x . So by definition, we havecos j k = x T x k and sin j k = k V T x k k . Evidently, x k → x if and only if tan j k →
0, i.e., sin j k → x k = k f k k ≤ g k sep ( , A ) min ( x k ) in INI, it holds that (cid:12)(cid:12) A − f k (cid:12)(cid:12) ≤ g k x k . Therefore, we have ( − g k ) x k ≤ x k + A − f k ≤ ( + g k ) x k . (32)As A − k ≥
0, it follows from the above relation that ( − g k ) A − k x k ≤ y k + ≤ ( + g k ) A − k x k . (33)Using the above relation, we obtaintan j k + = sin j k + cos j k + = k V T x k + k x T x k + = k V T y k + k x T y k + = k V T A − k ( x k + A − f k ) k x T A − k ( x k + A − f k ) = k L − k V T (cid:0) x k + A − f k (cid:1) k (cid:0) e − k x T − e − k c T L − k V T (cid:1) ( x k + A − f k )= k L − k V T (cid:0) x k + A − f k (cid:1) k e − k x T x k − e − k c T L − k V T x k + e − k x T A − f k − e − k c T L − k V T A − f k ≤ k L − k k e k (cid:0) sin j k + k A − f k k (cid:1) cos j k − c T L − k V T x k − k A − f k k − k c kk L − k kk A − f k k . (34)Note that if we solve the inner linear system exactly, i.e., f k =
0, we recover NI and gettan j k + ≤ k L − k k e k − c T L − k V T x k / cos j k tan j k : = b k tan j k . (35)Since L. Elsner [9] proved the quadratic convergence of the proposed Noda iteration, for k largeenough we must have b k = O ( tan j k ) → . Copyright c (cid:13)
Numer. Linear Algebra Appl. (2010)
Prepared using nlaauth.clsnlaauth.cls
Prepared using nlaauth.clsnlaauth.cls
DOI: 10.1002/nlaNEXACT NODA ITERATION b k < b < b <
1. Therefore, we havetan j k + < b tan j k for k ≥ N with N large enough. Theorem 3 (Main Theorem)Let A be an irreducible monotone matrix. If the sequence { l k } is generated by INI with therelaxation strategies (11) and (12), then { l k } is monotonically decreasing and lim k → ¥ l − k = s min ( A ) . Proof
From Lemma 3, the sequence n l k o is bounded and monotonically decreasing, and we must haveeither lim k → ¥ l k = r ( A − ) or lim k → ¥ l k = a > r ( A − ) . Next we prove by contradiction that, forINI, lim k → ¥ l k = r ( A − ) must hold. Suppose that lim k → ¥ l k = a > r ( A − ) . It follows (iii) of Lemma 4 show that1cos j k ≤ j k sin j k m = tan j k m . (37)From (ii) of Lemma 4, we have lim k → ¥ min ( x k ) = . This implies the inner tolerance k f k k →
0, i.e., k A − f k k is suitably small. In addition, by Lemma 5,it holds that cos j k is uniformly bounded below by min ( x ) , therefore, ( + k c kk L − k ) k A − f k k / cos j k < − c T L − k V T x k / cos j k (38)for k large enough.Using (34), (37) and (38), we obtaintan j k + ≤ k L − k k e k (cid:0) tan j k + k A − f k k / cos j k (cid:1) − c T L − k V T x k / cos j k − ( + k c kk L − k k ) k A − f k k / cos j k ≤ k L − k k e k ( tan j k + g k min ( x k ) tan j k / m ) − c T L − k V T x k / cos j k − ( + k c kk L − k k ) g k min ( x k ) tan j k / m ≤ k L − k k e k ( + g k min ( x k ) / m ) (cid:0) − c T L − k V T x k / cos j k (cid:1) − ( + k c kk L − k k ) g k min ( x k ) tan j k / m tan j k . Define b ′ k = k L − k k e k ( + g k min ( x k ) / m ) (cid:0) − c T L − k V T x k / cos j k (cid:1) − ( + k c kk L − k k ) g k min ( x k ) tan j k / m . Note that b ′ k is a continuous function with respect to min ( x k ) for 0 < g k <
1. Then it holds that b ′ k → b k defined by (35) as min ( x k ) →
0. Therefore, from (36), for k large enough we can choose asufficiently small d such that b ′ k ≤ ( + d ) b k < b < ( x k ) sufficiently small. As a result, we havetan j k + ≤ b tan j k for k ≥ N with N large enough and min ( x k ) sufficiently small. This means that tan j k →
0, i.e.,sin j k → . From (iii) of Lemma 4, sin j k is uniformly bounded below by a positive constant. So sin j k → and sin j k ≥ m , a contradiction. Therefore the initial assumption “lim k → ¥ l k = a > r ( A − ) ”must be false. Copyright c (cid:13)
Numer. Linear Algebra Appl. (2010)
Prepared using nlaauth.clsnlaauth.cls
Prepared using nlaauth.clsnlaauth.cls
DOI: 10.1002/nla C.-S. LIU
Theorem 3 has proved the global convergence of INI, but the results are only qualitative and donot tell us anything about how fast the INI method converges. In this subsection, we will showthe convergence rate of INI with different relaxation factors g k . More precisely, we prove that INIconverges at least linearly with an asymptotic convergence factor bounded by g + g for 0 ≤ g k ≤ g < g k = l k − − l k l k − , respectively.From (12) we have e k + = e k (cid:18) − ( − g k ) min (cid:18) x k e k y k + (cid:19)(cid:19) = : e k r k . (39)Since l k − l k + < l k − r ( A − ) , from (39) and (12), we have r k = − ( − g k ) min (cid:18) x k e k y k + (cid:19) = − l k − l k + l k − r ( A − ) ! < . Theorem 4
For INI, we have r k ≤ g k + g k < . Moreover, if g k ≤ g < , then lim k → ¥ r k ≤ g + g <
1, i.e., the convergenceof INI is at least globally linear. If lim k → ¥ g k = , then lim k → ¥ r k = , that is, the convergence of INI isglobally superlinear. Proof
From (32) and (33), we havemin (cid:18) x k e k y k + (cid:19) ≥ min x k ( + g k ) e k A − k x k ! = + g k min x k e k A − k x k ! . From (30), we get e k A − k x k = xx T x k − xc T L − k V T x k + e k V L − k V T x k . (40)From Theorem 3, we know that lim k → ¥ x k = x and lim k → ¥ l k = r ( A − ) , from which it follows that e k → L − k → (cid:0) r ( A − ) I − L (cid:1) − . On the other hand, since L − k → ( r ( A − ) I − L ) − andlim k → ¥ V T x k = V T x =
0, from (40) we get lim k → ¥ e k A − k x k = x . Consequently, we obtainlim k → ¥ min (cid:18) x k e k y k + (cid:19) ≥ + g k min lim k → ¥ x k e k B − k x k ! = + g k min (cid:16) xx (cid:17) = + g k > , leading to r k ≤ − − g k + g k = g k + g k < . (41) Copyright c (cid:13)
Numer. Linear Algebra Appl. (2010)
Prepared using nlaauth.clsnlaauth.cls
Prepared using nlaauth.clsnlaauth.cls
DOI: 10.1002/nlaNEXACT NODA ITERATION g k is small then INI must ultimately converge quickly. AlthoughTheorem 4 has established the superlinear convergence of INI, it does not reveal the convergenceorder. Our next concern is to derive the precise convergence order of INI. This is more informativeand instructive because it lets us understand how fast INI converges. Theorem 5
If the inner tolerance x k in INI satisfies condition (11) with the relaxation factors g k = l k − − l k l k − , (42)then INI converges quadratically (asymptotically) in the form of e k ≤ e k − (43)for k large enough, where the relative error e k + = e k / r ( A − ) . Proof
Since l k − > l k > r ( A − ) , we have g k = l k − − l k l k − ≤ l k − − r ( A − ) r ( A − ) = e k − r ( A − ) . (44)From (39), (41) and (44), we have e k = e k − r k − ≤ e k − g k + g k = e k − + g k ≤ e k − + r ( A − ) e k − = e k − e k − + r ( A − ) ≤ r ( A − ) e k − . Dividing both sides of the above inequality by r ( A − ) , we get e k = e k r ( A − ) ≤ r ( A − ) e k − = e k − .
5. THE MODIFIED INEXACT NODA ITERATIONIn this section, we propose a modified Noda iteration (MNI) for a non-negative matrix, and showthat MNI and NI are equivalent. Thus, by combining INI (Algorithm 2) with MNI we can proposea modified inexact Noda iteration for a monotone matrix
When b l k I − B tends to a singular matrix, the Noda iteration requires us to solve a possibly ill-conditioned linear system (3). Hence, we propose a rank one update technique for the ill-conditionedlinear system (3), i.e., (cid:18) B − b l k I − x k − x Tk (cid:19) (cid:18) D y k d k (cid:19) = (cid:20) ( b l k I − B ) x k (cid:21) , (45) Copyright c (cid:13)
Numer. Linear Algebra Appl. (2010)
Prepared using nlaauth.cls
DOI: 10.1002/nla C.-S. LIU where D y k = x k + − x k . Let r k = ( b l k I − B ) x k . In general, the linear system (45) is a well-conditioned linear system, unless B has the Jordan form corresponding to the largest eigenvalue,which contradicts the Perron–Frobenius theorem.From (45), 0 = (cid:16) B − b l k I (cid:17) ( x k + − x k ) − d k x k − r k = (cid:16) B − b l k I (cid:17) x k + − (cid:16) B − b l k I (cid:17) x k − d k x k − r k = (cid:16) B − b l k I (cid:17) x k + − d k x k . Hence, we have the following linear system (cid:16)b l k I − B (cid:17) (cid:18) x k + − d k (cid:19) = x k , or hb l k I − B i y k + = x k , with y k = − d k x k + . Thus, from (3) and (45), we have the new iterative vector x k + = y k + k y k + k = x k + D y k k x k + D y k k . (46)This means the Noda iteration and the modified Noda iteration are mathematically equivalent. Basedon (45) and (46), we state our algorithm as follows. Algorithm 3
Modified Noda Iteration (MNI)1. Given b l , x > with k x k = tol > for k = , , , . . . if k B x k + − b l k x k + k > √ tol
4. Solve ( b l k I − B ) y k + = x k .5. Normalize the vector x k + = y k + / k y k + k .6. else if
7. Solve (cid:18) B − b l k − x k − x Tk (cid:19) (cid:18) D y k d k (cid:19) = (cid:20) b l k x k − B x k (cid:21) .
8. Normalize the vector x k + = ( x k + D y k ) / k x k + D y k k .9. end
10. Compute b l k + = max (cid:16) B x k + x k + (cid:17) . until convergence: k B x k + − b l k x k + k < tol .Note that the sequence { b l k I − B } tends to a singular matrix, meaning that { b l k } converges to aneigenvalue of B , and (3) becomes an ill-conditioned linear system. Based on practical experiments,we propose taking k B x k + − b l k x k + k ≤ √ tol and switching from (3) to (45). For a monotone matrix A , we replaced B by A − in (45). The linear system (45) can be rewritten as (cid:18) I − b l k A − A x k − x Tk (cid:19) (cid:18) D y k d k (cid:19) = (cid:20) b l k A x k − x k (cid:21) . (47)Based on MNI, by combining INI (Algorithm
2) with equation ( ), we can propose a modifiedinexact Noda iteration for a monotone matrix, which is described as Algorithm 4. Copyright c (cid:13)
Numer. Linear Algebra Appl. (2010)
Prepared using nlaauth.clsnlaauth.cls
Prepared using nlaauth.clsnlaauth.cls
DOI: 10.1002/nlaNEXACT NODA ITERATION Algorithm 4
Modified Inexact Noda Iteration (MINI)1. Given b l , x > with k x k = tol > for k = , , , . . . if = k A x k + − l − k x k + k > √ tol
4. Run INI for monotone matrix A (Algorithm 2).5. else if
6. Solve (cid:18) I − l k A − A x k − x Tk (cid:19) (cid:18) D y k d k (cid:19) = (cid:20) l k A x k − x k (cid:21) exactly.7. Normalize the vector x k + = ( x k + D y k ) / k x k + D y k k .8. Compute l k + that satisfies condition (12).9. end until convergence: k A x k + − l − k x k + k < tol .6. NUMERICAL EXPERIMENTS In this section we present numerical experiments to support our theoretical results for INI, and toillustrate the effectiveness of the proposed MINI algorithms. All numerical tests were performed onan Intel (R) Core (TM) i7 CPU 4770@ 3 . a with themachine precision e = . × − under the Microsoft Windows 7 64-bit. I outer denotes the number of outer iterations to achieve the convergence, and I inner denotes thetotal number of inner iterations, which measures the overall efficiency of MNI and MINI. In viewof the above, we have the average number I ave = I inner / I outer at each outer iteration for our testalgorithms. In the tables, “Positivity”illustrates whether the converged Perron vector preserves thestrict positivity property. If “No”, then the percentage in the brace indicates the proportion thatthe converged Perron vector has the positive components. We also report the CPU time of eachalgorithm, which measures the overall efficiency too. We present an example to illustrate the numerical behavior of NI, INI 1 and INI 2 for monotonematrices. The approximate solution y k + of (14) satisfies ( l k A − I ) y k + = A x k + f k by requiring the following inner tolerances: • for NI: k f k k ≤ − ; • for INI 1: k f k k ≤ g k sep ( , A ) min ( x k ) with some 0 < g k < • for INI 2: k f k k ≤ l k − − l k l k − sep ( , A ) min ( x k ) for k ≥ l > r ( A − ) .We use the minimal residual method to solve the inner linear systems. For the implementations, weuse the standard Matlab function minres . The outer iteration starts with the normalized vector of ( , . . ., ) T , and the stopping criterion for outer iterations is k A x k − l − k x k k ( k A k k A k ¥ ) / ≤ − , where k · k and k · k ¥ are the one norm and the infinity norm of a matrix, respectively.Condition (11) ensures that the eigenvector in Lemma 3 does indeed preserve the strict positivityproperty. However, the formula in (
11) is not applicable in practice, because it uses sep ( , A ) , whichis unknown at the time it needs to be computed . Therefore, for practical implementations, wesuggest a relaxation strategy to replace sep ( , A ) by l − k . The quantity l − k is related to the lower Copyright c (cid:13)
Numer. Linear Algebra Appl. (2010)
Prepared using nlaauth.clsnlaauth.cls
Prepared using nlaauth.clsnlaauth.cls
DOI: 10.1002/nla C.-S. LIU −15 −10 −5 outer iterations r e s i dua l no r m s INI−1 with g k = 0.8INI−1 with g k = 0.5INI−2NI Figure 1. The outer residual norms versus outer iterations in Examples 1.Table I. The total outer and inner iterations in Example 1
Method I outer I inner I ave CPU time PositivityINI 1 with g = . g = . A , i.e., s min ( A ) ≥ l − k . For all examples, the stopping criterionfor the inner iteration is set at k f k k ≤ max { g k min ( x k ) / l k , − } for INI 1and k f k k ≤ max { l k − − l k l k − l k min ( x k ) , − } for INI 2. Example 1
We consider the finite-element discretization of the boundary value problem in [3, Example 4.2.4] − u xx − u yy = g ( x , y ) in W = [ , a ] × [ , b ] , a , b > , u = f ( x , y ) on ¶ W , using piecewise quadratic basis functions on the uniform mesh of p × m isosceles right-angled triangles. This is a matrix of order n = ( p − )( m − ) = ,
041 with p =
400 and m = g k = . . g k = . , . I ave for INI 1 isonly half I ave of INI 2. There are two reasons for this. First, since the approximate eigenvalues areobtained from the relation (12), then the parameter g k will lead to a difference in the convergencerates, as is seen from Figure
1. Second, from (11), INI 2 solves the inner linear systems more and more accurately as k increases . In contrast, the inner tolerance used by INI 1 is fixed except for thefactor min ( x k ) , which also makes the average number of the iterations of INI 1 only about half ofthose for INI 2. Copyright c (cid:13)
Numer. Linear Algebra Appl. (2010)
Prepared using nlaauth.clsnlaauth.cls
Prepared using nlaauth.clsnlaauth.cls
DOI: 10.1002/nlaNEXACT NODA ITERATION In the above section, INI 2 was considerably better than NI and INI 1 for overall efficiency.Therefore, in this subsection, we use MINI (INI 2 combined with MNI) to find the smallest singularvalue and the associated eigenvector of an M -matrix, and confirm the effectiveness of MINI andthe theory we presented in Sections 3 and 4. For MINI, the stopping criteria for inner and outeriterations are the same as those for monotone matrices. In the meantime , we compare MINIwith the algorithms JDQR [27] and JDSVD [11] and the Matlab function svds ; none of theseare positivity preserving for approximate eigenvectors. We show that the MINI algorithm alwaysreliably computes positive eigenvectors, while the other algorithms generally fail to do so.Since JDQR and JDSVD use the absolute residual norms to decide the convergence, then weset the stopping criteria “TOL = − ( k A k k A k ¥ ) / ” for outer iterations, and then we will getthe same stopping criteria as used for MINI. We set the parameters “sigma=SM” for JDQR,“opts.target=0” for JDSVD, and the inner solver “OPTIONS.LSolver=minres”. All the other optionsuse defaults. We do not use any preconditioning for inner linear systems. For svds , we set thestopping criteria “OPTS.tol = − ( k A k k A k ¥ ) / , and take the maximum and minimum subspacedimensions as 20 and 2 at each restart, respectively. Suppose that we want to compute the smallest singular value, and the corresponding singularvector, of a real n × n M-matrix M . This partial SVD can be computed by using equivalent eigenvaluedecomposition, that is, the augmented matrix A = (cid:20) MM T (cid:21) . Obviously, such a matrix A is no longer an M -matrix but will indeed be monotone. Example 2
We consider a symmetric M -matrix of the form M = s I − B , where B is the non-negative matrix rgg n 2 19 s0 from the DIMACS10 test set [8]. This matrix is a random geometric graph with2 vertices. Each vertex is a random point in the unit square and edges connect vertices whoseEuclidean distance is below 0 .
55 (log ( n ) / n ). This threshold is chosen in order to ensure that thegraph is almost connected. This matrix is a binary matrix with n = = ,
288 and 6 , , svds compute the desired eigenvalue, but the converged eigenvectors are not positive. More precisely,Table II indicates that for these algorithms roughly 50% of the components of each convergedeigenvector are negative.As far as overall efficiency is concerned, MINI is the most efficient in terms of I inner , I outer andthe CPU time. JDQR and svds require at least five times the CPU time of MINI; they are also more expensive than JDSVD in terms of the CPU time. Table II. The total outer and inner iterations in Example 2
Method I outer I inner I ave CPU time PositivityMINI 6 331 55 30 YesJDQR 25 4068 162 243 No (52%)JDSVD 34 1432 42 58 No (51%) svds
140 —– 140 144 No (57%)7. CONCLUSIONS
We have proposed an inexact Noda iteration method for computing the smallest eigenpair of alarge irreducible monotone matrix, and have considered the convergence of the modified inexact
Copyright c (cid:13)
Numer. Linear Algebra Appl. (2010)
Prepared using nlaauth.clsnlaauth.cls
Prepared using nlaauth.clsnlaauth.cls
DOI: 10.1002/nla C.-S. LIU
Noda iteration with two relaxation factors. We have proved that the convergence of INI is globallylinear and superlinear, with the asymptotic convergence factor bounded by g k + g k . More precisely,the modified inexact Noda iteration with inner tolerance x k = k f k k ≤ g k sep ( , A ) min ( x k ) convergesat least linearly if the relaxation factors meet the condition g k ≤ g <
1, and superlinearly if therelaxation factors meet the condition g k = l k − − l k l k − , respectively. The results for INI clearly showhow the accuracy of the inner iterations affects the convergence of the outer iterations.In the experiments, we also compared MINI with Jacobi–Davidson type methods (JDQR,JDSVD) and the implicitly restarted Arnoldi method ( svds ). The contribution of this paper istwofold. First, MINI always preserves the positivity of approximate eigenvectors, while the otherthree methods often fail to do so. Second, the proposed MINI algorithms have been shown to bepractical and effective for large monotone matrix eigenvalue problems and M -matrix singular valueproblems. REFERENCES1. A. S. A
LFA , J. X
UE AND
Q. Y E , Accurate computation of the smallest eigenvalue of a diagonally dominantM-matrix , Math. Comp., 71, pp. 217–236 (2002).2. J. A
BATE , L. C
HOUDHURY AND
W. W
HITT , Asymptotics for steady-state tail probabilities in structured Markovqueueing models , Commun. Statist-Stochastic Models, 1, pp. 99–143 (1994)3. O. A
XELSSON , L. K
OLOTILINA , Monotonicity and discretization error estimates , SIAM J.Numer.Anal., 27 (1990),1591-1611.4. O.A
XELSSON , S.V.G
OLOLOBOV , Monotonicity and discretization error estimates for convection-diffusionproblems,
Technical Report, Department of mathematic university of Nijmegen (2003).5. A. B
ERMAN AND
R. J. P
LEMMONS , Non-negative Matrices in the Mathematical Sciences , Vol. 9 of Classics inApplied Mathematics, SIAM, Philadelphia, PA (1994)6. J. B
ERNS -M ¨
ULLER , I. G. G
RAHAM AND
A. S
PENCE , Inexact inverse iteration for symmetric matrices , LinearAlgebra Appl., 416, pp. 389–413 (2006)7. L. C
OLLATZ , Numerical Treatment of Differential Equations , 3rd ed., Springer, berlin, (1960)8.
DIMACS10 test set and the University of Florida Sparse Matrix Collection . Available at .9. L. E
LSNER , Inverse iteration for calculating the spectral radius of a non-negative irreducible matrix , LinearAlgebra Appl., 15, pp. 235–242 (1976)10. G. H. G
OLUB AND
C. F. V AN L OAN , Matrix Computations , The John Hopkins University Press, Baltimore,London, 4th ed. (2012)11. M.E. H
OCHSTENBACH , A Jacobi-Davidson type SVD method , SIAM J. Sci. Comp. 23(2), pp. 606–628 (2001)12. R. A. H
ORN AND
C. R. J
OHNSON , Matrix Analysis , The Cambridge University Press, Cambridge, UK (1985)13. Z. J IA , On convergence of the inexact Rayleigh quotient iteration with MINRES , J. Comput. Appl. Math., 236,pp. 4276–4295 (2012)14. ,
On convergence of the inexact Rayleigh quotient iteration with the Lanczos method used for solving linearsystems , technical report, September 2011. http://arxiv.org/pdf/0906.2239v4.pdf .15. Z. J IA , W.-W. L IN , AND
C.-S. L IU , A positivity preserving inexact Noda iteration for computing the smallesteigenpair of a large irreducible M-matrix , Numer. Math., accepted, 2014.16. C. L EE , Residual Arnoldi method: theory, package and experiments , Ph D thesis, TR-4515, Department ofComputer Science, University of Maryland at College Park, (2007)17. Y.-L. L AI , K.-Y. L IN , AND
W.-W. L IN , An inexact inverse iteration for large sparse eigenvalue problems , Numer.Linear Algebra Appl., 4 (1997), pp. 425–437.18. T. N
ODA , Note on the computation of the maximal eigenvalue of a non-negative irreducible matrix , Numer. Math.,17 (1971), pp. 382–386.19. D. D. O
LESKY AND P. VAN DEN D RIESSCHE , Monotone positive stable matrices , Linear Algebra Appl., Volume48, pp. 381–401 (1983).20. A. M. O
STROWSKI , On the convergence of the Rayleigh quotient iteration for the computation of the characteristicroots and vectors. V. (Usual Rayleigh quotient for non-Hermitian matrices and linear elementary divisors) , Arch.Rational Mech. Anal., 3 (1959), pp. 472–481.21. B. N. P
ARLETT , The Symmetric Eigenvalue Problem , Society for Industrial and Applied Mathematics (SIAM),Philadelphia, PA, 1998.22. Y. S
AAD , Numerical Methods for Large Eigenvalue Problems , Manchester University Press, Manchester, UK,1992.23. J. S
CHROEDER , M-Matrices and generalizations using an operator theory approach , SIAM Rev. 20, 213-244(1978).24. P. S
HIVAKUMAR , J. W
ILLIAMS , Q. Y
E AND
C. M
ARINOV , On two-sided bounds related to weakly diagonallydominant M-matrices with applications to digital circuit dynamics , SIAM J. Matrix Anal. Appl., 17, pp. 298–312(1996)25. K.C. S
IVAKUMAR , A new characterization of nonnegativity of Moore–Penrose inverses of Gram operators ,Positivity 13 (1) (2009), 277–286.Copyright c (cid:13)
Numer. Linear Algebra Appl. (2010)
Prepared using nlaauth.cls
DOI: 10.1002/nlaNEXACT NODA ITERATION
26. G. L. G. S
LEIJPEN AND
H. A.
VAN DER V ORST , A Jacobi–Davidson iteration method for linear eigenvalueproblems , SIAM J. Matrix Anal. Appl., 17 (1996), pp. 401–425.27. G. L. G. S
LEIJPEN , JDQR and JDQZ codes . Available at
28. G. W. S
TEWART , Matrix Algorithms. Vol. II , Society for Industrial and Applied Mathematics (SIAM), Philadelphia,PA, 2001.29. J. X UE , Computing the smallest eigenvalue of an M-matrix , SIAM J. Matrix Anal. Appl., 17 (1996), pp. 748–762.Copyright c (cid:13)
Numer. Linear Algebra Appl. (2010)
Prepared using nlaauth.clsnlaauth.cls