Self-Corrective Iterations (SCI) for generalized diagonally dominant matrices

Accepted Manuscript Self-Corrective Iterations (SCI) for generalized diagonally dominant matrices Jinrui Guan, Linzhang Lu, Ren-Cang Li, Rongxia Shao PII: DOI: Reference:

S0377-0427(16)30068-1 http://dx.doi.org/10.1016/j.cam.2016.02.021 CAM 10518

To appear in:

Journal of Computational and Applied Mathematics

Received date: 6 June 2015 Revised date: 22 January 2016 Please cite this article as: J. Guan, L. Lu, R.-C. Li, R. Shao, Self-Corrective Iterations (SCI) for generalized diagonally dominant matrices, Journal of Computational and Applied Mathematics (2016), http://dx.doi.org/10.1016/j.cam.2016.02.021 This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

Self-Corrective Iterations (SCI) For Generalized Diagonally Dominant Matrices Jinrui Guana , Linzhang Lub,a,1 , Ren-Cang Lic,2,∗, Rongxia Shaoa a School of Mathematical Science, Xiamen University, Fujian Province, P. R. China of Mathematics and Computer Science, Guizhou Normal University, Guizhou Province, P. R. China c Department of Mathematics, University of Texas at Arlington, P.O. Box 19408, Arlington, TX 76019, USA b School

Abstract A suggestive indicator is proposed for predicting whether a given (complex or real) square matrix A is or isn’t a generalized diagonally dominant matrix (GDDM) by which we mean if A can be brought into a strictly diagonally dominant matrix by post-multiplying some diagonal matrix D. Based on the indicator, three self-corrective iterative algorithms (SCI) are presented for determining if an irreducible A is or is not a GDDM and at the same time delivering the matrix D in case when A is a GDDM. The basic idea is to push A towards being (strictly) diagonally dominant when the indicator suggests that A is likely a GDDM or towards being off-diagonally dominant otherwise. Among the three algorithms, each has their own feature: one takes the fewest number of iterations but the most amount of work per iteration, one takes the most number of iterations but the least amount of work per iteration, and the third one strikes a balance between the two extremes. It is shown that the algorithms will terminate in finite many steps under the assumption that A has no zero entries and the comparison matrix of A is nonsingular. Comparing with existing methods, new algorithms are more efficient, as demonstrated on known difficult examples in the literature, newly designed random matrices, as well as matrices from the University of Florida Sparse Matrix Collection. Keywords: Generalized diagonally dominant matrix, GDDM, M-matrix, H-matrix, self-corrective iteration, SCI 2010 MSC: 15B99, 65F10, 65F35

1. Introduction A (complex or real) square matrix A is called a generalized diagonally dominant matrix (GDDM), also called a nonsingular H-matrices, if there is a diagonal matrix D such that AD is strictly diagonally dominant. GDDMs play important roles in numerical analysis, matrix theory, control ∗ Corresponding

author Email addresses: [email protected] (Jinrui Guan), [email protected] (Linzhang Lu), [email protected] (Ren-Cang Li), [email protected] (Rongxia Shao) 1 Supported in part by NNSFC grant 10261012 and 11428104. 2 Supported in part by NSF grants DMS-1115834, DMS-1317330, CCF-1527104, and NNSFC grant 11428104.

Preprint submitted to J. of Comput. and Applied Math.

February 18, 2016

systems, and mathematical economics, to name a few (see, e.g., [1, 2, 3]). In certain applications, it’s critically important to know whether a given matrix is a GDDM or not. In general this is not an easy task. In recent years, there are quite few studies on this subject, and several methods some of which are quite efficient have been proposed for determining if a given A is or isn’t a GDDM and possibly a diagonal matrix D in case when it is GDDM (see, e.g., [4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15]). These methods can be divided into two categories: direct methods and iterative ones. It was argued in Alanelli and Hadjidimos [4] that direct methods have high computational complexities, while iterative ones are usually much cheaper and more efficient, especially for large scale sparse matrices. Among studies on iterative methods, we mention these papers [4, 5, 10, 11, 12, 13, 14, 15]. Here in this paper, our interest is also in iterative methods. Throughout the rest, A ≡ [aij ] is always an n × n (complex or real) matrix with aij being its (i, j)th entry. Let N = {1, 2, · · · , n}, and define ri (A) =

X

N∋j6=i

|aij |, ti (A) =

ri (A) |aii |

for i ∈ N

(1.1)

which are the off-diagonal absolute-row sums of A and the ratios of the sums over the corresponding diagonal entries, respectively. As a convention, we let ti (A) = +∞ if aii = 0 (even when ri (A) = 0). For any particular i, if ti (A) < 1, then the ith row is diagonally dominant, and the smaller ti (A) is, the stronger the corresponding diagonal dominance will be. If ti (A) < 1 for all 1 ≤ i ≤ n, i.e., [maxi ti (A)] < 1, then we say that A is strictly diagonally dominant; if, however, [maxi ti (A)] ≤ 1, A is diagonally dominant. For the opposite, if ti (A) > 1, then the ith row is off-diagonally dominant, and the bigger ti (A) is, the stronger the corresponding off-diagonal dominance will be. If ti (A) > 1 for all 1 ≤ i ≤ n, i.e., [mini ti (A)] > 1, then we say that A is strictly off-diagonally dominant; if, however, [mini ti (A)] ≥ 1, then we say that A is off-diagonally dominant. Often mini ti (A) ≤ 1 ≤ maxi ti (A) for a given matrix A. In terms of the notation ti ( · ), A is a GDDM if and only if maxi ti (AD) < 1 for some positive diagonal D which is usually unknown and can only be obtained by nontrivial computations. In the case that A is not a GDDM, such a D doesn’t exist, but it can be shown that there exists a positive diagonal D such that mini ti (AD) ≥ 1, provided A is irreducible (Theorem 2.2 below). In summary, when A is irreducible, we can always find some positive diagonal D such that either maxi ti (AD) < 1 when A is a GDDM or mini ti (AD) ≥ 1 when A is not a GDDM. The goal of this paper is to seek a positive diagonal D such that either mini ti (AD) ≥ 1 or maxi ti (AD) < 1 (or just less than or equal to 1). This is a realizable goal for a GDDM A or an irreducible non-GDDM A, as we just mentioned in the previous paragraph. To achieve this goal, we iteratively either decrease maxi ti (A) or increase mini ti (A) through updating A by postmultiplying it by some positive diagonal matrices. But when should we work to decrease maxi ti (A) or increase mini ti (A)? For that, we propose to use [min ti (A)] · [max ti (A)] (1.2) i

i

as an indicator as to how likely A can be transformed to AD that is diagonally dominant or offdiagonally dominant. If the product is less than or equal to 1, then we may think that mini ti (A) “dominates” maxi ti (A) and thus it is reasonable for us to predict that A tilts towards being a 2

GDDM. In such a case, we should work to decrease [maxi ti (A)]. On the other hand, if the product is bigger than 1, then mini ti (A) is “dominated” by maxi ti (A) and thus it is reasonable for us to predict that A tilts towards being off-diagonally dominant. In such a case, we should work to increase [mini ti (A)]. Like at any circumstance, prediction can go wrong. When it does, the conventional thinking is that any work before the prediction changes is waste of effort. Miraculously, this is not the case for our self-corrective algorithms in this paper. As our convergence analysis will show, our algorithms continue to make progress towards the final answer even during iterations under incorrect predictions. For example, suppose the indicator (1.2) goes from above 1 to under 1 and then back above 1 again. This means that we need to switch gear from pushing A towards being offdiagonally dominant, i.e., increasing mini ti (A), to pushing it towards being diagonally dominant, i.e., decreasing maxi ti (A) and then to pushing it towards being off-diagonally dominant again. A natural question is that: do we degrade mini ti (A) during the time when we work to decrease maxi ti (A)? If we did, the iterative process could oscillate without making progress. Luckily, such a scenario doesn’t happen for our algorithms. We will prove that mini ti (A) at the end of the step when the indicator switches back to above 1 is strictly bigger than the one right before the indicator goes under 1. The similar statement can be said for the case when the indicator goes from under 1 to above 1 and then back to under 1 again. Therefore our self-corrective algorithms always make progress regardless what the indicator says. That is the most favorable feature of our self-corrective algorithms – never stopping making progress towards the final answer. Bringing a GDDM A to strictly diagonally dominant AD has important numerical consequences, especially for linear systems Ax = b with a GDDM A. Such linear systems can be transformed into strictly diagonally dominant linear systems. With care, the latter can be solved more accurately [16, 17, 18] than not knowing they are strictly diagonally dominant systems. Symmetric diagonally dominant linear systems can be provably solved in nearly-linear time [19]. The rest of this paper is organized as follows. In section 2, we introduce some standard definitions and basic results. In section 3, we discuss two recent iterative methods that were tested to work better than other existing ones according to the literature. We propose our self-corrective iterative algorithms in section 4 and analyze their convergence behavior in section 5. In section 6, we present numerical examples and comparison results to demonstrate the effectiveness of our algorithms. Finally, we give a few concluding remarks in section 7. 2. Preliminaries First, we introduce a few standard definitions [1, 3]. For A ≡ [aij ] ∈ Cn×n , the set of all n × n complex matrices, set N− (A) = {i ∈ N : |aii | < ri (A)},

N0 (A) = {i ∈ N : |aii | = ri (A)}, N+ (A) = {i ∈ N : |aii | > ri (A)}.

(2.1a) (2.1b) (2.1c)

Evidently, N = N0 (A) ∪ N+ (A) ∪ N− (A). In , or simply I if its size is clear from the content, is the n × n identity matrix. For a vector x, diag(x) denotes the diagonal matrix with the entries of x on its main diagonal. For a matrix B and a vector x, B > 0 and x > 0 mean that they are entrywise positive, and similarly, B ≥ 0 and x ≥ 0 mean that they are entrywise nonnegative. 3

Definition 2.1. Suppose A ∈ Cn×n .

1. A is diagonally dominant if N0 (A)∪N+ (A) = N, and strictly diagonally dominant if N+ (A) = N. A is off-diagonally dominant if N0 (A) ∪ N− (A) = N, and strictly off-diagonally dominant if N− (A) = N. 2. A is said to be a generalized strictly diagonally dominant matrix (GDDM) if there exists a positive diagonal matrix3 D such that AD is strictly diagonally dominant. 3. A is reducible if there exists a n × n permutation matrix Π such that ΠAΠ =

k

T

n−k



k

n−k

A11 0

A12 A22



,

for some 1 ≤ k < n, where Π T is the transpose of Π. It is irreducible if it is not reducible. 4. A is an irreducibly diagonally dominant matrix if it is irreducible, N0 (A) ∪ N+ (A) = N, and N+ (A) 6= ∅, i.e., |aii | ≥ ri (A) for all i ∈ N and at least one of them is strict. Definition 2.2. Suppose A := [aij ] ∈ Rn×n , the set of all n × n real matrices.

1. A is called a Z-matrix if aij ≤ 0 for i 6= j. 2. Any Z-matrix A can be expressed as A = cI − B, where c > 0 and B ≥ 0. A is called an M-matrix if c ≥ ρ(B), and a nonsingular M-matrix if c > ρ(B), where ρ(B) is the spectral radius of B.

Definition 2.3. The comparison matrix of A ∈ Cn×n is denoted by M(A) := (mij ) and defined by ( |aij |, for i = j, mij = −|aij |, for i 6= j. A ∈ Cn×n is called an H-matrix if M(A) is an M-matrix, and a nonsingular H-matrix if M(A) is a nonsingular M-matrix. Next, we summarize several basic results that will be useful later. They may be found in the classical books [1, 20, 3]. We prove them here for self-containedness. Theorem 2.1. Suppose A ∈ Cn×n . (a) A is a GDDM if and only if AD is a GDDM for any given positive diagonal D. (b) A is a GDDM if and only if it is a nonsingular H-matrix. (c) An irreducible diagonally dominant matrix is a GDDM. (d) If A is a GDDM, then N+ (A) 6= ∅, i.e., A has at least one strictly diagonally dominant row. (e) Suppose A is irreducible and N− (A) = ∅. A is a GDDM if and only if N+ (A) 6= ∅. (f) If A is a GDDM, then aii 6= 0 for all i ∈ N. 3 The

requirement that D has positive diagonal entries is not essential, but has been imposed in the literature.

4

Proof. Items (a) and (f) are rather obvious by definition. For item (b), if A is a GDDM, then there exists a positive diagonal D such that AD is strictly diagonally dominant. So M(AD) = M(A)D is strictly diagonally dominant and thus M(A) is a nonsingular M-matrix, or equivalently, A is a nonsingular H-matrix. On the other hand, if A is a nonsingular H-matrix, then M(A) is a nonsingular M-matrix. Therefore there exists an n-vector x > 0 such that M(A)x > 0. It can be verified that AD is strictly diagonally dominant, where D = diag(x). Item (c) holds because any irreducible diagonally dominant matrix is nonsingular [3, p.23] and thus M(A) is a nonsingular M-matrix. So A is a nonsingular H-matrix. Now use item (b) to conclude the proof. For item (d), suppose that A is a GDDM, then AD is strictly diagonally dominant for some positive diagonal D = diag(d1 , . . . , dn ). Then M(A)x > 0, where x = [d1 , . . . , dn ]T . Let di be the smallest entry of x. Then |aii |di −

X j6=i

|aij |dj > 0

⇒

|aii | >

X j6=i

|aij |

X dj ≥ |aij |. di j6=i

That is i ∈ N+ (A). For item (e), if N+ (A) 6= ∅, then A is a GDDM by item (c). On the other hand, if A is a GDDM, then N+ (A) 6= ∅ by item (d). Theorem 2.1(b) leads to ways for checking if a given A is a GDDM by verifying if A is a nonsingular H-matrix. The basic idea is as follows. Since M(A) is a Z-matrix and can be written as M(A) = cI − B, where c > 0 and 0 ≤ B ∈ Rn×n . By the theory for nonnegative matrices [1, 21, 3], ρ(B) is an eigenvalue and it is among one of the largest ones in absolute value. Compute it and then check if c > ρ(B) to arrive at a conclusion. But this could be a very expensive approach, e.g., in the case when B is imprimitive, B has several eigenvalues equally spaced on the circle {z ∈ C : |z| = ρ(B)} [21, pp.675-680] and thus simple methods like the power method usually diverge. If we can find a positive diagonal D such that 1. N+ (AD) = ∅, i.e., mini ti (AD) ≥ 1, then A is not a GDDM by Theorem 2.1(d), or 2. N+ (AD) = N, i.e., maxi ti (AD) < 1, then A is a GDDM by definition. In the case when A is irreducible, the condition in the second item can be weakened to min ti (AD) < max ti (AD) ≤ 1. i

i

for claiming that A is a GDDM by Theorem 2.1(c,e). Basically, what most existing algorithms, and ours in this paper included, try to do in finding whether A is a GDDM or not is to seek a positive diagonal matrix D such that AD is either (strictly) diagonally dominant, i.e., maxi ti (AD) ≤ 1, or off-diagonally dominant, i.e., mini ti (AD) ≥ 1. The question is whether for any given A ∈ Cn×n it is always possible to find such a positive diagonal matrix D. In the case when A is a GDDM, by definition we can find a positive diagonal matrix D such that maxi ti (AD) < 1. What happens when A is not a GDDM? Is there a positive diagonal matrix D such that mini ti (AD) ≥ 1? The answer is yes when A is irreducible, but depends when A is reducible. 5

Theorem 2.2. Suppose A ∈ Cn×n is irreducible and has nonzero diagonal entries. If A is not a GDDM, then there exists a positive diagonal matrix D such that mini ti (AD) ≥ 1.

Proof. Write M(A) = cI − B, where c > 0 and 0 ≤ B ∈ Rn×n . B is irreducible because A is. By Perron-Frobenius theorem [1, 21], ρ(B) is a simple eigenvalue of B with the associated eigenvector x > 0, i.e., Bx = ρ(B) x. Thus M(A) x = [c − ρ(B)]x. That A is not a GDDM implies c ≤ ρ(B). It is not difficulty to verify that D = diag(x) is what we need. That A is irreducible in general cannot be removed from the conditions of this theorem. Consider the following block diagonal matrix   A1 with irreducible GDDM A1 , and ir(2.2) A= A2 reducible non-GDDM A2 . This A cannot be a GDDM. Otherwise there would exist a positive diagonal matrix D such that maxi ti (AD) < 1. Partitioning D = diag(D1 , D2 ) in the same way as for A, we find maxi ti (A2 D2 ) < 1 which would imply that A2 would be a GDDM, a contradiction. So A is not a GDDM. But there exists no positive diagonal matrix D such that mini ti (AD) ≥ 1. The latter can also be proved by a contradictory argument. With additional conditions, we can extend Theorem 2.2 a little bit further. Theorem 2.3. Suppose A ∈ Cn×n is reducible, has nonzero diagonal entries, and admits   A11 A12 . . . A1m  A22 . . . A2m    Π T AΠ =  ..  , ..  . .  Amm

(2.3)

where Π is a permutation matrix, and all Aii are irreducible.

1. If A is a GDDM, then all diagonal blocks Aii are GDDMs. In other words, if one of the diagonal blocks Aii is not a GDDM, then A cannot be a GDDM. 2. There exists a positive diagonal matrix D such that mini ti (AD) > 1 in any one of the following two cases: (a) all diagonal blocks Aii are not GDDMs and all M(Aii ) are nonsingular; (b) Amm is not a GDDM, M(Amm ) is nonsingular, and none of Aim for 1 ≤ i ≤ m − 1 has zero rows. Proof. Without loss of generality, we may simply take Π = In . Since each Aii is a principal submatrix of A, item 1 is a corollary of [5, Lemma 2.6]. For item 2(a), by the proof of Theorem 2.2, there are positive vectors xi such that M(Aii )xi < 0 T T T for 1 ≤ i ≤ m. Let x = [xT 1 , x2 , . . . , xm ] and D = diag(x). It can be verified that M(A) x < 0 which implies AD is strictly off-diagonally dominant, i.e., mini ti (AD) > 1. For item 2(b), by the proof of Theorem 2.2, there is a positive vector xm such that M(Amm )xm < 0. Now pick arbitrarily positive vectors xi for 1 ≤ i ≤ m of sizes compatible with Aii . Since |Aim | xm > 0 for 1 ≤ i ≤ m − 1, we can pick β > 0, sufficiently large such that M(Aii )xi −

m−1 X j+1

|Aij |xj − β|Aim | xm < 0 for 1 ≤ i ≤ m − 1, 6

T T T where |Aij | is interpreted as taking entrywise absolute values. Let x = [xT 1 , x2 , . . . , βxm ] and D = diag(x). It can be verified that M(A) x < 0 which implies mini ti (AD) > 1.

3. Existing algorithms By the discussion we had in section 2, we know 1. if N+ (A) = ∅, i.e., mini ti (A) ≥ 1, then A is not a GDDM, or 2. N+ (A) = N, i.e., maxi ti (A) < 1, then A is a GDDM. In the case when A is irreducible, we only need min ti (A) < max ti (A) ≤ 1. i

i

in order to claim that A is a GDDM. But for the general case when mini ti (A) < 1 < maxi ti (A), i.e., A has some rows that are strictly diagonally dominant and some rows that are strictly offdiagonally dominant, we cannot tell immediately whether such an A is a GDDM or not by simply examining its diagonal entries and absolute-row sums ri (A) of off-diagonal entries. Some nontrivial actions must be performed on A in order to make a determination. There are several existing iterative algorithms that were designed to do the job [4, 5, 10, 15, 22, 3]. Among them, the two algorithms in [4, 10] are perhaps the best in efficiency and robustness. In this section, we will state and briefly discuss the algorithms and later in section 6 we will compare our algorithms against them. We remark that an H-matrix in [10] as well as in [4] was really meant a nonsingular H-matrix, or equivalently a GDDM. To be consistent with our later algorithms, in Algorithm 1 as well as other algorithms later we will always use the word “GDDM” rather than “H-matrix” as in [4, 10]. Remark 3.1. In Algorithm 1 and those algorithms below, A is constantly updated in place by positive diagonal matrices. If we denote A at entry and exit by Ain and Aout , respectively, then, except for the trivial case at line 2, Aout = Ain D, where D is the last D at line 7. In case A is declared a GDDM, the last D at line 7 also makes that AD is strictly diagonally dominant. In case A is declared not a GDDM, either A has a zero diagonal entry or A is off-diagonally dominant and the last D at line 7 makes that AD is off-diagonally dominant. It is possible that the while-loop may be a dead loop for some input, i.e., its loop condition cannot be unsatisfied, e.g., for the matrix (2.2). So in actual implementation, one may set a maximum allowed number of the loop iterations to avoid a dead loop situation. This last comment applies to all the algorithms in what follows. Theorem 3.1 ([10]). Suppose Algorithm 1 terminates after a finite number of iterations. Then 1. A is a GDDM if tmax < 1; 2. A is not a GDDM if tmin ≥ 1.

In 2002, L. Li [12] proposed a method for determining whether a given matrix is a GDDM or not. In 2003, Ojiro, Niki, and Usui [15] gave another method for the same task. But in 2006, Alanelli and Hadjidimos [4] came up with three counterexamples. The method in [12] fails on two of them while the algorithm in [15] fails on the third. Specifically, the algorithm in [12] does not converge for     1 0 −0.5 1 0 −0.5 0  , B = −2 1 0 , A = −0.5 1 (3.1) 0 −2 1 0 −2 1 7

Algorithm 1 [10, Algorithm B] Input: irreducible A ≡ [aij ] ∈ Cn×n ; Output: whether A is or isn’t a GDDM; 1: ti = ti (A) for i ∈ N, and tmin = mini ti ; 2: if tmin ≥ 1 or aii = 0 for some i ∈ N then 3: return that A is not a GDDM; 4: else 5: D = In , tmax = maxi ti ; 6: while tmax ≥ 1 and tmin < 1 do 7: D1 = diag(ti + 1), D = DD1 , A = AD1 ; 8: ti = ti (A) for i ∈ N; 9: tmin = mini ti , tmax = maxi ti ; 10: end while 11: if tmax < 1 then 12: return that A is a GDDM; 13: else if tmin ≥ 1 then 14: return that A is not a GDDM; 15: end if 16: end if while the algorithm in [15] fails for  1 −2 −1 0  −2 1 0 −1  . (3.2) C =  0 −1/4 1 −1/2 −1/4 0 −1/2 1   0 0 0.5 It is easy to see A = M(A) = I − 0.5 0 0  =: I − E. It can be seen that E is irreducible 0 2 0 and has three eigenvalues ρ(E) exp(2πiι/3) for i = 0, 1, 2, all with the equal absolute value ρ(E) = .7937, where ι is the imaginary unit. The power method on E will diverge unless it starts with an eigenvector. The same can be said about the matrix B here. Later in section 6 we will revisit these three matrices. Having exposed the shortcomings of these two algorithms, Alanelli and Hadjidimos [4] proposed the following Algorithm 2, intended to improve them. Alanelli and Hadjidimos [4] also established Theorem 3.2 for Algorithm 2 which is essentially the power method. 

Theorem 3.2 ([4]). Suppose A ∈ Cn×n is irreducible.

1. Algorithm 2 always terminates, i.e., the while-loop from line 6 to line 10 will exit, except possibly when det(M(A)) = 0; 2. If Algorithm 2 terminates, then its answer is correct.

Later in [5], Alanelli and Hadjidimos presented a two-stage approach to extend Algorithm 2 for reducible A, basing on the fact that A ∈ Cn×n is not a GDDM if and only if there exists at least one 8

Algorithm 2 [4, Algorithm AH] Input: irreducible A ≡ [aij ] ∈ Cn×n ; Output: whether A is or isn’t a GDDM; 1: A = [diag(A)]−1 A and rmin = min ri (A); i∈N

2: 3: 4: 5: 6: 7: 8: 9: 10: 11: 12: 13: 14: 15: 16: 17: 18:

if rmin ≥ 1 or aii = 0 for some i ∈ N then return that A is not a GDDM; else D = In , rmax = max ri (A); i∈N

while rmin ≤ 1 and rmax ≥ 1 and rmin 6= rmax do 1+rn (A) 1 (A) D1 = diag( 1+r 1+rmax , . . . , 1+rmax ); D = DD1 , A = D1−1 AD1 ; rmin = min ri (A), rmax = max ri (A); i∈N

i∈N

end while if rmin > 1 then return that A is not a GDDM; else if rmax < 1 then return that A is a GDDM; else if rmax = rmin then return that M(A) is singular and thus A is not a GDDM; end if end if

principal submatrix of A that is not a GDDM [5, Lemma 2.6]. The basic idea is to run Algorithm 2 for up to certain number of the while-loop iterations, and if no determination can be made at the end, then it extracts out the principal submatrix corresponding to indices {i : ri (A) ≥ 1} and check if the principal submatrix is not a GDDM (by running the algorithm on the extracted principal submatrix). If it is not, then A is not a GDDM (otherwise still no determination can be made). Presumably this two-stage approach can help for some difficult irreducible non-GDDMs, too. 4. New algorithms In this section, we propose three new algorithms. Theoretic convergence analysis to support the algorithms will be presented in the next section. Recall our discussion in the first paragraph of section 3. Our intuitive idea is to “correct” A iteratively, towards either diagonal dominance or off-diagonal dominance. To implement this idea, we will propose a measure as a suggestive indicator to suggest when the underlying matrix A can likely be made towards being diagonally dominant or off-diagonally dominant. By how ti (A) is defined, it is quite natural to interpret ti (A) as something that quantifies the off-diagonal dominance in the ith row or, equivalently, 1/ti (A) as something that quantifies the diagonal dominance of the row. Let p = arg min ti (A),

q = arg max ti (A).

i∈N

i∈N

9

(4.1)

Along the same line, tp (A) and tq (A) quantify the least and the most off-diagonal dominance among all n rows. Our indicator measure is then simply the product (4.2)

tp tq .

When tp tq ≤ 1, it is decided that A tilts towards being a generalized diagonally dominant matrix, and as a result we will suppress the off-diagonal dominance of all or some off-diagonally dominant rows; likewise when tp tq > 1, it is decided that A tilts towards being a generalized off-diagonally dominant matrix, and as a result we will suppress the diagonal dominance of all or some diagonally dominant rows. The next theorem shed lights on how suppression should be done. For convenience, X(i,j) refers to the (i, j)th entry of the matrix X. Theorem 4.1. Suppose A ≡ [aij ] ∈ Cn×n has nonzero diagonal entries and that N+ (A) 6= ∅ and N− (A) 6= ∅. Let J ⊂ N and set B := AD, where D is diagonal with diagonal entries ( tj (A), for j ∈ J, D(j,j) = 1, for j 6∈ J. 1. If J ⊆ N+ (A), then

( 1, ti (B) ≤ ti (A),

In particular, maxi ti (B) ≤ max{1, maxi ti (A)}. 2. If J ⊆ N− (A), then ( 1, ti (B) ≥ ti (A), In particular, mini ti (B) ≥ min{1, mini ti (A)}.

for i ∈ J, for i ∈ 6 J.

for i ∈ J, for i ∈ 6 J.

Proof. Suppose J ⊆ N+ (A). Then for each i ∈ J, ti (A) < 1, and therefore for i ∈ J, P P P P j6∈J |aij | J∋j6=i tj (A)|aij | + j6∈J |aij | J∋j6=i |aij | + ti (B) = ≤ = 1. ti (A)|aii | ti (A)|aii |

For i 6∈ J,

ti (B) =

P

j∈J tj (A)|aij |

+

|aii |

P

J6∋j6=i

|aij |

≤

P

j∈J

|aij | +

This proves item 1. Item 2 can be proven in the same way.

P

|aii |

J6∋j6=i

|aij |

= ti (A).

(4.3)

This theorem can be understood in this way. In the case of item 1, the off-diagonal dominance for all rows j ∈ J ⊆ N+ (A) is suppressed, or equivalently the diagonal dominance for these rows are increased, although possibly not strictly (examining (4.3) for conditions for strict inequality there). In the case of item 2, the opposite outcome is achieved. Thus together with the indicator measure (4.2), we propose our Self-Corrective Iteration (SCI) detailed in Algorithm 3 for determining whether a given square matrix is a GDDM or not and, as by-product, yielding a diagonal matrix D by which AD is diagonally dominant in the case when A is a GDDM. 10

Algorithm 3 Self-Corrective Iteration (SCI) Input: irreducible A ≡ [aij ] ∈ Cn×n ; Output: that A is a GDDM and diagonal D that makes AD diagonally dominant, or that A isn’t a GDDM Pn and diagonal D that makes AD off-diagonally dominant; 1: γi = j=1 |aij | and ti = γi /|aii | − 1 for i = 1, 2, . . . , n; 2: p = arg min ti ; i∈N

3:

4: 5: 6: 7: 8: 9: 10: 11: 12: 13: 14: 15: 16: 17: 18: 19: 20: 21: 22: 23: 24: 25:

if tp ≥ 1 or aii = 0 for some i ∈ N then return that A is not a GDDM; else D = In , q = arg max ti ; i∈N

while tp < 1 < tq do if tp · tq ≤ 1 then J = {i : ti < 1}; else if tp · tq > 1 then J = {i : ti > 1}; end if for each j ∈ J do γi = γi + (tj − 1) · |aij | for i = 1, 2, . . . , n; D(j,j) = D(j,j) · tj , A(:,j) = tj · A(:,j) ; end for ti = γi /|aii | − 1 for i = 1, 2, . . . , n; p = arg min ti and q = arg max ti ; i∈N

i∈N

end while if tp ≥ 1 then return that A is not a GDDM, and diagonal D that makes AD off-diagonally dominant; else if tq ≤ 1 then return that A is a GDDM, and diagonal D that makes AD diagonally dominant; end if end if

Remark 4.1. A few comments on Algorithm 3 are in order. 1. The input A is required to be irreducible. This makes sure all ti > 0 and in particular, tp > 0 always. This condition is not restrictive in applications. In general, if A is reducible, we can always preprocess A by symmetrically permuting its rows and columns to the form (2.3) such that each Aii is irreducible. Suppose that we are interested in solving some linear system Ax = b. It can be seen that the system is now degenerated into m linear systems Aii yi = bi −

m X

Aij yj

j=i+1

that can be solved sequentially from j = m, m − 1 down to 1, where we have partitioned T T T T T Π T x = [y1T , y2T , . . . , ym ] and Π T b = [bT accordingly. The question is how 1 , b2 , . . . , bm ] to solve each system accurately. As we commented before, if each Aii can be brought into 11

a strictly diagonally dominant matrix Aii Di for some positive diagonal Di , then all linear systems can be solve more accurately than without transforming Aii to Aii Di . 2. The arg min and arg max at lines 2, 6, 18 may not be unique. When that is the case, simply picking any one index among the many will do just fine. 3. The updating formula for γi at line 14 can be seen from that the absolute row sum for the ith row n X X X X |aij |. tj |aij | + |aij | = (tj − 1)|aij | + j∈J

j6∈J

j∈J

j=1

4. As before, A is constantly updated in place by positive diagonal matrices. Denote A at entry and exit by Ain and Aout , respectively. Then, except for the trivial case at line 3, Aout = Ain D, where D is the one at exit. The matrix A at lines 21 and 23 refers to Ain . At line 22, it is possible that tq = 1. This means that A is a GDDM but Aout is not strictly diagonally dominant. In applications where this is not enough, a simple modification in line 7 to tp < 1 ≤ tq as the new stopping criteria will make Aout strictly diagonally dominant. Doing so usually increases the number of iterations by 1. 5. It is possible to combine the two-stage approach in [5] and the main body of Algorithm 3 to deal with some, but not all, reducible non-GDDMs such as the one in (2.2). Since the combination is rather straightforward, we will not dwell on it in this paper but leave it to the reader.

The cost per iterative step of Algorithm 3 varies and depends on the cardinality of J for each particular step. The larger the cardinality is, the bigger the cost will be. We now introduce two variations of the algorithm, aiming at reducing per step cost. In the first variation, we force J to contain just one index, i.e., p or q, depending on the situation, and thus only one column is updated per iterative step. We name the resulting algorithm as Algorithm 3a. In Algorithm 3, replace line 9 by J = {p} and line 11 by J = {q}, respectively. for future reference. In this case, the condition at line 22 can be satisfied only for tq = 1. In the second variation, we strive to balance out the cost and the amount of suppression at lines 13-16 per iterative step. Notice that the amount of suppression is proportional to relative difference of tj from 1. Therefore, it is reasonable to limit J to contain those indices j such that tj relatively further away from 1. With this guideline in mind, we formulate Algorithm 3b. In Algorithm 3, replace line 9 by J = {i : ti tq ≤ 1} and line 11 by J = {i : tp ti > 1}, respectively. 5. Convergence analysis In the following, we present a convergence analysis of Algorithm 3. We point out that every convergence result in this section is valid for both Algorithms 3a and 3b, too. To save spaces, we will not repeat them for the two variations. Our first theorem reminds us of the Jacobi method for the symmetric eigenvalue problem: it produces the eigen-decomposition for any 2 × 2 symmetric matrix [23, 24] in one rotation. Here it is shown that Algorithm 3 will tell whether a 2 × 2 matrix is a GDDM in at most one iterative step. Note that a reducible A ∈ C2×2 is excluded from all eligible inputs of Algorithm 3. 12

Theorem 5.1. For any irreducible A ≡ [aij ] ∈ C2×2 with nonzero diagonal entries, Algorithm 3 outputs a correct answer in at most one iterative step. Proof. Without loss of generality, we may assume that A > 0 and t1 =

a12 a21 ≥ t2 = . a11 a22

Thus p = 2 and q = 1. If t2 ≥ 1 or t1 ≤ 1, then Algorithm 3 will skip lines 7-19 and go straight to line 20 and give a correct answer. Suppose that t1 > 1 > t2 . Let t˜i denote the new ti (after one iteration). Then according to t1 t2 ≤ 1 or t1 t2 > 1, we have a12 a21 t˜1 = = t1 t2 ≤ 1, t˜2 = 1, or a11 a22 a12 a21 = t1 t2 ≥ 1. t˜1 = 1, t˜2 = a11 a22

(5.1a) (5.1b)

Thus Algorithm 3 will stop with a correct answer: A is a GDDM in the case of (5.1a) or A is not a GDDM in the case of (5.1b). For the ease of presentation, we introduce a new set of notations to trace the iterations in the algorithm. A(0) is the input A-matrix and D(0) = In , and A(k) and D(k) is the A- and D-matrix (k) in the kth iterative loop after executing the for-loop at lines 13–16. Define, accordingly, ti for 1 ≤ i ≤ n, and pk and qk with p0 being the p at line 2 and q0 the q at line 6. Theorem 5.2. If Algorithm 3 terminates4 , then its output is correct. Proof. Suppose Algorithm 3 exits after m iterative steps. There are two possibilities: 1) A is not a GDDM, and 2) A is a GDDM. 1. If that A is not a GDDM is claimed, then possible exits are at line 4 and line 21. If it is at line 4, then A(0) = A has either a 0 diagonal entry or N+ (A) = ∅. By Theorem 2.1(d,f), A (m) is not a GDDM. If it is at line 21, then we have tpm ≥ 1, which means N+ (A(m) ) = ∅. By (m) Theorem 2.1(a,d), A and thus A is not a GDDM. 2. Suppose that A is a GDDM is claimed. Then the only possible exit is at line 23. We have (m) (m) two cases: tqm < 1 or tqm = 1. (m) If tqm < 1, then for all i ∈ N, ti < 1, i.e., A(m) is strictly diagonally dominant. So A(m) and thus A is a GDDM by Theorem 2.1(a). (m) (m) If tqm = 1, then A is diagonally dominant. Since the case tpm ≥ 1 has been treated earlier (m) at line 21, we must have tpm < 1 when line 23 is reached. This means N+ (A(m) ) 6= ∅. Since the input A is irreducible, A(m) is also irreducible. Thus A(m) is an irreducible diagonally dominant matrix, and so it is a GDDM by Theorem 2.1(c) and so is A by Theorem 2.1(a). This completes the proof. 4 If

it terminates, it terminates in a finite number of iterative steps.

13

The next two theorems are about the “monotonicity” property of the tp - and tq -sequences (k) generated by Algorithm 3. Note that the opposite to the condition of either tqk ≥ 1 in item 1 of (k) Theorem 5.3 or tpk ≤ 1 in item 2 of Theorem 5.3 leads to immediate termination of the algorithm, and thus is welcome and not treated in the theorem. Theorem 5.3 is in fact a simple corollary of Theorem 4.1. (k−1)

(k−1)

Theorem 5.3. In Algorithm 3, suppose tpk−1 < 1 < tqk−1 . (k−1) (k−1)

(k)

(k−1)

(k)

(k−1) (k−1)

(k)

(k−1)

(k)

1. If tpk−1 tqk−1 ≤ 1 and if tqk ≥ 1, then tqk−1 ≥ tqk ;

2. If tpk−1 tqk−1 > 1 and if tpk ≤ 1, then tpk−1 ≤ tpk . (i) (i)

This theorem implies that if tpi tqi ≤ 1 during many consecutive iterations, then the corre(i) (i) (i) sponding tqi is monotonically decreasing (possibly not strictly though). Likewise, if tpi tqi > 1 (i) during many consecutive iterations, then the corresponding tpi is monotonically increasing (again possibly not strictly). Previously, we have mentioned that we would use tp tq as the indicator for suggesting how likely the input matrix A can be made towards being diagonally dominant or off-diagonally dominant, and accordingly we should suppress the diagonal or off-diagonal dominance of selected rows. But this indicator is only an indicator and it may predict incorrectly. When that happens, i.e., after certain number of iterations, tp tq may go from above (below) 1 to below (above) 1. Experience often tells us that any wrong prediction usually comes with waste of efforts in almost any circumstance. Luckily, what we have here is a counterexample to this lesson of experience. For example, suppose (i) (i) tpi tqi ≤ 1 during many consecutive iterations and therefore during each of these iterations, A is “corrected” towards being diagonally dominant, and suppose at the end of these consecutive (i) (i) iterations, tpi tqi becomes over above 1 and therefore A is “corrected” towards being off-diagonally dominant instead for the step. The next theorem says that A after this latest step is closer to being off-diagonally dominant than the one A before the many consecutive iterations for which (i) (i) tpi tqi ≤ 1. In other words, the work that has been done during all these consecutive iterations correcting A towards being diagonally dominant is not wasted, and somehow superstitiously the process knows that the indicator is making a wrong prediction during these consecutive iterations and continuously pushes A towards being off-diagonally dominant regardless. Theorem 5.4 illustrates how tp and tq behave when tp tq stays above (below) 1 and then goes below (above) 1 and then comes back to above (below) 1. (i)

(i)

Theorem 5.4. In Algorithm 3, suppose tpi < 1 < tqi for all i such that k − 1 ≤ i ≤ ℓ. (k−1) (k−1)

(i) (i)

(ℓ) (ℓ)

1. If tpk−1 tqk−1 ≤ 1, tpi tqi > 1 for k ≤ i ≤ ℓ − 1, and tpℓ tqℓ ≤ 1, then (ℓ) tq(k−1) ≥ t(k) qk > tqℓ ; k−1 (k−1) (k−1)

(5.2)

(ℓ) (ℓ)

(i) (i)

2. If tpk−1 tqk−1 > 1, tpi tqi ≤ 1 for k ≤ i ≤ ℓ − 1, and tpℓ tqℓ > 1, then (ℓ) tp(k−1) ≤ t(k) pk < tpℓ . k−1

(5.3)

Proof. First we prove (5.2). The first inequality there is a corollary of Theorem 5.3. Also due to this theorem, (k+1) (ℓ) t(k) pk ≤ tpk+1 ≤ · · · ≤ tpℓ . 14

tq :

❜

❵

❵

❵

❜

❜

(j) tq

>

❜

(k) tq

❵

❵

❵

❜

(ℓ)

❜

(s)

tq > tq

❜

❵

❵

❵ ✛ tp tq ≤ 1 ✲ ✛ tp tq > 1 ✲ ✛ tp tq ≤ 1 ✲ ✛ tp tq > 1 ✲ ✛ tp tq ≤ 1 r r r r r r r r r r r r r s−1 s s+1 k−1 k k+1 i j−1 j j+1 ℓ−1 ℓ ℓ+1

1:

❜

tp :

❵

❵

❵

(k)

❜

❜

(ℓ)

tp < tp

❜

❵

❵

❵

❜

❜

Figure 1: Behavior of tp and tq as tp tq goes through region above 1 and below 1 and back. As soon as one of tp and tq touches 1, the iterative process stops with a claim that A is either a GDDM or not a GDDM. (k)

(k)

(ℓ)

(ℓ)

(k)

(k)

(ℓ)

(ℓ)

Therefore tqk > 1/tpk ≥ · · · ≥ 1/tpℓ ≥ tqℓ , as expected. Next we prove (5.3). The first inequality there is a corollary of Theorem 5.3, too. Also due to this theorem, (k+1) (ℓ) t(k) qk ≥ tqk+1 ≥ · · · ≥ tqℓ . Therefore tpk ≤ 1/tqk ≤ · · · ≤ 1/tqℓ < tpℓ , as expected. Combining Theorems 5.3 and 5.4, we draw Figure 1 to illustrate how tp - and tq -sequences move during the iterations. In Figure 1, the line at the center corresponds to 1, tp is always below 1, and tq is always above 1. In the shaded area with opening-up, tq is above or at the bottom-line, possibly not monotonic, and possibly goes above its first value in the region (tp tq > 1). However the first value of tq immediately after the region is guaranteed smaller than the last value of tq right before the region. Similar statements applies to the shaded area with opening-down for tp . The trend is that tp moves up to 1 and tq moves down to 1. As soon as one of them touches or crosses 1, the iterative process stops with a claim. Theorem 5.5. If A ∈ Cn×n does not have any 0 entry and M(A) is nonsingular, then Algorithm 3 terminates in a finite number of steps. (k)

Proof. By the notations we introduced right before Theorem 5.2, we have A(k) ≡ [aij ] = A(k−1) D(k) . (0)

Without loss of generality, we may assume that A is entrywise positive, i.e., aij > 0 for all i, j, (k)

and thus aij > 0 for all i, j, k. Assume, to the contrary, that Algorithm 3 doesn’t terminate in a finite number of steps, and (k) (k) thus tpk < 1 < tqk for k = 0, 1, . . .. There are three possible cases: (k) (k)

(a) there is an integer K ≥ 0 such that tpk tqk ≤ 1 for all k ≥ K; (k) (k) (b) there is an integer K ≥ 0 such that tpk tqk > 1 for all k ≥ K; (k) (k) (k) (k) (c) each of tpk tqk ≤ 1 and tpk tqk > 1 occurs for infinitely many k. 15

In both cases (a) and (b), no generality is lost if we let K = 0. Consider case (a). Let Dk ≡ diag(dk,1 , dk,2 , . . . , dk,n ) = D(1) D(2) · · · D(k) .

(5.4)

(k) (k)

Since tpk tqk ≤ 1 for all k, we know 0 < dk+1,i ≤ dk,i ≤ 1. Thus dk,i converges as k → ∞, and let di = limk→∞ dk,i and D := diag(d1 , d2 , . . . , dn ) = limk→∞ Dk . Set dk,max = max dk,i , 1≤i≤n

b k ≡ diag(dˆk,1 , dˆk,2 , . . . , dˆk,n ) = Dk /dk,max . D

(5.5)

It can be seen that 0 < dk+1,max ≤ dk,max ≤ 1; so dk,max converges, too, as k → ∞, and let b k converges to D/dmax =: D. b Write D b = diag(dˆ1 , dˆ2 , . . . , dˆn ). It can be dmax = limk→∞ dk,max . D ˆ ˆ proved that max1≤i≤n di = 1 since max1≤i≤n dk,i = 1 for all k. We claim that dˆi > 0 for 1 ≤ i ≤ n. To this end, we notice A(k) = ADk . Therefore P ˆ−1 ˆ aij0 j6=i dk,i aij dk,j −1 −1 b (k) (k) (0) b tp0 ≥ tpk ≥ ti (A ) = ti (Dk ADk ) = ti (Dk ADk ) = , (5.6) ≥ ˆ aii dk,i aii

where j0 is the integer that achieves dˆk,j0 = dˆk,max = 1. This implies minj aij dˆk,i ≥ (0) tp0 aii

⇒

minj aij dˆi = lim dˆk,i ≥ (0) > 0, k→∞ tp0 aii

(5.7)

b ≡ [ˆbij ] = limk→∞ AD b k = AD. b By what we just proved, we know B b > 0. Let as expected. Let B δ = min ˆbij ,

γ = max aii .

(5.8)

1≤i≤n

i6=j

(k)

By Theorems 5.3 and 5.4 (see Figure 1), the sequence {tqk } is non-increasing and strictly bounded (k) (k) by 1 from below. So tqk converges. But we don’t know for certain if tpk converges. Let α = lim t(k) qk ≥ 1,

β = lim sup t(k) pk .

k→∞

(5.9)

k→∞

(k)

Since tpk < 1 for all k, we have β ≤ 1. We claim β < 1. There are two possibilities: (k) (k)

(k)

(k)

(i) α > 1. Since tpk tqk ≤ 1, we have tpk ≤ 1/tqk ≤ 1/α < 1 for all k, implying β < 1; (k) (ii) α = 1. Assume, on the contrary, that β = 1. There is a subsequence of {tpk } that converges b = limj→∞ = AD b kj for which to β. Let the indices of the subsequence be {kj }. Then also B b b ti (B) = 1 for 1 ≤ i ≤ n which imply that M(B) is singular, a contradiction, since M(A) is assumed nonsingular.

So β < 1. Note

(k)

b k = ADk /dk,max = A(k) /dk,max = [a ]/dk,max AD ij

b Let which converges to B. η=

(1 − β)δ > 0, γ

ǫ1 =

1−β , 2β

ǫ2 = 16

η , 2(1 + η)

ǫ = min {ǫ1 , ǫ2 } < 1.

(5.10)

There is an integer K > 0 such that α + ǫ > t(K+1) qK+1 > α,

α + ǫ > t(K) qK > α, (K) aij /dK,max (K) aii /dK,max

≥ δ(1 − ǫ) for i 6= j,

t(K) pK ≤ β(1 + ǫ),

(5.11a) (5.11b) (5.11c)

≤ aii /dK,max ≤ γ/dK,max .

We have, by (5.10) and (5.11), t(K+1) qK+1 =

X

j6=qK+1

≤

(K+1)

aqK+1 j (K+1)

aqK+1 qK+1

X

j∈{qK+1 ,pK }

(K)

aqK+1 j (K)

aqK+1 qK+1

(K)

+ t(K) pK

aqK+1 pK (K)

aqK+1 qK+1

(5.12)

(K)

≤ tqK+1 (A(K) ) − (1 − t(K) pK )

aqK+1 pK (K)

aqK+1 qK+1

(K)

(K) ≤ t(K) qK − (1 − tpK )

aqK+1 pK (K)

aqK+1 qK+1 δ(1 − ǫ2 ) ≤ α + ǫ2 − [1 − β(1 + ǫ1 )] γ = α − η/2, (k)

a contradiction, since tqk ≥ α for all k. We remark that the equality at (5.12) holds for Algorithm 3a. Case (b) can be proved in a similar way. (k) Consider case (c). By Theorems 5.3 and 5.4 (see Figure 1), {tpk } has a convergent non(kℓ ) (kℓ ) (kℓ ) b k the same as in (5.4) increasing subsequence {tpkℓ } with tpkℓ tqkℓ ≤ 1 for all ℓ. Define Dk and D b b and (5.5). {Dkℓ } has a convergent subsequence since entrywise 0 ≤ Dk ≤ In . Without loss of b ≡ [ˆbij ] = limℓ→∞ AD b k = AD. b It generality, we may assume the subsequence is itself. Now let B ℓ can be seen that we will still have (5.6) – (5.9) with all k replaced by kℓ and k → ∞ by ℓ → ∞. In the same way as for case (a) after (5.9), we will reach a contradiction. 6. Numerical Experiments Before we get into our numerical comparisons, we present Table 6.1 which displays the flop counts per iterative step for each algorithm we presented so far, where it is assumed that A is dense. It is clear that Algorithm 3a (for which nJ = 1) has the most favorable flop counts, linearly in n, while Algorithms 1 and 2 have O(n2 ) flops. But since nJ changes at the runtime for Algorithms 3 and 3b, their per step cost is not known a priori. In Experiment 6.5 for random matrices, the average nJ for Algorithm 3 is about n/2 and about 4 for Algorithm 3b on the testing matrices there. So for the experiment, Algorithm 3 costs about 1.5n2 per step while it is only 14n for Algorithm 3b. 17

algorithm

1

flops

2n

2 2

3n

3 2

3a

3b

nJ × 3n + 2n

Table 6.1: Flop counts per iterative step, where nJ is the cardinality of J and nJ = 1 for Algorithm 3a always.

Now we present our comparison results in terms of the numbers of iterations and correctness of Algorithms 1, 2, 3, 3a, and 3b on five experiments. The first four are on “difficult” matrices in the literature, and their dimensions are between 3 and 6. The fifth experiment is for randomly generated dense and sparse matrices with dimension n = 1000. Finally in the sixth experiment, the testing matrices are taken from the University of Florida Sparse Matrix Collection [25] with sizes up to almost 90, 000. Among all tests, Algorithm 3 always takes the fewest numbers of iterations. For the large random matrices in Experiment 6.5, we calculate that Algorithm 3b actually use the fewest flops overall. Experiment 6.1. This first experiment is on two matrices A and B given in (3.1) from [4]. Our algorithms work beautifully on both matrices while the algorithm in [12] fails (see [4]). This A is also the testing matrix for Example 1 in Experiment 6.4, and B is the one for Example 2 there. All Algorithms 1, 2, 3, 3a, and 3b generate the correct claim. See Examples 1 and 2 in Table 6.2 for detail. ✸ Experiment 6.2. In this experiment, we show that Algorithm 1 may malfunction, i.e., producing an erroneous claim. Consider the matrix     2 −1 −0.5 2 −1 0 2 −1  , B = −2 7 −2 . A =  −1 −1.5 −3 3.5 −2 −2 1

A is a slight modification of [4, Example 7] in changing its third row to three times as much. Test on [4, Example 7] will be reported in Table 6.2. It can be verified that this A is singular, and hence it cannot be a GDDM. But if Algorithm 1 is applied to A, at the end of 35th iteration, it is claimed that A is a GDDM. On the other hand, A is correctly identified as not a GDDM by all other algorithms. For B, all algorithms reach the same conclusion that B is not a GDDM, but Algorithms 1 and 2 take many more iterations than Algorithm 3 and its variations do, as shown in the following table algorithm 1 2 3 3a 3b A B

28

33 29

31 3

32 3

31 3

where “-” indicates that Algorithm 1 fails to make a correct determination. Experiment 6.3. This experiment is  1 0.05  0.01 A= 0.01  0.02 0.07

for the 6 × 6 matrix

0.01 0.02 0.01 0.03 1 0.10 0.02 0.01 0.01 1 1.01 0.01 0.03 1.002 1 0.01 0.01 0.02 0.01 1 0.01 0.01 0.01 0.01 18

 0.01 0.01  0.01  0.02  0.10 1

✸

Example [4] 1 2 3 4 5(i) 5(ii) 6(i) 6(ii) 7 8 9

Alg. 1 3 1 1 5 7 8 31 36 1 0

Alg. 2 4 4 2 6 8 9 32 37 33 2 0

Alg. 3 1 1 1 4 7 8 16 19 32 1 0

Alg. 3a 2 1 2 11 7 8 30 26 33 1 0

Alg. 3b 1 1 1 10 7 8 21 19 32 1 0

GDDM? Yes No No Yes Yes No Yes No No Yes No

Table 6.2: Experiment 6.4: the numbers of iterations on examples in [4, section 5]

taken from [5]. Algorithms 1 and 2 need 14 and 15 iterations, respectively, to reach the conclusion that A is not a GDDM, while Algorithms 3, 3a, and 3b take 5, 11, and 6 iterations, respectively. ✸ Experiment 6.4. In the fourth experiment, we compare the numbers of iterations of the five algorithms on all the examples in [4, section 5]. The outcomes are given in Table 6.2, where the second to sixth columns starting from the second row downwards list the numbers of iterations by the respective algorithms to reach the conclusions in the seventh column5 . Table 6.2 clearly demonstrates that Algorithm 3 uses the least numbers of iterations among all. But it is interesting to observe that Algorithm 2 performs worse than Algorithm 1, except for [4, Example 7], despite that Algorithm 2 was designed with an intention to improve algorithms that predated it. In fact, this observation is often true in our numerous tests. ✸ Experiment 6.5. In this experiment, we perform random testing on both dense and sparse matrices with n = 1000. We basically generate Z-matrices and make sure that they are either M-matrices (thus GDDMs) or not M-matrices (thus not GDDMs). Specifically, dense GDDMs are generated by, in MATLAB notation, A = abs(randn(n)); s = max(abs(eig(A))); A = (s + .01) ∗ eye(n) − A; Change s+.01 to s-.01 for dense non-GDDMs. Similarly, sparse GDDMs are generated by A = abs(sprandn(n, n, .05)); s = max(abs(eig(full(A)))); A = (s + .01) ∗ speye(n) − A; Again change s+.01 to s-.01 for sparse non-GDDMs. For each type, we generate 10 random matrices to test. Table 6.3 lists the average numbers of iterations of each algorithm on the 10 random matrices of each type. It clearly shows that Algorithm 3 consistently outperforms Algorithms 1 "

# √   1 2 2(1 + i) −4 2 matrix for [4, Example 8] is A = . It can be verified that M(A) = √ 1 −4 4 2(1 − i 3) which is a nonsingular matrix, but it was claimed that “M(A) IS SINGULAR” in [4]. One plausible explanation would be that A was misprinted in [4] since in our testing, Algorithm 2 ran correctly on this example. 5 The

√ 1+i 3 √ 4 2(1−i) 8

19

Alg. 1 13.0 13.1 13.5 13.9

A dense dense sparse sparse

Alg. 2 14.0 14.1 14.5 14.9

Alg. 3 11.5 11.1 11.9 11.7

Alg. 3a 3022.5 3092.8 3564.3 4801.8

Alg. 3b 989.4 997.3 1112.0 1124.5

GDDM? Yes No Yes No

Table 6.3: Experiment 6.5: the average numbers of iterations on random examples over 10 runs for each matrix type

A

Alg. 3

dense(Y) dense(N)

1.89 · 10 1.75 · 107

sparse(Y) sparse(N)

7

Alg. 3a

Alg. 3b

1.51 · 10 1.54 · 107

1.27 · 107 1.37 · 107

7

1.33 · 107 1.38 · 107

1.51 · 107 2.03 · 107

1.19 · 107 1.20 · 107

Table 6.4: Experiment 6.5: the average total costs in flops for each matrix type

and 2 in the numbers of iterations taken for decision-making. Since Algorithm 3 does not update all columns per iteration, unlike Algorithms 1 and 2, its per cost step, estimated at 1.5n2 for the dense case (average nJ ≈ n/2), is no higher than the latter two. Therefore Algorithm 3 uses fewest flops as well. On the other hand, the average numbers of iterations by Algorithms 3a and 3b are substantially larger than the other three algorithms. But we argue that these large average numbers of iterations can be deceptive in the sense that they do not necessarily translate into Algorithms 3a and 3b being more expensive. This is due to their low per step cost because of smaller nJ (see Table 6.1 for dense matrices). To get a sense of how small nJ can be, we list in the following table the average numbers nJ by Algorithms 3a and 3b, where “Y” and “N” indicates that the corresponding columns are for GDDMs or non-GDDMs. algorithm

dense(Y)

dense(N)

sparse(Y)

sparse(N)

3 3b

547.2 3.6

525.7 3.9

503.3 3.9

528.0 3.9

By data in this table, Tables 6.1 and 6.3, we can calculate the average total costs in flops for Algorithms 3, 3a, and 3b for the dense case. For the sparse case, we may estimate the average √ number of nonzero entries per column as .05 n and thus we can use this simple model √ (2 + .05 )nJ n + 2n for average flops per iterative step, where .05 is the sparse density we used in generating the sparse examples. Now we can calculate the average total costs in flops in Table 6.4. According to this table, Algorithm 3b outperforms Algorithms 3 and 3a. Moments ago, we argued that Algorithm 3 is cheaper than Algorithms 1 and 2. Putting all together, we conclude that Algorithm 3b stands out on the top in this random matrix experiment. We observed that often matrices AD at exit by Algorithms 3 and 3b are strictly (off-)diagonally dominant, i.e., either mini ti (AD) > 1 or maxi ti (AD) < 1 while those by Algorithm 3a are just (off-)diagonally dominant, i.e., either mini ti (AD) ≤ 1 or maxi ti (AD) ≤ 1 with equality. ✸ 20

A Trefethen 2000 Trefethen 20000b finan512 662 bus pesa poisson3Db rdb5000 rail 20209

dim 2,000 19,999 74,752 662 11,738 85,623 5,000 20,209

Alg. 1 4 2 4 115 209 2 17 37

Alg. 2 5 3 116 219 18 38

Alg. 3 4 2 3 58 72 2 18 25

Alg. 3b 6 3 16 60 127 2 54 26

GDDM? Yes Yes Yes Yes No No No No

Table 6.5: Experiment 6.6: matrices from the University of Florida Sparse Matrix Collection [25], where “-” indicates that the process takes too long (exceeding the preset maximum number of iterations). On this set, Alg. 3a takes too long to finish.

Experiment 6.6. In this experiment, we test matrices, with sizes up to almost 90, 000, from the University of Florida Sparse Matrix Collection [25]. They include matrices that are GDDM and those that are not. Table 6.5 reports typical results on several examples. What we can see here is that Algorithms 1, 3, and 3b work rather well on real life examples, large and small, GDDM and no GDDM. Among them, Algorithm 3 needs the fewest numbers of iterations while Algorithm 3b performs fairly well, too. Results by Algorithm 3a are disappointing on this set of examples, taking too long to finish. This suggests that updating just one column per iteration, although works fairly well on the small artificial examples in Experiments 6.1 – 6.4 and the random examples in Experiment 6.5, can hinder convergence in large practical examples. 7. Concluding remarks We presented three self-corrective iterative (SCI) algorithms for determining if a square matrix is a generalized strictly diagonally dominant matrix or not. Although they share the same goal, by design Algorithms 3 and 3a are opposite in per step cost: the former costs more due to that all possible columns which, if updated, can help to reduce maxi ti (A) or increase mini ti (A) are updated, while the latter costs less due to choosing just one significant column to update among all possible columns which, if updated, can help to achieve the same objective. Algorithm 3b strikes a balance between the two. Consequently, as we observed in numerous tests, Algorithms 3 takes the fewest iterations among our three algorithms, Algorithm 3a takes the most, and Algorithm 3b takes the number in between. Because of the differences in the per step cost, the number of iterations is no longer a good suggestion in telling which algorithm is cheaper and which is more expensive. For the random testings in Experiment 6.5, we calculated that Algorithm 3b uses the fewest flops among the three. However, the performance of Algorithm 3a on Experiment 6.6 is rather disappointing. We compared our algorithms against two existing ones from [4, 10] – Algorithms 1 and 2. They are perhaps the most efficient and robust iterative methods in the literature. But our tests suggest that Algorithm 2 often performs worse than Algorithm 1, despite that Algorithm 2 was created with an intention to improve algorithms that predated it. Overall, Algorithm 3 performs the best and Algorithm 3b performs comparably. Our analysis reveals that our algorithms have an elegant convergence behavior as illustrated by Figure 1. We also showed that our algorithms will terminates in finite many steps if M(A) is 21

nonsingular and A has no zero entries. We conjecture that A having no zero entries is too strong and possibly could be removed, but a proof eludes us. We point out that all our numerical tests including those not reported here, many of which have zero entries, terminated successfully in finite numbers of steps. While our Algorithms 3 and 3b work so well in our tests, beating the best ones in the literature, we don’t have good quantitative estimates on their rates of convergence. This will be one of the subjects of our future research. The definition of GDDM and Theorem 2.2 guarantee that the exit conditions in our algorithms are achievable. However, if they are somehow applied to reducible A regardless, the algorithms may not exit, e.g., for those reducible GDDM such as the one in (2.2). For the types of reducible GDDMs such as the ones in Theorem 2.3(2), the exit conditions can be satisfied, however. Hence extra care must be taken in dealing with reducible matrices A for which the two-stage approach in [5] can be helpful (see our discussion at the end of section 3). References [1] A. Berman, R. J. Plemmons, Nonnegative Matrices in the Mathematical Sciences, SIAM, Philadelphia, 1994. [2] R. A. Horn, C. R. Johnson, Topics in Matrix Analysis, Cambridge University Press, Cambridge, 1991. [3] R. Varga, Matrix Iterative Analysis, Springer-Verlag, Berlin, 2000. [4] M. Alanelli, A. Hadjidimos, A new iterative criterion for H-matrices, SIAM J. Matrix Anal. Appl. 29 (1) (2007) 160–176. [5] M. Alanelli, A. Hadjidimos, A new iterative criterion for H-matrices: the reducible case, Linear Algebra Appl. 428 (11-12) (2008) 2761–2777. [6] R. Bru, I. Gimnez, A. Hadjidimos, Is A ∈ Cn,n a general H-matrix?, Linear Algebra Appl. 436 (2) (2012) 364–380. [7] T.-B. Gan, T.-Z. Huang, Simple criteria for nonsingular H-matrices, Linear Algebra Appl. 374 (2003) 317–326. [8] M. Harada, M. Usui, H. Niki, An extension of the criteria for generalized diagonally dominant matrices, Intern. J. Comput. Math. 60 (1-2) (1996) 115–119. [9] T.-Z. Huang, A note on generalized diagonally dominant matrices, Linear Algebra Appl. 225 (1995) 237–242. [10] T. Kohno, H. Niki, H. Sawami, Y. Gao, An iterative test for H-matrix, J. Comput. Appl. Math. 115 (12) (2000) 349–355. [11] B. Li, L. Li, M. Harada, H. Niki, M. J. Tsatsomeros, An iterative criterion for H-matrices, Linear Algebra Appl. 271 (13) (1998) 179–190. [12] L. Li, On the iterative criterion for generalized diagonally dominant matrices, SIAM J. Matrix Anal. Appl. 24 (1) (2002) 17–24. 22

[13] L. Li, H. Niki, M. Sasanabe, A nonparameter criterion for generalized diagonally dominant matrices, Intern. J. Comput. Math. 71 (2) (1999) 267–275. [14] J. Liu, A. He, An interleaved iterative criterion for H-matrices, Applied Mathematics and Computation 186 (1) (2007) 727–734. [15] K. Ojiro, H. Niki, M. Usui, A new criterion for the H-matrix property, J. Comput. Appl. Math. 150 (2) (2003) 293–302. [16] M. Dailey, F. Dopico, Q. Ye, A new perturbation bound for the LDU factorization of diagonally dominant matrices, SIAM J. Matrix Anal. Appl. 35 (3) (2014) 904–930. [17] M. Dailey, F. Dopico, Q. Ye, Relative perturbation theory for diagonally dominant matrices, SIAM J. Matrix Anal. Appl. 35 (4) (2014) 1303–1328. [18] P. Koev, F. Dopico, Perturbation theory for the LDU factorization and accurate computations for diagonally dominant matrices, Numer. Math. 119 (2011) 337–371. [19] D. A. Spielman, S.-H. Teng, Nearly-linear time algorithms for preconditioning and solving symmetric, diagonally dominant linear systems, SIAM J. Matrix Anal. Appl. 35 (3) (2014) 835–885. [20] M. Fiedler, Special Matrices and Their Applications in Numerical Mathematics, 2nd Edition, Dover Publications, Inc., Mineola, New York, 2008. [21] C. D. Meyer, Matrix Analysis and Applied Linear Algebra, SIAM, Philadelphia, 2000. [22] R. S. Varga, On recurring theorems on diagonal dominance, Linear Algebra Appl. 13 (12) (1976) 1–9. [23] J. Demmel, Applied Numerical Linear Algebra, SIAM, Philadelphia, PA, 1997. [24] B. N. Parlett, The Symmetric Eigenvalue Problem, SIAM, Philadelphia, 1998. [25] T. Davis, Y. Hu, The University of Florida sparse matrix collection, ACM Trans. Math. Software 38 (1) (2011) 1:1–1:25.

23