Scaling by Binormalization - Semantic Scholar

Scaling by Binormalization Oren E. Livne

∗†

Gene H. Golub

∗‡

June 3, 2003 Abstract We present an iterative algorithm (BIN) for scaling all the rows and columns of a real symmetric matrix to unit 2-norm. We study the theoretical convergence properties and its relation to optimal conditioning. Numerical experiments show that BIN requires 2–4 matrix-vector multiplications to obtain an adequate scaling, and in many cases significantly reduces the condition number, more than other scaling algorithms. We present generalizations to complex, non-symmetric and rectangular matrices. Key words. Optimal scaling, BIN, Gauss-Seidel-Newton, relaxation methods.

1

Introduction

We consider the problem of binormalizing a real symmetric algebraic system of equations Ax = b, (1) where A = AT is an n × n matrix. Namely, we would like to diagonally scale A by D := diag(d1 , . . . , dn ) such that the obtained matrix A˜ := DAD will remain symmetric and satisfy X A˜2 = 1, i = 1, . . . , n. (2) ij

j

∗

The Program for Scientific Computing and Computational Mathematics, Computer Science Department, Gates 2B, Stanford University, Stanford, CA 94305-9025. Phone: +1-650-725-4163. Fax: +1-650-725-7411. † Email address: [email protected] ‡ Email address: [email protected]

1

BINORMALIZATION

2

In many matrix problems arising in scientific and engineering applications, binormalization is beneficial and even necessary. The main motivating reasons are as follows: • Scaling may be regarded as a form of pre-conditioning: the condition ˜ may be much smaller than κ(A). That is, the numerinumber κ(A) cal stability properties of the matrix are improved [4, 12, 13, 18, 22]. For instance, a theorem by Forsythe and Straus [13] ensures optimal conditioning by diagonal scaling for a wide class of matrices. • Different rows are now “comparable”: let x˜ be an approximation to the solution x of (1), e := x − ξ the error vector, r := b − Aξ = Ae ˜ the normalized residual vector ; then the residual vector, and r˜ := Ae k˜ rk indicates the true magnitude of “algebraic connections” between unknown variables, whereas krk does not [10, p. 4]. This is of interest for instance in geometric and algebraic multigrid algorithms [7],[9, §11],[10]. In particular, “algebraically smooth” errors e (that are approximated on coarse levels of these algorithms) are characterized by the relation k˜ r k  kDek. Binormalization algorithms have been developed and investigated in many papers (e.g. [4, 13, 18, 22] and their references). Most studies refer to scaling of the rows and/or columns in the 1-norm by some iterative “equilibrium” process. In this paper we set for ourselves three main goals, as follows: 1. Develop a fast iterative binormalization algorithm that calibrates the 2-norm of A’s rows (i.e., solves (2)), test it against practical numerical examples, and observe its convergence and conditioning properties. 2. Theoretically study the algorithm’s convergence properties and its possible relation to optimal conditioning of A. 3. Develop analogous algorithms for non-symmetric and rectangular matrices. Regarding the first (and main) aim, we propose an algorithm for efficiently solving (2), the BIN algorithm. BIN is based on the Gauss-Seidel-Newton method, whose convergence properties have been well analyzed [21, §7.4]. Moreover, BIN’s asymptotic convergence factor is in many empirical experiments between 0.1–0.25; hence, only a few iterations are needed to obtain an accurate binormalization.

BINORMALIZATION

3

BIN can efficiently handle sparse matrices, and yields a binormalization up to a tolerance ε in O(nz log( 1ε )) operations, where nz is the number of non-zeros in A. The algorithm’s implementation is quite simple, and it is easily and efficiently parallelizable. Furthermore, it turns out that for many matrices BIN can efficiently reduce κ(A). BIN’s main drawback is that it is 2–4 times more expensive than “trivial” ˜ tends to be smaller for BIN, (e.g. diagonal scaling) algorithms. However, κ(A) and since most binormalizations are done once-and-for-all and need not be repeated when resolving the problem (e.g. with a different b), BIN turns out useful. Second, we show that BIN converges locally, provided that a solution to (2) exists and that the initial guess is sufficiently close to it. Although we have not proved global convergence and there is no quantitative bound on the convergence factor, in practice, the first few iterations converge much faster and needed for approximate binormalization, which is by all means adequate for most problems. We also discuss possible reasons for BIN’s reduction of the condition number; theorems on optimal conditioning by diagonal scaling [13, 22] may indicate similar properties of BIN. Third, we generalize BIN to complex-valued, non-symmetric and rectangular matrices. The resulting algorithm (denoted NBIN) is typically asymptotically slower; however, it seems to retain all other “good properties” of BIN: fast convergence of the first few iterations and reducing κ(A). Our conclusion is that BIN is a novel iterative binormalization algorithm that may be useful in many applications. Further research may be directed to rectifying its theoretical convergence analysis, its precise relation to optimal conditioning and developing non-symmetric/rectangular theory for NBIN. The paper is organized as follows. In §2 we formulate the binormalization equations to be solved. In §3 we describe BIN, an iterative algorithm for solving these equations. Then, in §4 we discuss theoretical convergence analysis for BIN and a discussion of its relation to optimal conditioning. In §5, we discuss BIN’s generalizations to complex-valued, non-symmetric and rectangular matrices. Finally, §6 includes concluding remarks and suggestions for future research. Appendixes A and B.2 contain MATLAB codes and numerical experiments for both the symmetric and non-symmetric algorithms.

4

BINORMALIZATION

2 2.1

Formulation Notation

Matrices are denoted by capital roman letters (e.g. A, B), vectors by lowercase roman (e.g. f, g). Greek letter denote vector functions (e.g. β = β(x)); calligraphic letters denote matrix functions (e.g. B = B(x)). Let A be an n × n matrix. We assume that n ≥ 2, that A is real symmetric and that A has no zero rows, i.e. for any i there exists j such that Aij 6= 0 (otherwise (2) is insoluble). We denote D := diag(d1 , . . . , dn ) for any vector d ∈ Rn ; The l2 norm of a vector d ∈ Rn is defined by kdk2 := P

n X i=1

!1 2

|di |2

.

P

Unless otherwise stated, j denotes nj=1 ; and the matrix B refers to the componentwise-square of the A under consideration, i.e. the Schur product A ◦ A [5, p. 30]. Namely, Bij := A2ij , We also refer to

P

j

i, j = 1, . . . , n.

Bij as the ith row sum of B. P

Definition 1 A is called scaled if j A2ij are equal for all i = 1, . . . , n. Equivalently, B is called scaled if its row sums are equal.

2.2

The Binormalization Equations

We would like to transform A to a scaled symmetric matrix A˜ by scaling rows and corresponding columns. Thus, A˜ = DAD for some d1 , . . . , dn ∈ R. The unknown row scaling factors di are chosen to satisfy the equations X j

A˜2ij =

X j

d2i A2ij d2j ≡ const.,

i = 1, . . . , n.

(3)

Now, let xi := d2i > 0 and B(x) := diag(β1 (x), . . . , βn (x)),

βi = βi (x) :=

X j

Bij xj ,

(4)

5

BINORMALIZATION

for i = 1, . . . , n. Then we seek the positive solution x = x˜ of the autonomous system B(x)x = be, e := (1, . . . , 1)T , (5) |

{z n

where b > 0 is arbitrary (e.g. b = 1). Equivalently, 

}



1 1 In − eeT B(x)x = be − eeT be = 0. n n

(6)

Definition 2 A is called scalable if there exists a positive solution to (5), or equivalently, to (6). In most of this paper (§2.3 - §4) we will refer to the form (6) of the equations. Given an approximation x to x˜, we denote by s(x) :=

2 1 X xk βk − β n k

!1

2

,

(7)

the standard deviation of row sums, where β := (β1 , . . . , βn )T and β := 1 T x β. n

2.3

Existence and Uniqueness of Solutions

Equation (6) is an n×n system of quadratic equations in x1 , . . . , xn for which we seek a positive solution. A solution to (6) does not always exist. For instance, the matrix A=

1 1 1 0

!

(8)

is unscalable. Indeed, B = A and (6) can be written explicitly as x21 + x1 x2 = 1 x1 x2 = 1, which does not have a solution. (8) can be generalized to n × n unscalable matrices, in view of the following theorem. Definition 3 A matrix A is called diluted if there exists 1 ≤ i ≤ n such that Aii 6= 0 and Akj = 0 for all k 6= i, j 6= i.

6

BINORMALIZATION

Note that for a diluted matrix Aij 6= 0 for all i 6= j, otherwise A has a zero row. Theorem 1 If a matrix is scalable, it is either (a) not diluted, or (b) n = 2 and ! 0 a A= (9) a 0 for some a 6= 0.

Proof. First, if n = 2 and A is diluted and not of the form (9), then without loss of generality A11 6= 0, A22 = 0. (2) implies that A˜12 = A˜21 , hence row 1 cannot be scaled. Secondly, let n > 2 and let A be diluted with a non-zero row i and j 6= i. Then Aij 6= 0; moreover, Aij is the only non-zero entry in row j, therefore A˜ij = 1. Thus X A˜2ij = n − 1 + A˜2ii > n − 1 > 1, j

contradiction. Thus A is unscalable.  Even if a positive solution to (6) exists, it need not be unique. Note that if x solves (6) then cx is also a solution for any c > 0. Even if we look at 2 (2), the matrix (9) can be scaled by any x of the form (x1 , xa2 ). Similar n × n examples can be easily constructed. However, a unique solution may exist for a wide class of matrices. Evidence and intuition may be gained by the numerical experiments of App. A. Whenever there is a solution, be it unique or not, our main goal is its efficient numerical computation, described in §3. Finally, the notion of scalability can be extended to “ε-scalability” (see §6) so that a yet larger set of matrices can be (approximately) binormalized.

3 3.1

The BIN Algorithm GSN Method

The system (6) is a special case of a general system Φ(x) = 0 of n nonlinear equations in n unknowns, φ1 (x1 , . . . , xn ) = 0 .. . φn (x1 , . . . , xn ) = 0.

(10)

7

BINORMALIZATION

We focus on the Gauss-Seidel-Newton (GSN) iterative method for numerically solving (10) [21, §7.4], although other methods can be considered as well; see §5. Let x(k) be the approximation to the solution of (10) after k GSN iterations. A GSN iteration (“sweep”) consists of n steps; at step i, we solve (k+1) (k+1) (k) φi (x1 , . . . , xi−1 , xi , xi+1 , . . . , xn(k) ) = 0 (11) (k+1)

for xi , and set xi = xi . Equation (11) can be solved by one or a few steps of Newton’s method.

3.2

General Outline

The BIN algorithm is an application of the general GSN method to solving (6). The outline for the BIN code follows. ALGORITHM BIN Input: A, n × n real symmetric matrix; d(0) ∈ Rn , an initial guess for d; TOL: tolerance parameter ˜ the scaled A; d ∈ Rn , the scaling factors vector. Output: A, 1. Set B ← A ◦ A;

x←



 (0) 2 d1

,...,



d(0) n

2  T

.

2. Compute s(x) from (7). 3. While ((s(x) > T OL · β(x)) and (maximum #sweeps not reached)) do steps 4-6. 4. For i = 1, . . . , n do step 5. 5. Solve equation i of (6) (hereafter denoted (6)i ) for xi , keeping all other xj , j 6= i fixed at their current values. 6. Compute s(x) from (7). 7. Set di ← β

− 12

xi , i = 1, . . . , n,

A˜ ← DAD.

In other words, in each BIN sweep we loop over rows (and corresponding columns). Each row and corresponding column are scaled so that after the scaling, the row sum equals the average of all row sums (the corresponding

8

BINORMALIZATION

statement automatically holds for columns). This process will be referred to as a “scaling step”, whereas an entire iteration consisting of n scaling steps will be denoted by a “scaling sweep”. The sweeps are terminated when the standard deviation of the row sums xi βi is small compared with their average β. The scaling accuracy is gauged by the tolerance parameter T OL.

3.3

Scaling Step Specification

In each scaling step we satisfy (6)i , thereby replacing the current xi with its updated value, denoted by xi . Let xi , βi , β be the values before the ith scaling step. Note that (6)i can be written as (n − 1)xi βi =

X

xk βk .

k6=i

Substituting xi for xi , we get 

(n − 1)xi Bii xi +

X j6=i



Bij xj  =

X k6=i



xk 

X j6=i



Bkj xj + Bki xi  .

Rearranging the terms yields a quadratic equation in xi , namely, c2 x2i + c1 xi + c0 = 0, (12) 2 c2 := (n − 1)di , c1 = (n − 2)(βi − Bii xi ), c0 = −Bii xi + 2βi xi − nβ. Hence, in step 5 of the BIN algorithm we replace xi by the positive solution xi (12) . In the following section we show when such a solution must exist and how to handle exceptions.

3.4

Detecting Unscalability

First, c0 ≤ 0. Indeed, c0 = βi xi − Bii x2i + (nβ − βi xi ) = xi =

X j6=i

=

X j6=i

Bij xi xj − Bij xi xj −

XX

X j6=i

Bij xj −

X

xk βk

k6=i

Bkj xk xj

k6=i j

XX k6=i j6=i

Bkj xk xj −

X k6=i

Bki xk xi = −

XX k6=i j6=i

Bkj xk xj ≤ 0,

9

BINORMALIZATION

where the last equality used the symmetry of B. The case c0 = 0 occurs if and only if Bkj = 0 for all k 6= i, j 6= i, i.e. when B is diluted. Note that if c2 = 0 then c1 > 0, otherwise B had a zero row in contradiction to our assumption in §2.1; thus there exists a solution xi to (12) that vanishes if and only if c0 = 0. If c2 > 0, then the discriminant ∆ := c21 − 4c0 c2 is non-negative. ∆ = 0 if and only if c1 = 0 and c0 = 0, but that implies B to have zero entries except its i, i element, again contradicting §2.1. Thus, ∆ > 0 and the two distinct solutions of (6) are √ √ −c + −c − ∆ ∆ 1 1 x1i = x2i = , . 2c2 2c2 Clearly, x2i is strictly negative, since c1 ≥ 0 and ∆ > 0; x1i is strictly positive, since ∆ = c21 − 4c0 c2 > c21 . We conclude that a positive solution xi to (12) exists if and only if the matrix B is not diluted. The correct numerical computation of xi is xi =

−2c0 √ , c1 + ∆

which works √ even when c2 = 0, and does not suffer from round-off errors when c1 ≈ − ∆ [15]. Note that unlike the general GSN method [21, §7.4], we need not explicitly perform Newton steps on (6)i . Consequently, if A is scalable, a positive solution to (12) always exists, as implied by Theorem 1. Note that Theorem 1 does not guarantee that if BIN can be carried out, A is necessarily scalable. At any rate, if A is diluted, the BIN procedure will encounter c0 = 0 for this i and immediately detect unscalability. A small c0 value can be used to flag near-unscalability of A.

4

Theoretical Properties

In §4.1 we prove that the BIN algorithm converges to a solution of (6) for any scalable matrix, provided with a good enough initial guess. In §4.2 we describe a procedure to compute the asymptotic convergence factor. Although a qualitatively sharp upper bound on the convergence factor does not exist, we experimentally showed that for many test matrices convergence is rapid (see §A). Finally, §4.3 discusses the relation of BIN to optimal conditioning of a matrix.

10

BINORMALIZATION

4.1 4.1.1

Convergence Linearization of General Nonlinear Systems

The general theory of generalized linear iterative solvers for (10) [21, §10.3] states that the local convergence properties of the GSN method for (10) are identical to those of the Gauss-Seidel (GS) relaxation [26, Chap. 3] for the linearized system ∇Φ(˜ x)y = 0, (13) where x˜ is a solution of (10). Here  

∇Φ(˜ x) =  

∂φ1 (˜ x) ∂x1

...

.. . ∂φn (˜ x) . . . ∂x1

∂φ1 (˜ x) ∂xn



 .. . .  ∂φn (˜ x) ∂xn

Specifically, we have the following ([21, Theorem 10.3.3] for ω = 1): Theorem 2 Let x˜ such that Φ(˜ x) = 0. If the asymptotic convergence factor of GS for (13) is σ < 1, then there exists an open neighborhood S of x˜ such that for any initial x(0) ∈ S, GSN for (10) converges to x˜ with asymptotic convergence factor σ. 4.1.2

The Linearized (6)

Applying this theory to (6), assume A is scalable; hence, a solution to the binormalization equations (6) exists and will be denoted by x˜. For simplicity we assume that x˜ has been scaled so that all βi x˜i = 1. We define 

φ(x1 , . . . , xn ) := In −



1 T ee B(x)x. n

(14)

It can be easily shown that 



B11 x˜1 + (1 − n2 )β1 B12 x˜1 − n2 β2 . . . B1n x˜1 − n2 βn   .. .. .. .. , ∇Φ(˜ x) =  . . . .   Bn2 x˜n − n2 β2 . . . Bnn x˜n + (1 − n2 )βn Bn1 x˜n − n2 β1 (15) where βi are given by (4) for x = x˜. By Theorem 2, the convergence properties of GSN depend on the eigenvalues of the (singular) matrix H := (D − L)−1 U, where D, −L, and −U are the diagonal, strictly lower-, and

11

BINORMALIZATION

strictly upper-triangular parts of ∇Φ(˜ x). Note that if we multiply ∇Φ(˜ x) by any non-singular matrix M on the left, H is transformed into (M D − M L)−1 M U = (D − L)−1 M −1 M U = H, i.e. H remains unchanged. Note that multiplying by M on the right of ∇Φ(˜ x) −1 transforms H to M HM , whose eigenvalues are identical to H’s. We utilize ˜ := diag(˜ this fact and multiply (15) by X x1 , . . . , x˜n ) on the left. H is replaced by ˜ G := X∇φ(˜ x) =  B11 x˜21 + (1 − n2 )β1 x˜1 B12 x˜2 x˜1 − n2 β2 x˜2 . . . B1n x˜1 x˜2 − n2 βn x˜n  .. .. .. ..  . . . .  2 2 2 Bn1 x˜n x˜1 − n β1 x˜1 Bn2 x˜n x˜1 − n β2 x˜2 . . . Bnnx˜n + (1 − n2 )βn x˜n

˜ := XB ˜ X ˜ to obtain We substitute βi x˜i = 1 and B 

˜11 + 1 − B  .. G= .  ˜ Bn1 − 2 n

2 n

˜12 − 2 . . . ˜1n − 2 B B n n .. .. .. . . . 2 ˜ ˜ Bn2 − n . . . Bnn + 1 −

 2 n

  

= In −



 . 

2 T ˜ (16) ee + B. n

˜ is a symmetric stochastic matrix : all its entries are non-negative, Note that B the sum of each row is one and it has no zero rows. Thus, the eigenvalues λi (B) of B satisfy |λi (B)| ≤ 1, by Gershgorin’s theorem [24][p. 249]. Moreover, by the Perron-Frobenius Theorem [23, p. 27], if B has p eigenvalues equal in modulus to 1, then they must be all different and satisfy λp = 1. But here B is also symmetric, hence its eigenvalues are real. Thus p ≤ 2 and B may have unit-modulus eigenvalues 1 and −1, or only 1. In any case, λ = 1 is a simple eigenvalue. ˜ is a semi positive definite matrix with row sums equal Next, C := In + B to 2. Its eigenvalues satisfy 0 ≤ λ ≤ 2, and the largest is λ = 2 and simple with the corresponding eigenvector w. Let λ 6= 2 be an eigenvalue of C with eigenvector x. Then λ

X

xi =

i

Consequently,

X

λxi =

X i

i

P

i

(Cx)i =

XX i

j

Cij xj =

X

xj

j

xi = 0. Moreover, (Gx)i = (Cx)i −

2X xj = λxi , n j

X i

Cij = 2

X j

xj .

12

BINORMALIZATION

hence (λ, x) is an eigenpair of G as well. Conversely, if λ 6= 0 is an eigenvalue P of G with eigenvector x, then Cx = λx + n2 j xj . It follows that 2

X

xi =

i

X i

= λ(

 

X

X i

j



Cij  xi =

xi +

XX i

Cji xi =

j

XX i

X 2X xj . xj ) = (λ + 2) n j j

j

Cij xj =

X

(Cx)i

i

P

Thus, j xj = 0 and (λ, x) is an eigenpair of C. We conclude that the eigenvalues of C and G are identical except for the largest eigenvalue of C, λ = 2, which corresponds to λ = 0 and the same eigenvector w, for G. In particular, G is semi positive definite. Furthermore, Gii > 0 for all i if n > 2, which will be assumed for the rest of this section (BIN’s convergence for n = 2 can be separately and easily proved). 4.1.3

The Ostrowski-Reich Theorem

The Ostrowski-Reich Theorem [26, p. 83] gives necessary and sufficient conditions for the convergence of the successive over-relaxation (SOR) and Jacobi relaxations for linear algebraic systems. For the GS relaxation, it reduces to the following result. Theorem 3 Let A be an n × n symmetric matrix. Then GS relaxation for A˜ x = b converges for any initial guess x(0) if and only if A is positive definite. To prove that GS relaxation converges for G, we need to generalize this theorem to semi -positive definite matrices. Although there is not a unique solution to G˜ x = 0 (or no solution at all for G˜ x = b,b 6= 0), we can prove an analogous result to Theorem 3. Theorem 4 Let G be an n × n semi-positive symmetric matrix with Gii > 0 for all i. Then GS relaxation for G˜ x = 0 converges to a solution, for any (0) initial guess x . Proof. Let x ∈ Rn . The GS sweep [26, p. 64–65] consists of n steps. At step i, x is replaced by x + δei , where ei is the ith unit vector and δ is chosen to satisfy X Gij xj + Gii (xi + δ) = 0, j6=i

13

BINORMALIZATION

where xj ,j 6= i are at their current values (e.g. “new” for j < i and “old” for j > i, for a lexicographic ordered GS); xi is the “old” value (before the δ-update). Hence, δ=−

ri , Gii

ri := (Gx)i =

X

Gij xj .

j

ri is the dynamic residual at the point i. It follows that (x + δei )T G(x + δei ) = xT Gx + δeTi (Gx) + δxT Gei + δ 2 eTi Gei = xT Gx + 2ri δ + Gii δ 2 = xT Gx + 2(−Gii δ)δ + Gii δ 2 = xT Gx − Gii δ 2 . Because Gii > 0, this indicates that the energy norm xT Gx is strictly decreased after any GS step. G is semi-positive definite, thus xT Gx ≥ 0 for all x. A bounded monotone sequence has a limit, thus for any initial vector x, xT Gx → 0 with the GS sweeps. At convergence, xT Gx = 0, i.e. x = 0 or Gx = 0. Either way, x is a solution of Gx = 0.  In fact, a stronger result (complex-valued version, “iff” theorem, Jacobi and SOR for 0 < ω < 2, and so on) can be proved, but Theorem 4 suffices for our purposes. Theorems 2 and 4 imply the local convergence of BIN. Our summarizing theorem is given next.

4.2

Asymptotic Convergence Factor

The matrix G defined by (16) satisfies the conditions of Theorem 2. Hence we have the following. Theorem 5 Let x˜ be a positive solution of (6). Then there exists an open neighborhood S of x˜ such that for any initial x(0) , the BIN algorithm converges linearly to x˜ with convergence factor σ, where σ is the second largest eigenvalue modulus of the matrix H := (D − L)−1 U , where D, −L, and −U are the diagonal, strictly lower-, and strictly upper-triangular parts of ˜ G = In − n2 eeT + B. Proof. BIN converges in a local neighborhood of x˜ with the same convergence factor of the linear GS for Gx = 0, thanks to Theorems 2,4. It remains to discuss the relation of σ to the eigenvalues of H.

14

BINORMALIZATION

Without loss of generality, let λ1 , . . . , λN be the simple eigenvalues of H, 1 = |λ1 | > . . . |λN |, and v1 , . . . , vN the corresponding eigenvectors. Let P x(0) = nk=1 ck vk be an initial guess. Then after k sweeps on Gx = 0 we obtain x

(k)

Gx

(k)

= c 1 v1 +

n X

λkl cl vl

l=2

= c1 Gv1 + λk2 c2 v2 + o(λk2 ).

At convergence, x(k) tends to a solution of Gx = 0; hence, x(k) → c1 v1 . Moreover, kx(k) − c1 v1 k2 = O(|λ2 |k ); hence, σ = |λ2 |k is the asymptotic convergence factor.  For debugging purposes, the eigenvalues of this G may be computed (if this task is feasible), and σ may be estimated; however, there is no a priori uniform bound on the convergence factor, in view of the example of App. A.2.7. 4.2.1

Initial Guess and Stopping Criterion

Assuming a solution to the original binormalization equations (6) (otherwise G might not be meaningful), the BIN algorithm converges for any initial guess in the vicinity of a solution of (6). A convenient start, unless more information is given on the problem, is xi = 1 for all i. The solution of (10), where Φ is given by (14), is not unique; in particular, it is determined up to an multiplicative constant. Hence we can observe convergence in practice by examining kx(k) − x(k−1) k2 or the residual norm kΦ(x)k2 , after every BIN sweep. Our stopping criterion is s(x) = √1 kΦ(x)k2 < ε, where ε is typically not too small, e.g. 10−4 to 10−1 for most n applications. Every BIN sweep costs O(nz ) computer operations, where nz is the number of non-zero entries in A; or equivalently, one matrix-vector multiplication (MVM). An additional MVP is needed at initialization (step 1). As we will see, BIN rapidly converges in the first few sweeps, thus we can obtain s(x) ≤ ε in O(nz log( 1ε )) operations; typically, 1–3 sweeps are required, i.e. 2–4 MVMs.

4.3

Optimal Conditioning

Given a real symmetric matrix A, how can we choose a diagonal matrix D so that κ(DAD) is minimized? This problem is often known as the problem of

15

BINORMALIZATION

“optimal scaling” (or “optimal conditioning”). Studies of this problem can be found for instance in [4, 13, 18, 22]. For example, the following theorem shows that for symmetric positive definite matrices, “straightforward” symmetric scaling of all diagonal entries to one is “sub-optimal”, and in some cases the optimal choice. ˜ := diag(A) and A˜D := Theorem 6 Let A be symmetric positive definite, D − 21 − 12 ˜ ˜ D AD . 1. (Forsythe–Straus, [13]) Let λ the minimal eigenvalue and Λ the maximal eigenvalue of A; let xλ , xΛ be the corresponding eigenvectors. If |xλi | = |xΛi |, then

i = 1, . . . , n,

(17)

κ(A˜D ) =

min κ(DAD). diagonal If λ, Λ are simple eigenvalues, (17) is also necessary. D

2. (van der Sluis, [22]) For any symmetric positive definite A, κ(A˜D ) ≤ n min κ(DAD). D diagonal

(18)

The factor n can be replaced by p when there are at most p non-zero elements per row of A. The Forsythe-Straus theorem implies that all matrices with Property A whose diagonal entries are all one, are optimally scaled [13, p. 344]. By contrast, BIN scales all the diagonal entries of AAT to one. It cannot be better than diagonal scaling when A has Property A; however, we can hope to still be close to the optimal scaling. For instance, all 2 × 2 matrices have Property A, and it is easy to show that the BIN scaling yields in this case A˜ = cA˜D for some constant c; hence, BIN is always optimal for n = 2 (see App. A.2.1). ˜ than κ(A˜D ) for Furthermore, BIN may give better conditioning κ(A) matrices that do not have Property A (as an indication, note that if Aii = 0 for some i, diagonal scaling is ill-defined but BIN is still possible to execute). ˜ An interesting problem would be to prove a bound similar to (18) for BIN’s A, possibly with a constant less than n or p. As it turns out from our numerical ˜ is smaller, sometimes far smaller, than κ(A˜D ). experiments, κ(A) Appendix A contains an illustrative MATLAB code of BIN and numerical examples.

16

BINORMALIZATION

5

BIN Extensions

5.1

Hermitian Matrices

BIN can be immediately generalized to Hermitian matrices by setting Bij := |Aij |2 in §2.1, the rest of the algorithm of §3 remaining intact. Our numerical experiments show similar convergence rates as for the real case.

5.2

Lp normalization

BIN can be used to scale the Lp -norms of the rows of A to 1 by replacing B = A ◦ A by Bij = |Aij |p , i, j = 1, . . . , n. (19) The scalability properties of A and convergence results (§2.3 - §4) carry over, for any 0 < p < ∞.

5.3

Non-Symmetric and Rectangular Matrices

It is possible to extend BIN to the case of a general m × n matrix A. We assume for simplicity that A is real; however, it may be complex-valued throughout this section. The scaling matrix has the form A˜ = F AG where F = diag(f1 , . . . , fm ) and G = diag(g1 , . . . , gn ) are the rows- and columnscaling factors, respectively. Setting Bij := |Aij |2 ,

xi := fi2 ,

yj := gj2 ,

βi :=

n X

j=1

Bij yj ,

γj :=

m X

Bij xi ,

i=1

Here, the analogous form of (5) is

for some b, c > 0, where

B(y)x = be C(x)y = ce,

B(y) := diag(β1 (y), . . . , βm (y)),

C(x) := diag(γ1 (y), . . . , γn (y));

(20) (21)

Note that c may be different from b, otherwise (20)–(21) may be insoluble, e.g., for m = 1, n > 1. Eq. (6) is generalized to   1 T Im − e˜e˜ B(y)x = 0 (22) m   1 In − eeT C(x)y = 0, (23) n

17

BINORMALIZATION

where e˜ := (1, . . . , 1)T . This is an (m + n) × (m + n) system of autonomous |

{z

}

m

quadratic equations in z := (x, y)T , where x = (x1 , . . . , xm ), y = (y1 , . . . , yn ). In block form,  

Φ(z) := 

Im −

1 e˜e˜T m

0



B(y)



0    In − n1 eeT C(x)

x y

!

=

0 0

!

. (24)

Because (24) is a non-symmetric system, we cannot use a GS type iteration. Instead, we can use one of the following ideas, or their combinations: • Kaczmarz relaxation. The linear Kaczmarz relaxation [6, §3],[7] is a distributive relaxation scheme suited for general non-symmetric systems of equations, say Ax = b. At any step i = 1, . . . , m we satisfy equation i by changing all the unknowns appearing in this equation with weights relative to their coefficients in that equation. Namely, a Kaczmarz step means xj ← xj + δAij ,

j = 1, . . . , n,

(25)

P

where δ is chosen so that the new xj s satisfy m j=1 Aij xj = bj . A Kaczmarz step is equivalent to a GS step for the system AAT x˜ = b, thus convergence is always ensured [7]. Analogously, we could define the Kaczmarz-Newton (KN) method for (24) by sweeping over all (m + n) equations and changing all xi and yj at any step. For instance, for the ith equation, i = 1, . . . , m we update 

xi ← x i + δ 1 − xk yj





1 βi m 

1 ← xk − δ βk , m ← yj − δ (Bij xi − γj ) ,

(26) ∀k 6= i ∀j = 1, . . . , n,

(27) (28)

where βi , γj are computed with the values of the xi ,yj prior to the (i.e. B(y) is “frozen” during the Kaczmarz step). δ is chosen to satisfy (22) for that i with the updated xi ’s and yj ’s. Note that the Kaczmarz sweep’s complexity is O(˜ nz ) where n ˜ z is the number of non-zero entries T in AA . Although this complexity can be much higher than for the symmetric case (especially when A is full), it is still feasible in many

BINORMALIZATION

18

cases. For instance, if A is a k-point stencil n × n discretization of some PDE, AAT may typically have n ˜ z ≈ k 2 n compared with nz = kn of A. For the case of m = n we can also try a “red black” ordering in the Kaczmarz sweep, namely, sweep on row 1, then column 1, then row 2, then column 2, etc. Numerical experiments show that the “lexicographic” (sweeping over all rows, then over all columns) and “red black” ordering, give similar convergence factors; see §B.3.

• Autonomous systems; composite relaxation. One can save a significant portion of the work in the Kaczmarz sweep by employing instead Composite Kaczmarz Newton (CKN) relaxation on (22)–(23). The notion of “composite relaxation” was introduced in [8],[19, p. 50],[20],[25, Chap. 4] for relaxing products of differential (or discrete) operators. In our case, a CKN on a current approximation (x(0) , y (0) ) to the solution of (22)–(23) is defined as follows: – Perform ν1 Kaczmarz sweeps on the equation B(y (0) )x = 0 on x only starting from x(0) , obtaining an updated x(1) . – Perform ν2 Kaczmarz sweeps on the equation C(x(1) )y = 0, starting from y (0) and obtaining y (1) . – The new approximation is (x(1) , y (1) ). Note that at each stage we perform Kaczmarz sweeps on a linear system of equations, which is cheaper than for the original non-linear system. Furthermore, at every step we modify fewer unknowns (e.g. only the x’s in the first stage). The convergence factor may be slightly worse than for the full Kaczmarz sweep, but the cost-effectiveness may still be higher [20]. • Non-average form; Jacobi relaxation. Instead of using (22)–(23), we go back to (20)–(21). Now we apply for instance KN or CKN to this system, or separately to (20) for x, and then to (21) for y. The latter form is in fact equivalent to Composite Jacobi (CJ) relaxation, since the corresponding M, N are diagonal. The CJ sweep on a current approximation (x, y) consists of the following steps: 1. Set xi ← 1/βi for all i = 1, . . . , m.

2. Update γ ← B T x.

19

BINORMALIZATION 3. Set yj ← 1/γj for all j = 1, . . . , n.

4. Update β ← By.

Note that this sweep costs twice a BIN sweep when m = n, i.e. 2 MVMs. Working with the non-average form (5) in the symmetric case is also possible, but it leads to a typically slower GSN convergence rate than with form (6), in our numerical experiments. • Newton as an outer iteration. Instead of using the relaxation scheme as the primary (outer) iteration and Newton’s method (or the explicit quadratic formula) as the inner iteration, we can reverse their roles [21]. That is, we solve the system (22) with Newton’s method where a few relaxation sweeps are applied to the solution of −∇Φ(z)δz = −Φ(z), where z is our current approximation and δz is the Newton displacement vector. This approach needs to be further investigated, for both the symmetric and non-symmetric cases. • Secular Equations Zero finders. An alternative approach to solving (6) (similarly (24)) is by casting it in the form x = B −1 (x), for which Newton’s method is applied. This form proves highly beneficial for solving certain “secular equations” for constrained least-squares problems [14]. Although a direct application of this approach to (6) leads to the Newton iteration x(k+1) = x(k) − B −1 B(x(k) )x(k) , which is feasible only when B is invertible, there may be other ways of incorporating this approach, in conjunction with the previous techniques mentioned in this section. 5.3.1

Numerical Experiments

So far, the best non-symmetric scaling algorithm seems to be the CJ variant (steps 1-4 mentioned previously in this section). An advantage of CJ over the other variants is that it automatically ensures convergence to positive x,y, for any positive initial guesses x(0) , y (0) . Jacobi relaxation is also ideal for parallel computing. For many cases, its asymptotic convergence is quite slow (around 0.95) compared with BIN’s asymptotic convergence; however, the condition number can still be reduced using only a few sweeps; these sweeps tend to converge

BINORMALIZATION

20

much faster, as in the symmetric case. The possible reason for this generally fast initial convergence is the ability of relaxation to reduce the residual norm efficiently as long as low-eigenvalue eigenvectors do not dominate the error [6, §1]; by the time these modes overtake convergence, the row and/or sums are already averaged well enough for most practical purposes. An illustrative MATLAB code of the non-symmetric BIN (NBIN) algorithm appears in Appendix B.1. Appendix B.2 contains some numerical examples.

6

Concluding remarks

6.1

Main Findings

In this paper we introduce the BIN algorithm for scaling general real symmetric matrices. The algorithm scales the rows and columns to unit 2-norm, preserving the matrix’s symmetry. BIN is an iterative method based on the Gauss Seidel Newton (GSN) method for solving the quadratic system of “binormalization equations” (6). • The BIN algorithm converges for any scalable matrix A; we can also use it to detect (near-) unscalability. • Although there is no a priori bound on the convergence factor, convergence is fast for a wide class of matrices; this includes many matrices arising from PDE discretizations, and more. • For all practical purposes, only 1 − 3 sweeps of BIN are needed for quite accurate scalability. The algorithm complexity is then 2nz to 4nz operations, where nz is the number of non-zero entries in A. • BIN tends to reduce the condition number of A toward optimality, after a few sweeps. Although diagonal scaling must be optimal when A has Property A, BIN yields nearly optimal κ(A) in that case; moreover, when A does not have Property A, BIN seems to perform better in many cases. Its performance seems to be more robust than other scaling algorithms (e.g. diagonal scaling). • BIN’s main drawback is its higher complexity with respect to “trivial” scaling algorithms (diagonal or row scaling): all scale linearly with nz ,

21

BINORMALIZATION

but BIN’s complexity has larger constant. This is especially revealed in cases where its asymptotic convergence is poor and absolutely accurate scaling is desired, although the first few sweeps of interest still converge much faster. Moreover, simultaneous relaxation [6, §3] may help restore fast asymptotic convergence, in cases where only a few individual rows overtake (slow) convergence. At any rate, BIN may be of importance for various matrix scaling problems, especially when once-and-for-all conditioning or approximate scaling of the matrix are desired. • BIN can be extended to Hermitian matrix and to non-symmetric and rectangular matrices. Our current variant of the non-symmetric BIN (NBIN) is based on Jacobi relaxation, each sweep costs twice as a BIN sweep, and asymptotic convergence is much slower; yet, the condition number can still be reduced using only a few, faster-to-converge sweeps.

6.2

Open Problems

Further research may be focused in several directions. First, BIN is only one possible algorithm for solving the binormalization equations (2); other approaches are feasible and worth exploring; for example, Newton-GS type methods (see §5.3). Focusing on BIN, it may still be possible to bound the convergence factor for a certain class of matrices, even though it is not generally separated from one. Moreover, the question of BIN’s relation to optimal conditioning remains open: does there exist a factor α such that ˜ ≤α κ(A)

D

min κ(DAD), diagonal

(29)

where A˜ is the BIN-scaled matrix obtained from a real symmetric A? Such a result would correspond to van der Sluis’ Theorem [22] for diagonal scaling and α = n. It may be possible that the 2-norm-based-BIN’s α is less than n as it experimentally exhibits a more robust behavior than diagonal scaling. Further numerical and theoretical studies regarding the impact of BIN on the fill-in in sparse factorizations, are of highly practical importance. Finally, the non-symmetric and rectangular variant NBIN is still not fully developed, although it can be already used for practical numerical experiments. Convergence theory and practical use of NBIN (§5) need to be formulated and experimented in detail. In particular, an inequality like (29)

BINORMALIZATION

22

may hold for NBIN as well, as van der Sluis’ theorem [22] can be extended to non-symmetric and rectangular matrices. It may be also interesting to explore the relation of the result of NBIN to the best l2 -scaling of a matrix (see [18]). Again, other approaches for solving (22–23) exist. Finally, even if a matrix is unscalable (see the examples in §2.3), the matrix may be “ε-scalable”, namely, there may be a solution to (2) with an ε-norm residual. Conditions characterizing the set of ε-scalable matrices and estimates of the condition number, would be valuable.

A

Numerical Experiments: BIN

We tested BIN over several representative model problems, to check various aspects discussed in §4. For each of the following model problems, we examined the convergence properties of BIN and the effect of scaling on the condition number κ(A) [17]. §A.3 contains the results of BIN over the University of Florida Sparse Matrix Collection [3]. Moreover, it turned out that BIN improved the efficiency of the sparse LU factorization for symmetric matrices, which is part of the UMFPACK V4.0 package [2].

A.1

BIN MATLAB Code

function [C,f,avg_c,est_c] = bin(A,max_iter,verbose,x) % BIN.M - the BIN binormalization algorithm for a real % symmetric matrix A, based on the average % binormalization equations % Input : * NxN matrix A % * maximum number of sweeps (default = 10) % * verbose mode flag (0=no printout, default; % 1=printouts) % * x = (x1,...,xn),initial condition for fi^2. % Default: all xi = 1 % * plots = plots flag (0=false,Default; 1=true) % Output : * The normalized matrix C = F*A*F % * The scaling factors f = (f1...fN); F = diag(f) in % the previous line. % * avg_c = average convergence factor of the

23

BINORMALIZATION % % %

residuals of the binormalization equations. * est_c = estimated convergence factor (see Theorem 5).

% Define parameters n = size(A,1); if (nargin < 2) max_iter = 10; end if (nargin < 3) verbose = 0; end if (nargin < 4) x = ones(n,1); % initial guess; x(i) = f(i)^2. end if (nargin < 5) plots = 0; end TOL = 1e-10; % stopping threshold % Initializations B = d = beta = avg = e = std_initial =

abs(A).^2; diag(B); B*x; x’*beta/n; 1.; full(sqrt(sum((x.*beta-avg).^2)/n))/avg;

if (verbose) fprintf(’\nBINORMALIZTION - BIN Algorithm\n\n’); fprintf(’Sweep STD RATE DIFF-X RATE\n’); fprintf(’INITIAL %.3e\n’,std_initial); end % The main loop over BIN sweeps for r = 1:max_iter x_old = x;

BINORMALIZATION

24

% Reached Below tolerance ==> break if (std < TOL) break; end % step over all variables i=1,...,n and update each by turn % (Gauss Seidel) for i = 1:n, % Solve quadratic equation for the updated xi bi = beta(i); di = d(i); xi = x(i); c2 = (n-1)*di; c1 = (n-2)*(bi-di*xi); c0 = -di*xi^2 + 2*bi*xi - n*avg; if (-c0 < eps) fprintf(’Nearly unscalable: row %d, c0=%f\n’,i,c0); C = A; f = ones(n,1); avg_c = -1; G = -1; est_c = -1; return; else xi = (2*c0)/(-c1 - sqrt(c1*c1 - 4*c2*c0)); end % Update the auxiliary variables accordingly delta = xi - x(i); % xi_new - xi_old avg = avg + (delta*x’*B(:,i) + ... delta*bi + di*delta^2)/n; beta = beta + delta*B(:,i); x(i) = xi; end % Compute std_old = e_old = std =

and print the convergence statistics std; e; full(sqrt(sum((x.*beta-avg).^2)/n))/avg;

25

BINORMALIZATION e = norm(x-x_old)/norm(x); conv_factor = std/std_old; fprintf(’%3d %.3e %.3f %.3e r,std,conv_factor,e,e/e_old); end

%.3f\n’,...

end % Print convergence factor and display err if desired avg_c = (std/std_initial)^(1/r); if (verbose) fprintf(’\nAverage convergence factor : %.5f\n’,avg_c); end % Finalize f = F = C = beta = avg = norm_c = C = f =

by computing f x and scaling everything to 2-norm=1 sqrt(x); spdiags(f,0,n,n); F*A*F; sum(C.^2)’; full(sum(beta)/n); sqrt(avg); C*(1./norm_c); f*(1./sqrt(norm_c));

% Compute an estimate of the asymptotic convergence factor if (nargout >= 4) w = ones(n,1); G = eye(n) + C.^2 - 2*w*w’/n; H = -inv(tril(G,0))*triu(G,1); lambda = flipud(sort(abs(eig(H)))); est_c = lambda(2); fprintf(’Estimated asymp. convergence rate: %.5f\n’,... est_c); end if (verbose > 0) fprintf(’Original condition # : %f\n’,condest(A)); fprintf(’After scaling condition # : %f\n\n’,condest(C)); end

26

BINORMALIZATION

A.2 A.2.1

Numerical Experiments n=2

Let

1 1 1 100

A=

!

.

(30)

Table 1 shows the convergence history of BIN. We define the residual for a given approximation x to (6) by ri = x i

X j

Bij xj −

1X X xk Bkj xj , n k j

i = 1, . . . , n;

Let x(k) be the approximation after k iterations and r (k) the correspond(k) ing residual. The table contains the values of R(k) := √1n kr (k) k2 , RR := R(k) /R(k−1) , E (k) := sweep.

√1 kx(k) n

(k)

− x(k−1) k2 and RE := E (k) /E (k−1) , after each

Table 1: Convergence history of BIN for (30) #Sweep R(k) Initial 9.996 ∗ 10−1 1 0.000 ∗ 100 κ(A) = 1.0103 ∗ 102 ˜ = 1.2222 ∗ 100 κ(A)

(k)

RR − 0.000

E (k) − 9.900 ∗ 10−1

(k)

RE − −

Note that we reach exact scaling with one BIN sweep; this is fact true for all 2 × 2 scalable matrices. Moreover, the original large condition number reduced almost to 1. In addition, A˜ =

0.99504 0.099504 0.099504 0.99504

!

= 0.99504

1 0.1 0.1 1

!

= 0.99504A˜D

where A˜D is the diagonal scaling (see §4.3). Hence BIN’s A˜ is optimally conditioned.

27

BINORMALIZATION A.2.2

Perfectly Conditioned, Property A

The discrete Laplacian matrix 

2 −1  −1 2 −1  A=  ... ... ... −1 2

    

(31)

is perfectly conditioned by Forsythe–Straus Theorem [13, Theorem 4], since it satisfies property A and has all diagonal entries equal. Running BIN here can only increase κ(A), but it turns out that this increase is small, as shown in Table 2. Table 2: Convergence history of BIN for (31), n = 50 (k)

#Sweep R(k) RR E (k) INITIAL 1.788 ∗ 10−2 − − −3 1 1.028 ∗ 10 0.057 1.017 ∗ 10−2 2 1.478 ∗ 10−4 0.144 5.467 ∗ 10−4 3 1.258 ∗ 10−5 0.085 7.831 ∗ 10−5 4 1.765 ∗ 10−6 0.140 6.627 ∗ 10−6 5 1.869 ∗ 10−7 0.106 9.300 ∗ 10−7 6 2.456 ∗ 10−8 0.131 9.853 ∗ 10−8 7 2.511 ∗ 10−9 0.102 1.294 ∗ 10−8 8 3.362 ∗ 10−10 0.134 1.323 ∗ 10−9 9 3.199 ∗ 10−11 0.095 1.774 ∗ 10−10 κ(A) = 4.977325 ∗ 100 ˜ = 4.977346 ∗ 100 κ(A) Average convergence factor: 0.133 Estimated convergence factor: 0.114

(k)

RE − − 0.054 0.143 0.085 0.140 0.106 0.131 0.102 0.134

The estimated convergence factor (σ) was computed using Theorem 5. 1 The average convergence factor is defined by (kR (9) k2 /kR(0) k2 ) 9 . The two numbers are indeed close.

28

BINORMALIZATION A.2.3

Perfectly Scaled

Another check of our BIN code is to start with a scaled matrix and observe that BIN does not alter the matrix. Of course, if we start from xi = 1 for all i this readily happens, since (6) is already satisfied. Thus, we started from 1 random 0 ≤ xi ≤ 1, and run BIN for the matrix C := A 2 where A is given by (31). Note that C is full and does not have Property A. Since A = C 2 = CC T has all diagonal entries equal, by Theorem 6 (i) it follows that C is scaled. Indeed, BIN converged to C˜ = C with an asymptotic convergence rate of 0.104, starting from a random initial guess. A.2.4

Badly Conditioned, Property A

In contrast to (31), we consider the following badly scaled n×n sparse matrix 

2 −1   −1 2 −1   ... ... ...   −1 2 −1 A=   −1 1002 −1   ... ... ...  −1 1002

            

.

(32)

n×n

That is, A is tridiagonal with −1 on the upper and lower sub-diagonals; on the diagonal we have n/2 entries equal to 2 and n/2 entries equal to 1002. A corresponds to a discretization of the 1D differential operator −uxx + au, where a = a(x) is discontinuous (a = 1 in one half of the domain and a = 1000 in the other half). Clearly, A has Property A. Table 3 shows the convergence history of BIN for (32), n = 50. Note that A has Property A; thus, diagonal scaling provides the optimal condition number, κ(A˜D ) = 273.327, versus the original κ(A) ≈ 6.886 × 104 . ˜ = 273.349; interestingly, κ(A) ˜ − κ(A˜D ) = Alternatively, BIN yields κ(A) −2 2.132 × 10 . In other words, BIN yielded a nearly optimal scaling for this problem, although the diagonal entries of A˜ are not equal: they range between 0.81 and 1 (see Fig. 1). Moreover, A˜n/2,n/2+1 = A˜n/2+1,n/2 = −2.124 × 10−2 , correctly indicating that the matrix is nearly reducible to two separate tridiagonal matrices. Finally, BIN’s convergence factor does not deteriorate when n gets larger. For instance, for n = 500 it turns out to be 0.098.

29

BINORMALIZATION

Table 3: Convergence history of BIN for (32), n = 50 (k)

#Sweep (k) R(k) RR E (k) INITIAL 1.000 × 100 − − −1 1 1.179 × 10 0.118 9.973 × 10−1 −3 2 8.413 × 10 0.071 9.232 × 10−2 3 1.322 × 10−3 0.157 6.388 × 10−3 4 1.013 × 10−4 0.077 1.021 × 10−3 5 1.076 × 10−5 0.106 7.688 × 10−5 6 9.017 × 10−7 0.084 8.286 × 10−6 7 6.654 × 10−8 0.074 6.840 × 10−7 8 6.646 × 10−9 0.100 4.830 × 10−8 9 9.669 × 10−10 0.145 4.828 × 10−9 10 7.311 × 10−11 0.076 7.246 × 10−10 κ(A) = 6.88551 × 104 ˜ = 2.73327 × 102 κ(A) κ(AD ) = 2.73348 × 102 Average convergence factor: 0.097 Estimated convergence factor: 0.104

1

ii

A~

0.95

0.9

0.85

0.8 5

10

15

20

25 i

30

35

40

45

(k)

RE − − 0.093 0.069 0.160 0.075 0.108 0.083 0.071 0.100 0.150

50

Figure 1: The diagonal entries A˜ii versus i for (32).

30

BINORMALIZATION A.2.5

Semi Positive Definite

In many PDE problems we encounter semi positive definite matrices. This often happens with periodic boundary conditions. For example, consider the boundary value problem −∆u = 0,

(x, y) ∈ [0, 1] × [0, 1],

(33)

where u = u(x, y) is 1-periodic in x and y. (33) has non-zero solutions, u ≡ c. Discretizing (33) on a uniform grid with meshsize h yields the following discrete equations: 1 (4uij − ui−1,j − ui+1,j − ui,j−1 − ui,j+1 ) = 0, i, j = 1, . . . , ng , (34) h2 where h = 1/n, with the boundary conditions ui+n,j = uij and ui,j+n = uij for any i, j. Ordering the variables lexicographically, the resulting matrix is semi positive definite and has the sparsity pattern shown in Figure 2. The matrix is n × n, n = n2g and has 5n2 non-zero entries. 0

10

20

30

40

50

60

0

10

20

30 nz = 320

40

50

60

Figure 2: The sparsity pattern of the matrix formed from (34) for ng = 8. (34) provides a non positive definite realistic test case for BIN. Since it is already scaled, we start BIN from random values for xi for this matrix, which is in fact already scaled. The average convergence factor over ten sweeps is given in Table 4, for various n.

31

BINORMALIZATION

Table 4: Convergence factors of BIN for the discretized periodic Laplacian matrix. nz is the number of non-zero entries in A. ng n nz µ 8 64 320 0.089 16 256 1280 0.093 32 1024 5120 0.096 64 4096 20480 0.092 For ng = 64: κ(A) = 9.8956 ∗ 102 ˜ = 4.91392 ∗ 100 κ(A) Average convergence factor: 0.110 Estimated convergence factor: 0.113

A.2.6

Random Matrices

We consider an n × n symmetric matrix A whose entries are uniformly distributed random numbers between −1 and 1. Table 5 summarizes the convergence of BIN for such matrices. For each n we generated 100 matrices and averaged BIN’s average convergence factor over all of them to get the value µ listed in the second column. µ0 refers to the same quantity, over 100 random matrices with zero diagonal. Table 5: Convergence factors of BIN for random symmetric matrices. n 8 16 32 64 128 256 512 1024

µ 0.187 0.137 0.117 0.113 0.114 0.122 0.118 0.114

µ0 0.264 0.150 0.118 0.111 0.115 0.121 0.119 0.121

32

BINORMALIZATION

It follows that BIN’s convergence does not depend on the size of the matrix, but on its specific “nature” (i.e. sparsity pattern and entries). See also §A.2.7 A.2.7

3 × 3, Slow To Converge

In view of the results so far, one might convict that BIN converges fast for any matrix. However, it is possible to construct a matrix with BIN convergence factor as close to 1 as desired. Consider Let 



0 1 1  A= 1 0 1   1 1 a

(35)

for some positive a. The binormalization equations (6) have the form x1 x2 + x 1 x3 = β x1 x2 + x 2 x3 = β x1 x3 + x2 x3 + a2 x23 = β,

(36)

where β = 2(x1 x2 + x1 x3 + x2 x3 ) + a2 x23 . The unique positive solution of (36) is x1 = x2 = 0.5(1 + d)x3 x3 =

1 + a2 − d a4 − 2a2

!1 2

,

1 ˜ G (see (16)) where d := (1 + 4a2 ) 2 . Thus, we can analytically compute A, and its eigenvalues. It turns out that G has the three eigenvalues

λ1 = 2 − 3a−1 + 4.5a−2 + O(a−3 ), λ2 = 0.5a−1 − 0.75a−2 + O(a−3 ), λ3 = 0. The convergence of Gauss-Seidel relaxation depends on the value of λ2 (λ3 should be discarded as it corresponds to the degree of freedom in scaling xi by any constant). Since λ2 = O(a−1 ) → 0 as a → ∞, the convergence factor is O(1 − a−1 ) can be made arbitrarily close to 1.

33

BINORMALIZATION

Table 6: Convergence history of BIN for (A.2.7), a = 104 . The last column shows the condition number of A˜ after each sweep. #Sweep INITIAL 1 2 3 4 5 10 1000

R(k) 1.421 ∗ 100 1.137 ∗ 10−1 7.104 ∗ 10−2 5.147 ∗ 10−2 4.029 ∗ 10−2 3.308 ∗ 10−2 1.740 ∗ 10−2 1.755 ∗ 10−4

(k)

RR − 0.080 0.625 0.725 0.783 0.821 0.905 0.999

E (k) − 1.000 ∗ 100 5.793 ∗ 10−1 3.160 ∗ 10−1 2.164 ∗ 10−1 1.643 ∗ 10−1 7.415 ∗ 10−2 6.480 ∗ 10−4

(k)

RE − − 0.579 0.546 0.685 0.759 0.890 0.999

˜ κ(A) 104 3.099 2.384 2.077 1.901 1.785 1.517 1.052

On the other hand, running BIN in practice (see Table 6) we observe that the first few sweeps convergence much faster; A˜ is already well conditioned after these few sweeps. . Consequently, the scaled matrix is also well conditioned, by using only 1 − 2 sweeps of BIN during which convergence is fast as well. A.2.8

Toeplitz Matrices

An interesting example is the Toeplitz matrix A = (Aij )ni,j=1 ,

Aij = ρ|i−j|

(37)

that describes a first order Markov chain [16, §4]. A does not have Property A but A−1 does since it is tridiagonal [16]. The optimal scaling turns out to be A˜G := DG ADG , DG = diag(c, 1, . . . , 1, c), (38) 1

where c := (1 + ρ2 )− 2 . Our numerical experiments were performed with n = 50. Interestingly, BIN’s scaling factors behave the opposite way: more weight is given to the extreme rows (see Figure 3) which may seem a more “natural” choice from a physical point of view (giving more weight to the initial condition and the last few observations, which may be more accurate). However, BIN’s scaling has a slightly larger condition number than the original one, where (38) has a smaller, but not considerably, condition number (see Figure 4).

34

BINORMALIZATION

6

10

5.6

10

5

10

5.5

10 4

10

3

10 0.965

0.97

0.975

0.98

0.985

0.99

0.995

1

1.005

0.9996

0.9996

0.9997

0.9997

0.9997

0.9997

0.9997

0.9998

0.9998

Figure 3: Left figure: the condition number of the original (37) versus the condition number after BIN (circles) and the optimal conditioning ((38), x-marks), versus ρ. Right figure: a zoom-in of the left figure near ρ = 1.

1.1

1.05

1

0.95

0.9

0.85

0.8

0.75

0.7

0

5

10

15

20

25

30

35

40

45

50

Figure 4: Scaling factors of BIN (solid line) and optimal scaling factors (dashed) versus row number, for ρ = 0.99.

BINORMALIZATION

A.3

35

The UF Sparse Matrix Collection

The UF Sparse Matrix Collection [2, 3] is maintained by Prof. Tim Davis and includes about a thousand sparse matrices representing real-life applications as well as esoteric matrices with bad numerical properties. We compared BIN versus existing scaling algorithms: • Diagonal scaling, symmetrically scales rows and columns to obtain all diagonal entries equal to one [13]. • Row and column scaling: first, sweep over the rows and scale by each the maximal element in each row; then sweep over columns and scale by the maximal element in each column. • MC29: a routine developed by Reid and Curtis [11], [24, p. 247]. • MC29 followed by row and column scaling, as recommended by the NAG Harwell Sparse Matrix Library Release 2 [1]. All of these four algorithms were tested on 209 real symmetric test matrices from the collection, and for each matrix we show the best-performing algorithm (denoted “BEST”). BIN was run for the 209 matrices as well. Figure 5 (left) shows the original condition number and the condition number after scaling with BEST and with 3 sweeps of BIN. Figure 5 (right) shows the reduction in condition number for BEST versus 3 sweeps of BIN, and Figure 6 shows the average convergence factor of BIN (over 10 sweeps) versus the matrix size. As can be seen, BIN performed as well as the other scaling algorithms and much better in many cases. Moreover, its convergence factor does not depend on the size of the problem and ranges between 0 and 0.55.

B B.1

Numerical Experiments: NBIN NBIN MATLAB Code

function [C,f,g] = nbin(A,max_iter,verbose) % NBIN.M - the NBIN binormalization algorithm % for a general matrix A % Input : * mxn matrix A

36

BINORMALIZATION

40

10

0

10

35

10

−5

30

10

10

25

10

−10

10

20

10

−15

10 15

10

−20

10 10

10

−25

10

5

10

0

10

−30

0

20

40

60

80

100 120 Test Case Number

140

160

180

10

200

0

20

40

60

80

100 120 Test Case Number

140

160

180

200

Figure 5: Left figure: the original condition number of each of the 209 test cases (black), the condition number κ1 after scaling with BEST (green), and the condition number κ2 after 3 sweeps of BIN (red). Right figure:κ2 /κ1 for each test case.

0.6

0.5

0.4

0.3

0.2

0.1

0

2

3

10

10

4

10

n

Figure 6: Average convergence factor of BIN for the 209 test cases versus the matrix size (n).

37

BINORMALIZATION % * maximum number of sweeps (default = 10) % * verbose mode flag (0=no printout, % default; 1=printouts) % Output : * The normalized matrix C = F*A*F % * The row-scaling factors f = (f1...fm), % column scaling factors g = (g1,...,gn). TOL if (nargin < 2) max_iter end if (nargin < 3) verbose end [m,n] x y B C beta gamma sum1 sum2 std1 std2 std_initial

= 1e-10; = 10;

= 0;

= = = = = = = = = = = =

size(A); ones(m,1); % initial guess; x(i) = f(i)^2. ones(n,1); % initial guess; y(i) = g(i)^2. abs(A).^2; B’; B*y; C*x; n; m; full(sqrt(sum((x.*beta-sum1).^2)/m))/sum1; full(sqrt(sum((y.*gamma-sum2).^2)/n))/sum2; sqrt(std1^2+std2^2);

if (verbose) fprintf(’\nBINORMALIZTION - NBIN Algorithm\n\n’); fprintf(’Sweep STD RATE\n’); fprintf(’INITIAL %.3e\n’,std_initial); end % The main loop over NBIN sweeps for r = 1:max_iter x_old = x; y_old = y;

BINORMALIZATION

38

% Below tolerance or convergence stalls due to round-off if (std < TOL) break; end x = sum1./beta; x(find(abs(beta) < eps)) = 1.0; gamma = C*x; y = sum2./gamma; y(find(abs(gamma) < eps)) = 1.0; beta = B*y; %%%%%% Compute convergence factors std_old = std; std = full(sqrt(sum((x.*beta-sum1).^2)/m))/sum1; conv_factor = std/std_old; if (verbose) fprintf(’%3d %.3e %.3f\n’,... r,std,conv_factor); end end if (verbose) fprintf(’\nAverage residual convergence rate: %.5f\n\n’,... (std/std_initial)^(1/r)); end f g F G C

B.2

= = = = =

sqrt(abs (x)); sqrt(abs (y)); spdiags(f,0,m,m); spdiags(g,0,n,n); F*A*G;

Numerical Experiments: NBIN

As a first example, we consider the matrix no. 632 in the UF Sparse Matrices Collection (unsymmetric matrix submitted by Francois Cachard, Grenoble,

39

BINORMALIZATION

March 1981). A is a 185 × 185 real non-symmetric matrix with sparsity pattern shown in Figure 7. 0

20

40

60

80

100

120

140

160

180 0

20

40

60

80 100 nz = 975

120

140

160

180

Figure 7: The sparsity pattern of the Grenoble matrix. Table 7 shows the convergence history of NBIN. We define the residual r for a given approximation x, y to (20) by ri = x i

X j

Bij yj − n,

i = 1, . . . , n.

(39)

Let x(k) be the approximation after k iterations and r (k) the corresponding (k) residual. The table contains the values of R(k) := √1n kr (k) k2 and RR := R(k) /R(k−1) after each sweep (note that the residual of (21) is zero at the end of each Jacobi sweep). Note that the condition number was reduced by a substantial factor, as for the symmetric case.

B.3

Different Ordering of Rows and Columns

If the Kaczmarz version of NBIN is used (see §5), we can order the rows and columns differently within the Kaczmarz sweep. For instance, sweep on row 1, then column 1, then row 2, then column 2, etc. This ordering is denoted “red black” ordering, versus the original “lexicographic” ordering

40

BINORMALIZATION

Table 7: Convergence history of NBIN for the Grenoble matrix #Sweep (k) R(k) Initial 9.981 × 10−1 1 1.708 × 10−1 2 1.480 × 10−1 3 1.338 × 10−1 4 1.230 × 10−1 5 1.141 × 10−1 6 1.065 × 10−1 7 9.984 × 10−2 8 9.394 × 10−2 9 8.862 × 10−2 10 8.378 × 10−2 κ(A) = 5.0878 × 105 ˜ = 4.2926 × 104 κ(A)

(k)

RR − 0.171 0.866 0.904 0.919 0.928 0.933 0.938 0.941 0.943 0.945

(sweeping over all rows, then over all columns in an NBIN sweep). This is only meaningful for m = n, but may be generalized to arbitrary m, n. It turns out that the two ordering methods yield similar convergence factors for NBIN. For different values of n, we computed the convergence factor with both methods, averaged over 10 matrices whose entries were random in [−1, 1]. The results are summarized in Table 8; note that the asymptotic convergence factor becomes better with n, but tends to a positive limit (about 0.1).

B.4

Rectangular Matrices

As a second example, we consider the following rectangular 3 × 6 matrix 



0.3891 0.9137 −0.6541 −0.4953 −0.7270 −0.6017  0.9595 0.7515 −0.9765 −0.4026  A =  0.2426 0.0452 . 0.5896 0.7603 −0.4571 0.4746 0.7878 0.3229 The convergence history of NBIN is summarized in Table 9.

(40)

41

BINORMALIZATION

Table 8: Convergence Factors of Lexicographic versus Red Black Ordering n 8 16 32 64 128

Lexicographic 0.6433 0.4362 0.2349 0.1374 0.1172

Red Black 0.6378 0.4542 0.2633 0.1498 0.1165

Table 9: Convergence history of NBIN for (40) #Sweep (k) R(k) Initial 8.610 × 10−1 1 6.402 × 10−2 2 1.424 × 10−2 3 3.284 × 10−3 4 7.624 × 10−4 5 1.773 × 10−4 6 4.124 × 10−5 7 9.593 × 10−6 8 2.232 × 10−6 9 5.192 × 10−7 10 1.208 × 10−7 κ(A) = 1.4829 ˜ = 1.2258 κ(A)

(k)

RR − 0.074 0.222 0.231 0.232 0.233 0.233 0.233 0.233 0.233 0.233

42

BINORMALIZATION The scaled matrix has the form 



0.8429 1.2765 −0.8615 −0.7961 −0.8186 −1.2705  ˜ 1.3487 1.2891 −1.1735 −0.9072  A =  0.5609 0.0674 , 1.4053 1.1688 −0.6625 0.8394 0.9761 0.7502

with the corresponding scaling factors

f = (1.4945, 1.5950, 1.6445)T , g = (1.4493, 0.9348, 0.8813, 1.0754, 0.7535, 1.4129)T . All row sums of squares of A˜ are 6 where the column sums are 3. In this example even the asymptotic convergence is fast. Finally, we tested NBIN on the non-symmetric and rectangular matrices in the UF Collection (the results are not shown here). NBIN converged for all matrices and reduced the condition number of the non-symmetric matrices by a substantial factor (1.5 to 1030 ). Interestingly, NBIN can sometimes reduce the condition number of a symmetric matrix, although the scaled matrix is no longer symmetric. In sum, NBIN is a useful tool for scaling and conditioning non-symmetric as well as rectangular matrices, although not as efficient as BIN for symmetric matrices. Finally, note that the matrices may have complex entries in all discussions of this section.

C

Acknowledgments

The first author is grateful to Prof. Achi Brandt from the Weizmann Institute of Science for his innovative ideas and suggestions regarding iterative derivation of diagonal scaling; and to Prof. Tim Davis from the University of Florida for providing the UF Sparse Matrix Collection and for many fruitful discussions. This work was supported by DOE grant DE-FG03-94ER25219.

References [1] http://www.nag.co.uk/ [2] http://www.cise.ufl.edu/research/sparse/umfpack/ [3] http://www.cise.ufl.edu/research/sparse/matrices/

BINORMALIZATION

43

[4] Bauer, F. L., Optimally scaled matrices, Num. Math., 5 (1963), pp. 73– 87. [5] Bellman, R. B., Introduction to Matrix Analysis, SIAM, Philadelphia, 1997. [6] Brandt, A., Guide to multigrid development, in Hackbusch, W. and Trottenberg, U. (Eds.), Multigrid Methods, Proceedings of Conference Held at K¨oln-Porz, November 23–27, 1981, Lecture Notes in Mathematics, Vol. 960, Springer-Verlag, Heidelberg, 1982. [7] Brandt, A., Algebraic multigrid theory: the symmetric case, Appl. Math. Comput., 19 (1986), pp. 23–56. [8] Brandt, A., Barriers to achieving textbook multigrid efficiency in CFD, ICASE Interim Report No. 32, NASA/CR-198-207647. Appears as Appendix C in [25]. [9] Brandt, A., Multiscale scientific computation: review 2001, in Barth, T.J., Chan, T.F. and Haimes, R. (Eds.), Multiscale and Multiresolution Methods: Theory and Applications, Lecture Notes in Computational Science and Engineering, Vol. 20, Springer-Verlag, Heidelberg, 2001. [10] Brandt, A. and Ron, D., Multigrid solvers and multilevel optimization strategies. In Multilevel Optimization and VLSICAD (J. Cong and J. R. Shinnerl, eds.), Kluwer Academic Publishers, Boston, 2002, pp. 1– 69 (in preparation). [11] Curtis, A. R. and Reid, J. K., On the automatic scaling of matrices for Gaussian elimination, J. Inst. Maths. Applics., 10 (1972) 118–124. [12] Eisenstat, S. C., Lewis, J. W., and Schultz, M. H., Optimal block diagonal scaling of block 2-cyclic matrices. Linear Algebra Appl., 44 (1982), pp. 181–186. [13] Forsythe, G. E. and Straus, E. G., On best conditioned matrices, Proc. Amer. Math. Soc., 6 (1955), pp. 340–345. [14] Gander, W., Golub, G. H. and von Matt, Urs, A constrained eigenvalue problem, Linear Algebra Appl., 114/115 (1989), pp. 815–839.

BINORMALIZATION

44

[15] Goldberg, D., What every computer scientist show know about floating-point arithmetic, ACM Comp. Surv., 23 (1991), pp. 5–48. [16] Golub, G. H., Comparison of the variance of minimum variance and weighted least squares regression coefficients, Ann. Math. Stat., 34 (1963), pp. 984–991. [17] Golub, G. H. and Van Loan, C., Matrix Computations, 3rd edition, Johns Hopkins, Baltimore, 1996. [18] Golub, G. H. and Varah, J. M., On the characterization of the best l2 -scaling of a matrix, SIAM J. Numer. Anal., 11 (1974) pp. 472–479. [19] Livne, O. E., Multiscale Eigenbasis Algorithms, Ph. D. Thesis, Weizmann Institute of Science, Rehovot, 2000. [20] Livne, O. E. and Brandt, A., Local mode analysis of multicolor and composite relaxation schemes, Comp. Math. Appl. (in press). [21] Ortega, J. M. and Rheinboldt, W. C., Iterative Solution of Nonlinear Equations in Several Variables, SIAM, Philadelphia, 2000. [22] van der Sluis, A., Condition numbers and equilibration of matrices, Num. Math., 14 (1969), pp. 14–23. [23] Stewart, G. W., Matrix Algorithms, Volume I: Basic Decompositions, SIAM, Philadelphia, 1998. [24] Stewart, W. J., Introduction to the Numerical Solution of Markov Chains, Princeton University Press, Princeton, 1994. ¨ller, A., Multi[25] Trottenberg, U., Oosterlee, C.W., and Schu grid , Academic Press, London, 2000. [26] Varga, R. S., Matrix Iterative Analysis, Springer Verlag, Berlin, 2000. [27] Wilkinson, J. H., The Algebraic Eigenvalue Problem, Oxford University Press, Oxford, 1965.