A Polyakov
email: andrey.polyakov@inria.fr

J I Chairez Oria

A New Homogeneous Quasi-Continuous Second Order Sliding Mode Control

A new homogeneous quasi-continuous second order sliding mode control algorithm is developed. The Lyapunov function for finite-time stability analysis is constructed. The homogeneity property of the algorithm is analyzed. The chattering attenuation procedure using Lyapunov function of the closed-loop system is discussed. The scheme for control parameters tuning based on linear matrix inequalities is elaborated. The theoretical results are supported by numerical simulations.

INTRODUCTION

Introduced 50 years ago sliding mode is the oldest robust control methodology (see, for example, [START_REF] Utkin | Variable structure systems with sliding modes[END_REF] and references therein). The main advantage of sliding mode technique is the insensitivity of the closed-loop system to the so-called matched disturbances and uncertainties, see [START_REF] Edwards | Sliding Mode Control: Theory And Applications[END_REF], [START_REF] Utkin | Sliding Mode Control in Electro-Mechanical Systems[END_REF], [START_REF] Shtessel | Sliding Mode Control and Observation[END_REF].

The concept of the second order sliding mode algorithm was introduced in [START_REF] Levant | Sliding order and sliding accuracy in sliding mode control[END_REF] in order to develop an novel approach to sliding mode control design for systems with relative degree 2 and to reduce of chattering effect (see, [START_REF] Utkin | Sliding Mode Control in Electro-Mechanical Systems[END_REF]) of sliding mode systems. The so-called quasi-continuous sliding mode algorithm (see, Levant [2005a]) looks the most perspective in this context, since it is discontinuous only on the surface {x :∈ R n : s(x) = ṡ(x) = 0}, where s(x) = 0 is a conventional sliding surface, s : R n → R. However, in contrast to other second order sliding algorithms the Lyapunov function-based stability analysis of quasi-continuous sliding mode algorithms have never been presented before. The Lyapunov functions for "twisting", "super-twisting" and "nested" second order sliding mode algorithms are designed in [START_REF] Polyakov | Lyapunov function design for finite-time convergence analysis: "twisting" controller for second-order sliding mode realization[END_REF], [START_REF] Moreno | Strict lyapunov functions for the super-twisting algorithm[END_REF], [START_REF] Polyakov | Unified Lyapunov function for a finite-time stability analysis of relay second-order sliding mode control systems[END_REF]). This paper elaborates the new quasi-continuous second order sliding mode control algorithm. The main advantage of this algorithm is the simplicity of the tuning of control parameters based on the linear matrix inequalities. The stability of the closed-loop system is analyzed using non-smooth Lyapunov function method (see, for example, [START_REF] Polyakov | Stability notions and lyapunov functions for sliding mode control systems[END_REF]). The finite-time convergence of system trajectories to the sliding set is proven by means of implicit Lyapunov function (ILF) technique (see, [START_REF] Polyakov | Finite-time stabilization using implicit lyapunov function technique[END_REF]. The explicit analytical representation of the Lyapunov function is also derived in order to present explicitly the feedback law. The control algorithm and the Lyapunov function are designed using the framework of homogeneous vector fields [START_REF] Zubov | On systems of ordinary differential equations with generalized homogenous right-hand sides[END_REF], [START_REF] Hermes | Nilpotent approximations of control systems and distributions[END_REF] et al), which is rather popular in the sliding mode control theory( see, for example, [START_REF] Orlov | Finite time stability and robust control synthesis of uncertain switched systems[END_REF], Levant [2005a], [START_REF] Bernuau | On homogeneity and its application in sliding mode control[END_REF]).

NOTATION

Through the paper the following notations will be used:

• R + = {x ∈ R : x > 0}, R -= {x ∈ R : x < 0}, where R is the set of real number; • • is the Euclidian norm in R n , i.e. x = x 2 1 + . . . + x 2 n for x = (x 1 , . . . , x n ) T ∈ R n ; • diag{λ i } n
i=1 is the diagonal matrix with the elements λ i on the main diagonal;

• for a matrix P ∈ R n×n , which has the real spectrum, the minimal and maximal eigenvalues are denoted by λ min (P ) and λ max (P ), respectively; • if P ∈ R n×n then the inequality P > 0 (P ≥ 0, P < 0, P ≤ 0) means that P is symmetric and positive definite (positive semidefinite, negative definite, negative semidefinite).

• The set of all subsets of a set M ⊆ R n is denoted by 2 M . • ∇V (x), x ∈ R n is the gradient of the continuously differentiable function V : R n → R.

PROBLEM STATEMENT

The paper deals with the standard problem of the second order sliding mode control design. Let the single input control system be defined as follows:

ż(t) = f (t, z) + g(t, z)u(t), t ∈ R + , (1)
where z ∈ R n is the state vector, u ∈ R is the control input, f : R n+1 → R n , g : R n+1 → R n are smooth functions. Let the C 2 function s : R n → R be an auxiliary output of the control system defining the sliding manifold. Following the classical methodology of the second order sliding mode design [START_REF] Levant | Construction principles of 2-sliding mode design[END_REF]), the control have to constrain the system trajectories to evolve onto the augmented sliding surface:

S := {z ∈ R n : s(z) = ṡ(z) = 0} .
In addition, the reaching phase of the surface S is usually required to be realized in a finite-time [START_REF] Levant | Construction principles of 2-sliding mode design[END_REF], [START_REF] Polyakov | Stability notions and lyapunov functions for sliding mode control systems[END_REF]).

Let us assume that the relative degree of s with respect to control input is equal to 2, then the output dynamics has the form s(t) = d(t, z) + c(t, z)u(t), where d : R n+1 → R and c : R n+1 → R. The functions d and c are traditionally assumed to be unknown but bounded by some numbers as follows:

|d(t, z)| ≤ d max , 0 < c 1 ≤ c(t, z) ≤ c 2 , ∀t ∈ R + , z ∈ R n . (2)
Denote x 1 = s, x 2 = ṡ and consider the following extended differential inclusion

ẋ(t) ∈ Ax(t) + b ([c 1 , c 2 ]u(t) + [-d max , d max ]) , (3)
where x = x 1 x 2 ∈ R 2 , A = 0 1 0 0 , b = 0 1 .
Evidently, the control law, which will stabilize the origin of the system (3), will provide the required second order sliding mode on the surface S to the system (1).

Finally, let us assume that only the output s and its derivative ṡ can be measured and used for the control purposes.

The problem is to design a feedback control u = u(s, ṡ) that stabilize the origin of the system (3) in a finitetime. In addition, the origin of the system (3) must be the unique discontinuity point of the control law to be designed (similarly to the so-called quasi-continuous sliding mode algorithm Levant [2005b]).

PRELIMINARIES

Finite-time stability

Consider the system of the form ẋ

(t) ∈ F (x(t)), t ∈ R + , x(0) = x 0 , (4)
where x ∈ R n is the state vector, F : R n → 2 R n is a nonlinear set-valued vector field.

Assume that F satisfies the standard conditions, i.e. F (x) is nonempty, compact and convex for any x ∈ R n and F is upper semi-continuous ( see [START_REF] Filippov | Differential equations with discontinuous right-hand sides[END_REF]).

An absolutely continuous function x(•, x 0 ) defined on [0, t * ), t * ≤ +∞, is called the solution to the Cauchy problem associated to (4) if x(0, x 0 ) = x 0 and x(•, x 0 ) satisfies the differential inclusion (4) for almost everywhere on (0, t * ).

Assume that the origin is an equilibrium point of the system (4), i.e. 0 ∈ F (0).

The paper studies only strong global stability properties of the origin of the system (4). We will omit the words strong and global in all later definitions and considerations. Definition 1. [START_REF] Orlov | Finite time stability and robust control synthesis of uncertain switched systems[END_REF], [START_REF] Polyakov | Stability notions and lyapunov functions for sliding mode control systems[END_REF]). The origin of system ( 4) is said to be uniformly finitetime stable if it is asymptotically stable and there exists a locally bounded function T : R n → R + ∪ {0} such that x(t, x 0 ) = 0 for all t > T (x 0 ), where x(•, x 0 ) is any solution of the Cauchy problem (4).

The scalar first order sliding mode control system is the classical example of the finite-time stable differential inclusion: ẋ(t) ∈ -2sign[x(t)]+[-1, 1], where x ∈ R, t > 0, x(0) = x 0 and sign[•] is the set-valued extension of the sign function based on the standard Fillipov regularization procedure (see, for example, [START_REF] Filippov | Differential equations with discontinuous right-hand sides[END_REF], [START_REF] Polyakov | Stability notions and lyapunov functions for sliding mode control systems[END_REF]). Obviously, the settling time function of this system is T (x 0 ) = |x 0 |. Theorem 2. [START_REF] Roxin | On finite stability in control systems[END_REF], [START_REF] Polyakov | Stability notions and lyapunov functions for sliding mode control systems[END_REF]). Let a function V : R n → R be continuous, radially unbounded, continuously differentiable outside the origin and sup

y∈F (x) ∇ T V (x)y ≤ -1 for x ∈ R n \{0}. (5)
Then the origin of the system (4) is finite-time stable and

T (x 0 ) ≤ V (x 0 ) for x 0 ∈ R n , ( 6
)
where T is a settling-time function.

The function V satisfying the conditions of Theorem 2 is known as the finite-time Lyapunov function.

Note that the settling-time function may coincide with a finite-time Lyapunov function for the system (4). For instance, the settling-time function T (x) = |x| from the example considered above is evidently the finitetime Lyapunov function for the corresponding first order sliding mode control system. The Lyapunov function to be designed in this paper will have the same property.

Frequently, finite-time stability correlates with the socalled homogeneity property of the system (see, for example, [START_REF] Rosier | Homogenous lyapunov function for homogenous coninuous vector field[END_REF], [START_REF] Bhat | Finite-time stability of continuous autonomous systems[END_REF], [START_REF] Orlov | Finite time stability and robust control synthesis of uncertain switched systems[END_REF], Levant [2005a], [START_REF] Andrieu | Homogeneous approximation, recursive obsertver and output feedback[END_REF]).

Homogeneity

Homogeneity is an intrinsic property of an object which remains consistent with respect to some scaling: level sets (resp. solutions) are preserved for homogeneous functions (resp. vector fields).

Let λ > 0, r i > 0, i ∈ {1, . . . , n} then one can define: [START_REF] Zubov | On systems of ordinary differential equations with generalized homogenous right-hand sides[END_REF], Levant [2005a]). A vector field F : R n → 2 R n is said to be r-homogeneous of degree m iff for all λ > 0 and for all x ∈ R n we have

• the vector of weights r = (r 1 , . . . , r n ) T , • the dilation matrix D r (λ) = diag{λ ri } n i=1 , (7) note that D r (λ)x = (λ r1 x 1 , . . . , λ ri x i , . . . , λ rn x n ) T . Definition 3.
λ -m D -1 r (λ)F (D r (λ)x) = F (x).
For a given x ∈ R n , the set (D r (λ)x) λ>0 is a curve on R n . An object is homogeneous iff the behavior of this object is symmetric along these particular curves. Theorem 4. (Levant [2005a], [START_REF] Bernuau | On extension of homogeneity notion for differential inclusions[END_REF]). Let a set-valued vector field F : R n → R n is r-homogenous of negative degree and the origin of the system (4) is asymptotically stable, then it is finite-time stable.

Homogeneity automatically provides additional robustness (Input-to-State Stability ) properties to the closed-loop control system (see [START_REF] Bernuau | On homogeneity and its application in sliding mode control[END_REF]). That is why, this paper designs the homogenous sliding mode control.

Homogeneous Implicit Lyapunov Function

The theorem given below refines Theorem 2 for implicit definition of Lyapunov function. Theorem 5. [START_REF] Polyakov | Finite-time stabilization using implicit lyapunov function technique[END_REF]). Let there exists a continuous function

Q : R + × R n → R (V, x) → Q(V, x) such that C1 ) Q is continuously differentiable for all x ∈ R n \{0} and V ∈ R + ; C2 ) for any x ∈ R n \{0} there exist V ∈ R + such that Q(V, x) = 0; C3 ) for Ω = (V, x) ∈ R n+1 : Q(V, x) = 0 we have lim x→0 (V,x)∈Ω V = 0 + , lim V →0 + (V,x)∈Ω x = 0, lim x →∞ (V,x)∈Ω V = +∞; C4 ) the inequality -∞ < ∂Q(V,x) ∂V < 0 holds for all V ∈ R + and x ∈ R n \{0}; C5 ) the inequality sup y∈F (x) ∂Q(V,x) ∂x ≤ ∂Q(V,x)
∂V holds for all (V, x) ∈ Ω. Then the origin of system (4) is globally finite-time stable with the following settling time estimate

T (x 0 ) ≤ V 0 , where V 0 ∈ R + such that Q(V 0 , x 0 ) = 0.
The proof of this theorem is based on combination of Theorem 2 and the classical Implicit Function Theorem [START_REF] Courant | Introduction to calculus and analysis[END_REF] providing

∂V ∂x = -∂Q ∂V -1 ∂Q ∂x .
Let us consider the following ILF candidate [START_REF] Nakamura | Stabilization of homogeneous systems using implicit control lyhapunov functions[END_REF], [START_REF] Polyakov | Finite-time stabilization using implicit lyapunov function technique[END_REF]) Q(V, x) := x T D r (V -1 )P D r (V -1 )x -1, (8) where V ∈ R + , x ∈ R n , D r (•) is the dilation matrix of the form (7) with r = (r 1 , r 2 , ..., r n ) T ∈ R n + , P ∈ R (n+1)×(n+1) is a symmetric positive definite matrix.

The proof of the next corollary can be found in Appendix. Corollary 6. Let the set-valued vector field F : R n → R n be r-homogeneous of degree -1 with the vector of weights r = (r 1 , r 2 , ..., r n ) T ∈ R n + . If there exists a matrix P ∈ R n×n , which satisfy matrix inequalities diag(r)P + P diag(r) > 0, P > 0 (9) and the inequality sup

y∈F (z)
z T P y + y T P z ≤ -z T (diag(r)P + P diag(r))z, ( 10) holds for all z ∈ R n such that z T P z = 1, then the system (4) is finite-time stable and the settling time function T : R n → R + ∪ {0} admits the following estimate:

T (x 0 ) ≤ V 0 (11) where x 0 ∈ R n and V 0 ∈ R + is such that Q(V 0 , x 0 ) = 0. If the inequality (10) becomes equality then T (x 0 ) = V (x 0 ).
This corollary has the simple geometric interpretation: if the homogeneous system (4) admits a Lyapunov function with at least one ellipsoidal level set, then it has a homogeneous ILF defined by the formula (8).

The present paper develops such quasi-continuous homogenous second order sliding mode control algorithm that provides to the closed-loop system a homogeneous ILF with ellipsoidal level sets ("quadratic ILF"). Corollary 6 will be used for the ILF design.

DESIGN OF NEW QUASI CONTINUOUS SECOND ORDER SLIDING MODE CONTROL

Originally, the second order sliding mode control was designed based on some geometrical constructions (Levant [2005a], [START_REF] Levant | Construction principles of 2-sliding mode design[END_REF]). Next, the traditional stability analysis based on Lyapunov function approach has been developed [START_REF] Polyakov | Lyapunov function design for finite-time convergence analysis: "twisting" controller for second-order sliding mode realization[END_REF], [START_REF] Moreno | Strict lyapunov functions for the super-twisting algorithm[END_REF]. This paper introduces "the reverse procedure", which designs the Lyapunov function before the control law. Note that the scheme presented below completely differs from the Control Lyapunov Function approach.

Design of the Lyapunov function

Let us define the ILF candidate for the system (3) in the form ( 8), with dilation matrix D r (λ) = λ 2 0 0 λ for λ > 0. Such selection of the dilation matrix is motivated by the structure of the system (3). The right-hand side of the open-loop system (3) satisfies Definition 3 with homogeneity degree -1.

For n = 2 the equation Q(V, x) = 0 can be rewritten as follows

V 4 -p 22 x 2 2 V 2 -2p 12 x 1 x 2 V -p 11 x 2 1 = 0, ( 12
) where x = (x 1 , x 2 ) T ∈ R 2 and {p ij } are elements of the matrix P > 0. In order to find the finite-time Lyapunov function, the matrix P ∈ R 2×2 should satisfy the conditions of Corollary 6. In this case, there exists a unique smooth (outside the origin) positive definite function V satisfying the equation ( 12) for all x ∈ R 2 (see the proof of Corollary 6). For n = 2 the function V can be found analytically. We use the procedure from [START_REF] Carpenter | On the solution of the real quartic[END_REF] for this purpose (see, Sections 8.2-8.4 for the details) Proposition 7. If the matrix P ∈ R 2×2 satisfies the linear matrix inequalities (9) with r = (2, 1) T , then 1) the smooth (outside the origin) positive definite function V : Q(V, x) = 0, x ∈ R n \{0} admits the following explicit representation:

V (x) =                          2p22x 2 2 + 4p 12 x 1 x 2 √ z 0 (x) -z0(x)+ √ z0(x) 2 if p 12 x 1 x 2 > 0, 2p22x 2 2 - 4p 12 x 1 x 2 √ z 0 (x) -z0(x)- √ z0(x) 2 if p 12 x 1 x 2 < 0, p22x 2 2 + √ p 2 22 x 4 2 +4p11x 2 1 2 if p 12 = 0, |x 1 | √ p 11 if x 2 = 0, |x 2 | √ p 22 if x 1 = 0, ( 13
) where z 0 (x) = 2p22x 2 2 -C1(x)-C2(x) 3 , C i (x) = 3 (-1) i √ ∆ 2 1 (x)-4∆ 3 0 (x)+∆1(x) 2 , i = 1, 2, ∆ 0 (x) = p 2 22 x 4 2 -12p 11 x 2 1 , ∆ 1 (x) = 2p 22 x 2 2 (∆ 0 + 6Dx 2 1
) and D = 8p 11 -9p 2 12 /p 22 > 0; 2) the function V defined by ( 13) with P = X -1 is homogeneous, positive definite, radially unbounded, continuous in R 2 and continuously differentiable outside the origin.

The proof is given in the Appendix.

Feedback Design

Based on the ideas of the paper [START_REF] Polyakov | Finite-time stabilization using implicit lyapunov function technique[END_REF] the quasi-continuous homogeneous second order sliding mode control for the system (3) is designed as follows u

(x) = kD r V -1 (x) x, (14)
where x ∈ R 2 , D r (•) is the dilation matrix for r = (2, 1) T , V is defined by ( 13) for P ∈ R 2×2 satisfying (9) and k = (k 1 , k 2 ) ∈ R 1×2 is a constant row vector. Theorem 8. Let the linear matrix inequalities (LMI) ÃX + X ÃT + c i (by

+ y T b T ) + d max (X + bb T ) ≤ 0, Xdiag(r) + diag(r)X > 0, X > 0, i = 1, 2 (15
) be feasible for some X = X T ∈ R n×n , y ∈ R 1×n and à = 2 1 0 1 then 1) the control function ( 14) with k = yX -1 is continuous outside the origin and globally bounded as follows

u 2 (x) ≤ k T k λmin(P )
for all x ∈ R 2 ;

2) the closed-loop system (3), ( 14) is homogenous and finite-time stable with the settling-time function estimate:

T (x) ≤ V (x).
Proof.

1) The equation Q(V, z) = 0 implies

x T D r (V -1 )P D r (V -1 )x = 1 ⇒ D r (V -1 )x 2 ≤ 1 λmin(P ) . and u(x) ≤ k • D r (V -1 )x ≤ k √ λmin(P )
.

2) Let us denote the right-hand side of the closed-loop system (3), ( 14) by

F (x) = Ax + b ([c 1 , c 2 ]u(x) + [-d max , d max ]
) . Taking into account homogeneity of the Lyapunov function function we easily derive u(D r (λ)x) = u(x) and

F (D r (λ)x) = λ -1 D r (λ)F (x), ∀λ ∈ R + .
Let z ∈ R n be such that z T P z = 1. The equality Q(V, z) = 0 implies V (z) = 1, u(z) = kz and sup y∈F (z)

y T P z + z T P y = z T (A T P + P A)z+ sup α∈[c1,c2] α(z T P bkz + z T k T b T P z)+ sup β∈[-dmax,dmax] β(z T P b + b T P z) = sup λz T (A T P + P A + c 1 (P bk + k T b T P ))z+ (1 -λ)z T (A T P + P A + c 2 (P bk + k T b T P ))z+ sup β∈[-dmax,dmax] β(z T P b + b T P z)
The LMI (15) implies for i = 1, 2

A T P + P A + c i (P bk + k T b T P ) ≤ -diag(r)P -P diag(r) -d max (P + P bb T P ) < 0.

Since -z T P z -zT P bb T P z ± (z T P b + b T P z) ≤ 0 for z T P z = 1 then sup y∈F (z)

y T P z + z T P y ≤ -z T (diag(r)P + P diag(r))z.

Applying Corollary (6) we finish the proof.

Note that if c = c 1 = c 2 and d max = 0 then replacing the first inequality in (15) with ÃX + X Ã + c(by + y T b T ) = 0 and repeating the proof of Theorem 8 it can be shown that T (x) = V (x), i.e. the Lyapunov function (13) of the closed-loop system (1), ( 14) coincides with the settlingtime function in the disturbance-free case. Remark 9. The condition u(x) ≤ u 0 , u 0 ∈ R + can be represented in the form of LMI

X y T y u 2 0 ≥ 0. ( 16
)
Indeed, the inequality u 2 =x T D r (V -1 )k T kD r (V -1 )x ≤ u 2 0 =u 2 0 x T D r (V -1 )P D r (V -1 )x holds for all V ∈ R + and x ∈ R n such that Q(V, x) = 0, then k T k ≤ u 2 0 P .

Remark on chattering reduction

Practical importance of the high order sliding mode control concept is usually motivated by essential reduction of the chattering phenomenon. In practice, any high gain control (even linear) may be subjected by chattering. The Lyapunov function method provides the additional tool for chatting attenuation. Following classical ideas of the chattering reduction (see [START_REF] Utkin | Sliding Mode Control in Electro-Mechanical Systems[END_REF]) we linearize our control close to sliding surface, i.e.

u(x) = kD(V -1 min )x for V (x) ≤ V min kD V -1 (x) x for V (x) ≥ V min (17)
where V min > 0 is a parameter. Note that the level set of the Lyapunov function defined by inequality V (x) ≤ V min is an ellipsoid in R n , which is positively invariant for the closed-loop system. So, in order to reduce the chattering we replace the quasi-continuous second order sliding mode control law with linear one inside the ellipsoid {x ∈ R n : V (x) ≤ V min }. The size of the ellipsoid is defined by V min .

NUMERICAL SIMULATIONS

For c 1 = 1, c 2 = 1.2, d max = 0.2 and u 0 = 2 we solve the LMI (15): P = 23.7436 10.0502 10.0502 4.7858 and k = (-5.9663, -3.6786).

We design the sliding mode control in the form (17) with V min = 0.05. Figure 1 presents the simulation results for the system

ẋ(t) = Ax(t) + b ((1 + 0.2 sin x 1 )u(t) + 0.2 sin x 2 ) .
The simulation has be done using explicit Euler method with the step size 0.01.

CONCLUSION

The homogeneous quasi-continuous second order sliding mode algorithm is developed. The LMI-based algorithm for parameter tuning is elaborated. The robust finite-time stability analysis based on Lyapunov function method is provided. The chattering reduction scheme based on Lyapunov function method is presented and demonstrated on numerical simulations. I. The function Q(V, x) defined by ( 8) satisfies the conditions C1 )-C3 ) of Theorem 5. Indeed, it is continuously differentiable for all V ∈ R + and ∀x ∈ R n .

Since P > 0 then the following chain of inequalities λmin(P ) x 2 max i {V 2r i } ≤ Q(V, x) + 1 ≤ λmax(P ) x 2 min i {V 2r i } implies that for any z ∈ R n \{0} there exist V -∈ R + and V + ∈ R + : Q(V -, x) < 0 < Q(V + , x). Due to continuity of Q the condition C2) also fulfilled.

If Q(V, x) = 0 then the same chain of inequalities gives

min i {V 2r i } λmax(P ) ≤ x 2 ≤ max i {V 2r i }
λmin(P ) . Therefore, the condition C3 ) of Theorem 5 holds.

Since ∂Q ∂V = -x T Dr(V -1 )(diag(r)P +P diag(r))Dr(V -1 )x V , then diag(r)P + P diag(r) > 0 implies ∂Q ∂V < 0 for ∀V ∈ R + and x ∈ R n \{0}. So, the condition C4 ) of Theorem 5 also holds.

We conclude that the equation Q(V, x) = 0 implicitly defines a proper positive definite function V : R n → R + ∪ {0} such that for any x ∈ R n we have Q(V (x), x) = 0. (z T P D r (V -1 )y + y T (t, x 0 )D r (V -1 )P z), where x ∈ R n \{0} and z = D r (V -1 )x. The set-valued vector field F is homogeneous, i.e. λ -m D -1 r (λ)F (D r (λ)x) = F (x), ∀λ ∈ R + and ∀x ∈ R n , where m ∈ R is the degree of homogeneity. Hence, V -1 D -1 r (V -1 )F (z) = F (x) for λ = V -1 and V (x) ≤ sup y∈F (z)

-∂Q ∂V

Fig

  Fig. 1. The simulation results 8. APPENDIX 8.1 Proof of Corollary 6

  II. Let us calculate the time derivative of the function V along the trajectories of the system (4

V -1 (z T P y + y T P z) = sup y∈F (z) z T P y+y T P z z T (diag(r)P +P diag(r))z . Since Q(V, x) = 0 then z T P z = 1 and the inequality (10) implies V (x) ≤ -1 for x = 0.

Cubic equation

Consider the cubic equation z 3 + pz 2 + qz + r = 0, p, q, r ∈ R. Introduce the following numbers

where

is the real root of the cubic equation.

Quartic equation

Let us consider the quartic equation V 4 + aV 2 + bV + c = 0, a, b, c ∈ R and the adjoint cubic equation

is a root of the cubic equation the roots of the quartic one are the following

, where D = 8p 11 -9p 2 12 /p 22 > 0. On the on hand, if ∆ 0 ≤ 0 then ∆ 2 1 -4∆ 3 0 ≥ 0. On the other hand, since

Therefor the real non-negative root of the cubic equation is the following:

where C 1 and C 2 are represented by (18).

III. Let us study the case z 0 (x 1 , x 2 ) = 0 or, equivalently, 2p 22 x 2 2 = C 1 + C 2 . The last equality can be rewritten as follows:

∆ 1 + 6∆ 0 p 22 x 2 2 = 8p 3 22 x 6 2 . Hence, z 0 (x 1 , x 2 ) = 0 is equivalent to x 11 x 22 = 0.

IV. Finally taking into account the formulas for the roots of the quartic equation we finish the proof of the formula (13).

V. Since the inequality (9) holds for P = X -1 , then the Lyapunov function candidate V is well-defined by the equation Q(V, x) = 0 (see, the proof of Corollary 6). Since Q(V, D r (λ)x) = Q(λV, x) and 0 = Q(V (D r (λ)x), D r (λ)x) then λV (x) = V (D r (λ)x), ∀λ ∈ R + , i.e. V is homogeneous function.


Philip A Knight
email: p.a.knight@strath.ac.uk

Daniel Ruiz
email: daniel.ruiz@enseeiht.fr

Bora Uçar
email: bora.ucar@ens-lyon.fr

AND BORA U C ¸ar

A Symmetry Preserving Algorithm for Matrix Scaling

Keywords: sparse matrices, matrix scaling, equilibration AMS subject classifications. 05C50, 65F35, 65F50

We present an iterative algorithm which asymptotically scales the ∞-norm of each row and each column of a matrix to one. This scaling algorithm preserves symmetry of the original matrix and shows fast linear convergence with an asymptotic rate of 1/2. We discuss extensions of the algorithm to the 1-norm, and by inference to other norms. For the 1-norm case, we show again that convergence is linear, with the rate dependent on the spectrum of the scaled matrix. We demonstrate experimentally that the scaling algorithm improves the conditioning of the matrix and that it helps direct solvers by reducing the need for pivoting. In particular, for symmetric matrices the theoretical and experimental results highlight the potential of the proposed algorithm over existing alternatives.

Introduction.

Scaling a matrix consists of pre-and postmultiplying the original matrix by two diagonal matrices. We consider the following scaling problem: given a large, sparse matrix A ∈ R m×n , find two positive diagonal matrices D and E such that all rows and columns of the scaled matrix A = DAE have the same length in the ∞-norm. We propose an iterative algorithm for this purpose which preserves symmetry (when present), and we investigate its convergence properties. We discuss the extension of the algorithm to the 1-norm, and describe how it can be used for other norms as well.

Scaling or equilibration of data in linear systems of the form Ax = b can provide significant computational benefits. In this setting, we form a scaled matrix A = DAE and solve the equation Ax = b, where x = E -1 x and b = Db, with the expectation that the scaled system has better conditioning and hence provides more accurate results. Scaling can also be helpful in designing algorithms for linear systems. For example, Bunch and Parlett [START_REF] Bunch | Direct methods for solving symmetric indefinite systems of linear equations[END_REF] prove the stability of their factorization method by assuming that the rows (hence the columns) of the symmetric input matrix have ∞-norm equal to 1.

There are several well-known algorithms (see, for example, [16, sect. 4.12], [20, sect. 3.5.2], and [START_REF] Schneider | A comparative study of algorithms for matrix balancing[END_REF]) for scaling matrices. The row and column scaling methods are amongst the simplest. The row scaling method divides each row in the original matrix by its norm. Different norms, such as the ∞-norm or the 1-norm, may be considered, depending on the application. Column scaling is equivalent to row scaling on A T . A more elaborate scaling method is proposed by Curtis and Reid [START_REF] Curtis | On the automatic scaling of matrices for Gaussian elimination[END_REF] and is implemented as MC29 in the HSL Mathematical Software Library. This method aims to make the nonzeros of the scaled matrix close to one by minimizing the sum of the squares of the logarithms of the moduli of the nonzeros. MC29 reduces that sum in a global sense and can provide benefits over a wide range of sparse matrices. The routine MC30 in the HSL library is a variant of MC29 for symmetric matrices. Scaling can also be combined with permutations (see [START_REF] Duff | The design and use of algorithms for permuting large entries to the diagonal of sparse matrices[END_REF] and the HSL routine MC64). A matrix is first permuted to maximize the product of the absolute values of the entries on the diagonal. Then the matrix is scaled so that the diagonal entries equal one and the off-diagonals are less than or equal to one. This has been shown [START_REF] Duff | On algorithms for permuting large entries to the diagonal of a sparse matrix[END_REF] to be useful for finding a good sequence of pivots for sparse direct solvers and for building good incomplete LU preconditioners for iterative methods.

In the 1960s, optimality properties were established in terms of condition numbers for scaled matrices all of whose rows or columns have norm one [START_REF] Bauer | Optimally scaled matrices[END_REF][START_REF] Bauer | Remarks on optimally scaled matrices[END_REF][START_REF] Van Der Sluis | Condition numbers and equilibration of matrices[END_REF]. In particular, the optimal scaling of a matrix A which minimizes the condition number in the ∞norm is characterized by both DAE and E -1 A -1 D -1 having equal row sums of absolute values [START_REF] Bauer | Optimally scaled matrices[END_REF]. Other optimality results for one sided scaling, i.e., DA or AE are also shown [START_REF] Van Der Sluis | Condition numbers and equilibration of matrices[END_REF], again based on the equivalence of the norms of rows or columns in the ∞and 1-norms.

This paper is organized as follows. In section 2 we introduce the proposed algorithm, derive a convergence analysis, and show some particular properties it yields for scaled matrices. In section 3, we consider the algorithm in norms other than the ∞norm. Following a discussion in [START_REF] Parlett | Methods for scaling to double stochastic form[END_REF], we establish precise conditions for convergence in the 1-norm, and we look at the rate of convergence, too. We present numerical results in section 4 and finish with some concluding remarks in section 5. While proving the convergence of the proposed algorithm in 1-norm, we deviate somewhat from the approach in [START_REF] Parlett | Methods for scaling to double stochastic form[END_REF]. The appendix includes a detailed explanation of the differences. , where • ∞ stands for the ∞-norm of a real vector (that is the maximum entry in absolute value; sometimes called the max norm). If a row (or a column) in A has all entries equal to zero, we replace the diagonal entry in R (or C, respectively) by 1. In the following, we will assume that this does not happen. One can scale the matrix A on both sides, forming the scaled matrix A as (2.2)

A = R -1 AC -1 .
The idea of the algorithm we propose is to iterate on that process, resulting in Algorithm 1. The algorithm converges when (2.3)

max 1≤i≤m |1 -r (k) i ∞ | ≤ ε and max 1≤j≤n |1 -c (k) j ∞ | ≤ ε
for a given value of error tolerance parameter ε > 0.

Algorithm 1 Simultaneous row and column scaling in the ∞-norm.

1: A (0) ← A 2: D (0) ← Im 3: E (0) ← In 4: for k = 0, 1, 2, . . . until convergence do 5:

R ← diag r (k) i ∞ i=1,...,m r (k) i
is the ith row of A (k) 6:

C ← diag c (k) j ∞ j=1,...,n c (k) j
is the jth column of A (k) 7:

A (k+1) ← R -1 A (k) C -1
8:

D (k+1) ← D (k) R -1
9:

E (k+1) ← E (k) C -1
2.1. A salient property. We highlight that the proposed iterative scaling procedure preserves the symmetry in the original matrix. In fact, this is one of our main motivations. If the given matrix A is symmetric, then the diagonal matrices R and C in (2.1) are equal and, consequently, matrix A in (2.2) is symmetric, as is the case for the matrices A (k) at any iteration in Algorithm 1. This is not the case for most scaling algorithms which alternately scale rows followed by columns or vice-versa.

In the case of unsymmetric matrices, one may consider the use of the Sinkhorn-Knopp iterations [START_REF] Sinkhorn | Concerning nonnegative matrices and doubly stochastic matrices[END_REF] with the ∞-norm in place of the 1-norm. This method simply normalizes all rows and then all columns of A, and iterates on this process until convergence. In ∞-norm, the convergence is achieved after one single step. Because of its simplicity, this method is very appealing. There are, however, differences. The Sinkhorn-Knopp iteration may provide different results when applied to A T . On the other hand, Algorithm 1 provides exactly the same results when applied to A or A T in the sense that the scaled matrix obtained on A T is the transpose of that obtained on A. We have quoted the Sinkhorn-Knopp method, in particular, because it has been originally proposed to obtain doubly stochastic matrices (that is nonnegative matrices with all rows and columns of 1-norm equal to one), and we shall come back to this issue with respect to Algorithm 1 in section 3.

Convergence rate.

The proposed algorithm converges quickly to a matrix whose rows and columns have ∞-norm equal to one as stated below.

Theorem 2.1. If A is an m × n matrix with no zero rows or columns, then Algorithm 1 produces a sequence of matrices A (0) , A (1) , . . . which converge linearly to a matrix all of whose rows and columns have infinity norm equal to one. The asymptotic rate of convergence is 1/2.

Proof. After the first iteration of the algorithm, all the entries in A (1) are less than or equal to one in absolute value. This is because all entries a ij in A are divided by the square roots of two numbers, r

(k) ip | = r (k) i ∞ and |a (k) j | = c (k) j ∞ .
With this notation, we can easily verify that both entries a (k+1) ip and a (k+1) j in the scaled matrix A (k+1) are greater, in absolute value, than the square root of the corresponding value at iteration k, and are still less than one. We can write

1 ≥ |a (k+1) ip | = |a (k) ip | r (k) i ∞ c (k) p ∞ = |a (k) ip | c (k) p ∞ ≥ |a (k) ip | since |a (k) ip | = r (k) i ∞ and c (k) p ∞ ≤ 1 for any k ≥ 1. Similarly, |a (k) j | ≤ |a (k+1) j | ≤ 1 for any k ≥ 1.
From this, we can write that the iterations in Algorithm 1 provide scaled matrices A (k) , k = 1, 2, . . . with the following properties:

(2.4) for all k ≥ 1, 1 ≤ i ≤ m, r (k) i ∞ ≤ |a (k+1) ip | ≤ r (k+1) i ∞ ≤ 1, and
(2.5) for all k ≥ 1, 1 ≤ j ≤ n, c (k) j ∞ ≤ |a (k+1) j | ≤ c (k+1) j ∞ ≤ 1,
which shows that both row and column norms must converge to 1. To conclude our demonstration, we just need to see that for all i

1 -r (k+1) i ∞ = 1 -r (k+1) i 2 ∞ 1 + r (k+1) i ∞ ≤ 1 -r (k) i ∞ 1 + r (k+1) i ∞ ,
and that similar equations hold for the columns as well.

A small example [START_REF] Ruiz | A Scaling Algorithm to Equilibrate Both Rows and Columns Norms in Matrices[END_REF] shows that this asymptotic rate is sharp. Let us consider the following 2 × 2 matrix with a badly scaled row:

A = α α 1 1 . If α 1, then iteration k (k ≥ 1)
of the algorithm provides the following matrices:

D (k) = α -(1-1 2 k ) 0 0 1 , A (k) = α 1 2 k α 1 2 k 1 1 , E (k) = I 2 .
When the algorithm converges, the scaled matrix A is the matrix with all ones, D has its first diagonal entry equal to α -1 , and E stays as the identity matrix. In this case, the linear rate of 1 2 is met.

Comparison with Bunch's algorithm.

As we have stated before, Algorithm 1 is well suited for symmetric scaling of symmetric matrices. For the ∞-norm case, Bunch [START_REF] Bunch | Equilibration of symmetric matrices in the max-norm[END_REF] also developed an efficient symmetric scaling algorithm. We highlight some differences between the two algorithms. Bunch's algorithm processes the rows of the lower triangular part of the given matrix in such a way that the largest entry seen in each row is made to be ±1 in the partially scaled matrix. At the end, the scaling found for the ith row is applied to the ith column of the whole matrix. Bunch's algorithm runs in O(nnz)-time, equivalent to one iteration of the proposed algorithm. Since the scaling value of a row is calculated using the scaled values of its nonzero entries in the lower triangular part, there is a dependency between the computation of the scaling values. The dependencies can be an obstacle to reduce the running time in a parallel computing environment. For example, Bunch's algorithm has to run sequentially to scale a tridiagonal matrix. On the other hand, Algorithm 1 is more amenable to parallelism (as was shown before [START_REF] Amestoy | A parallel matrix scaling algorithm[END_REF]). One may consider permuting a matrix to take advantage of parallelism. But care should be taken, as Bunch's algorithm is sensitive to symmetric permutations applied to the original matrix (any diagonal nonzero can be symmetrically permuted to the first position so that that entry is one in the scaled matrix). On the other hand, the proposed algorithm is independent of any permutations (not even necessarily symmetric) applied to the original matrix in the sense that the scaling matrices computed for the permuted matrix would be equal to the permuted scaling matrices computed for the original matrix.

Extensions to other norms.

A natural idea would be to change the norm used in Algorithm 1, and to try the 2-norm or the 1-norm because of the optimal properties they induce (see [START_REF] Parlett | Methods for scaling to double stochastic form[END_REF][START_REF] Van Der Sluis | Condition numbers and equilibration of matrices[END_REF]), and still expect convergence (where all rows and columns are of length 1 in the corresponding norm). We shall see, in the remainder of this section, that this will usually, but not always, work and we investigate the potential and limitations of such extensions. Compared to the ∞-norm, where different algorithms can raise very different solutions, scalings in the 1-norm or the 2-norm, when they exist, can be considered as unique as the scaled matrix is essentially unique.

With respect to the extension of Algorithm 1 to the scaling of matrices in other norms, the case of the 1-norm is central. Indeed, Rothblum, Schneider, and Schneider [26, page 13] showed that the problem of scaling a matrix A in the p -norm, for 1 < p < ∞, can be reduced to the problem of scaling the pth Hadamard power of A, i.e., the matrix A [p] = [a p ij ], in the 1-norm. We can apply that discussion to Algorithm 1 by replacing the matrix A with A [p] , and then by taking the Hadamard pth root, e.g., D

[1/p] 1 = [d 1/p
ii ], of the resulting iterates. For this reason, we shall analyze the convergence properties of Algorithm 1 for the 1-norm only, knowing that these will implicitly drive the conclusions for any of the p norms, for 1 < p < ∞.

Background.

The idea of equilibrating a matrix such that the 1-norm of the rows and columns are all 1 is not new, dating back to at least the 1930s (see some historical remarks by Knight [START_REF] Knight | The Sinkhorn-Knopp algorithm: Convergence and applications[END_REF]). Here, we briefly review some of the previous work.

Sinkhorn and Knopp [START_REF] Sinkhorn | Concerning nonnegative matrices and doubly stochastic matrices[END_REF] studied a method for scaling square nonnegative matrices to doubly stochastic form, that is a nonnegative matrix with all rows and columns of equal 1-norm. Sinkhorn [START_REF] Sinkhorn | A relationship between arbitrary positive matrices and doubly stochastic matrices[END_REF] originally showed the following: Any positive square matrix of order n is diagonally equivalent to a unique doubly stochastic matrix of order n, and the diagonal matrices which take part in the equivalence are unique up to scalar factors. Later, a different proof for the existence part of Sinkhorn's theorem with some elementary geometric interpretations was given [START_REF] Borobia | Matrix scaling: A geometric proof of Sinkhorn's theorem[END_REF].

Sinkhorn's result was further extended to the case of nonnegative, nonzero matrices [START_REF] Sinkhorn | Concerning nonnegative matrices and doubly stochastic matrices[END_REF]. A few definitions are necessary to state the result. A square n × n nonnegative matrix A ≥ 0 is said to have support if there exists a permutation σ such that a i,σ(i) > 0 for 1 ≤ i ≤ n. Note that matrices not having support are matrices for which no full transversal can be found (see [16, page 107]), i.e., there is no column permutation making the diagonal zero-free, and are thus structurally singular. A matrix A is said to have total support if every positive entry in A can be permuted into a positive diagonal with a column permutation. A nonnegative nonzero square matrix A of size n > 1 is said to be fully indecomposable if there does not exist permutation matrices P and Q such that PAQ is of the form ( Different contributions have also been made in the study of convergence of the Sinkhorn-Knopp method under various hypothesis. Geometric convergence rate for positive [START_REF] Sinkhorn | Diagonal equivalence to matrices with prescribed row and column sums[END_REF] and nonnegative matrices with total support [START_REF] Soules | The rate of convergence of Sinkhorn balancing[END_REF] have been shown. The converse of the second result has also been established [START_REF] Achilles | Implications of convergence rates in Sinkhorn balancing[END_REF], i.e., geometric convergence of the Sinkhorn-Knopp algorithm implies total support for a nonnegative matrix. The explicit rates of convergence for fully indecomposable matrices are given by Knight [START_REF] Knight | The Sinkhorn-Knopp algorithm: Convergence and applications[END_REF].

Parlett and Landis [START_REF] Parlett | Methods for scaling to double stochastic form[END_REF] present three iterative scaling algorithms with experimental evidence of better average and worst-case convergence behavior than the Sinkhorn-Knopp method (in at least one of the three algorithms). They also give a generalized version of the convergence theorem of Sinkhorn and Knopp [START_REF] Sinkhorn | Concerning nonnegative matrices and doubly stochastic matrices[END_REF], including a characterization of scaling algorithms that will converge to a doubly stochastic matrix when the starting matrix A has support. Such algorithms are called diagonal product increasing (DPI) algorithms. In the appendix, we recall and partially extend Parlett and Landis's results [START_REF] Parlett | Methods for scaling to double stochastic form[END_REF]Theorem 1]. These extensions are needed to specify generic properties that our algorithm fulfills which are sufficient to ensure convergence of scaling algorithms in the 1-norm in general.

In the remainder of this section, we first establish that Algorithm 1 converges in the 1-norm as well as in any p-norm, for any p < ∞, so long as A satisfies appropriate conditions. We also analyze the convergence rate of the algorithm in the 1-norm (separately for symmetric and unsymmetric matrices), allowing a fuller comparison with other balancing algorithms.

Convergence analysis.

Algorithm 1 produces a sequence of matrices diagonally equivalent to the starting matrix A = A (0) using the iteration (3.2) n, andc (k) j , j = 1, . . . , n, denote the 1-norm of rows and columns, respectively,

A (k) = a (k) ij = D (k) AE (k) , k = 1, 2, . . . , D (k) = diag d (k) 1 , . . . , d (k) n , E (k) = diag e (k) 1 , . . . , e (k) n , where D (0) = E (0) = I. For convenience, let r (k) i , i = 1, . . . ,
(3.3) r (k) i = r (k) i 1 = n j=1 |a (k) ij | , c (k) j = c (k) j 1 = n i=1 |a (k) ij | , so that d (k) i = 1/ r (k) i and e (k) j = 1/ c (k) j .
Since A and |A| scale in the same way, we can assume that A ≥ 0. Under this simplification, the 1-norms of the rows and columns reduce to the row and column sums, respectively, r

(k) i = n j=1 a (k) ij and c (k) j = n i=1 a (k)
ij , and to generalize our results to any matrix, one just needs to extend the definition of a doubly stochastic matrix so that the absolute value of the matrix under consideration is doubly stochastic.

Theorem 3.1. Given the sequence (3.2) of diagonal equivalents for A, in which

a (k+1) ij = a (k) ij r (k) i c (k) j , 1 ≤ i, j ≤ n, with r (k) i and c (k) j
given by (3.3), we have the following: 1. If A has support, then S = lim k→∞ A (k) exists and is doubly stochastic. 2. If A has total support, then both D = lim k→∞ D (k) and E = lim k→∞ E (k) exist and S = DAE. We will prove this theorem in a more general setting, following an approach by Parlett and Landis [START_REF] Parlett | Methods for scaling to double stochastic form[END_REF], by establishing that the results in Theorem 3.1 hold for any iterative scaling algorithm having the following properties:

c (k) j = 1, 1 ≤ i, j ≤ n.
The full consequences of these properties, not least that they imply convergence of DPI algorithms, are detailed in the appendix. We have made minor changes to the properties used by Parlett and Landis [START_REF] Parlett | Methods for scaling to double stochastic form[END_REF] so as to include Algorithm 1 in the analysis. In summary, our property (P1 * ) is more restrictive than that of Parlett and Landis who require the product of the products to be monotonically increasing. Property (P2) relaxes the two other hypothesis made by Parlett and Landis when defining DPI algorithms, but is sufficient, combined with (P1*), to obtain the conclusions in Theorem 3.1.

Before the proof, we note the following observations. By the arithmetic-geometric mean inequality, we have

n i=1 r (k+1) i ≤ 1 n n i=1 r (k+1) i n = ⎛ ⎜ ⎝ 1 n 1≤i,j≤n a (k) ij r (k) i c (k) j ⎞ ⎟ ⎠ n
for all k, with the same inequality for

n j=1 c (k+1) j
since the sum of the column sums must equal the sum of the row sums. Additionally, using the Cauchy-Schwarz inequality on the vectors a

(k) ij / r (k) i 1≤i,j≤n and a (k) ij / c (k) j 1≤i,j≤n , 1≤i,j≤n a (k) ij r (k) i c (k) j ≤ 1≤i,j≤n a (k) ij r (k) i 1≤i,j≤n a (k) ij c (k) j = √ n √ n ,
and thus, for all k ≥ 0, we have

(3.5) n i=1 r (k+1) i ≤ 1 n n i=1 r (k+1) i n ≤ 1 and n j=1 c (k+1) j ≤ ⎛ ⎝ 1 n n j=1 c (k+1) j ⎞ ⎠ n ≤ 1.
Notice that after the first iteration all the entries in matrix A (k) are less than or equal to 1, since both r

n i=1 d (k+1) i n i=1 d (k) i = 1 n i=1 r (k) i
(and similarly for the column scaling factors), we see that property (P1 * ) holds. Establishing (P2) is not quite so straightforward. From the arithmetic-geometric mean inequality and (3.5),

n i=1 r (k) i c (k) i ≤ n i=1 1 2n r (k) i + c (k) i 2n ≤ 1.
From (P1 * ) and the fact that all the elements a (k) ij for k ≥ 1 are bounded above by 1, we can conclude from Lemma A.2 of the appendix that lim k→∞ s k = ξ > 0 exists, where

s k = n i=1 d (k) i e (k)
i , and, consequently, that

lim k→∞ n i=1 r (k) i c (k) i = lim k→∞ s k s k+1 = 1.
We conclude that

(3.6) lim k→∞ n i=1 r (k) i c (k) i = 1 and lim k→∞ n i=1 1 2n r (k) i + c (k) i = 1.
Now, since all the elements in A (k) are less than 1 after the first iteration, each of the two sequences (r

(k) i ) k≥1 and (c (k) j ) k≥1 , for all 1 ≤ i, j ≤ n, are bounded. Let us introduce the sequence (v (k) ) k≥1 of the 2n-vectors v (k) = (r (k) 1 , . . . , r (k) n , c (k) 1 , . . . , c (k) n )
, which is also bounded in R 2n of finite dimension. Consider then any convergent subsequence (v (q) ) q, and let

x i = lim q→∞ r (q) i , 1 ≤ i ≤ n, and y j = lim q→∞ c (q) j , 1 ≤ j ≤ n. From (3.6), we can write n i=1 x i y i = n i=1 1 2n (x i + y i ) 2n = 1,
and since equality holds in the arithmetic-geometric mean inequality only if all the elements are equal, then

x 1 = • • • = x n = 1 = y 1 = • • • = y n .
Therefore, any convergent subsequence of the bounded sequence (v (k) ) k≥1 in finite dimensional space must have the same limit (a vector of ones), which implies that the sequence (v (k) ) k≥1 is necessarily convergent and that (3.4) holds.

Since our algorithm fulfills properties (P1 * ) and (P2), it converges. Additionally, since the computation of the scaling factors in Algorithm 1 is only based on the row and column sums in matrix A (k) at each iteration, we see that the algorithm is permutation insensitive, and it does not mix information from independent subsets of entries residing in different blocks of a decomposable matrix A. This can be used to show that the scaling matrices also exist (see Corollary A.4 and the accompanying discussion in the appendix) when A has total support.

Tracking the computations performed by Algorithm 1 it is straightforward to establish that it preserves symmetry.

Corollary 3.2.

1. If A is symmetric and has support, then Algorithm 1 in the 1-norm builds a sequence of symmetric scalings of A converging to a symmetric doubly stochastic limit. 2. If A is symmetric and has total support, then A is symmetrically equivalent to a symmetric doubly stochastic matrix S, and Algorithm 1 in the 1-norm builds a convergent sequence of diagonal matrices D (k) such that S = lim k→∞ D (k) AD (k) .

Rate of convergence for symmetric matrices.

Let e be a vector of ones (dimension should be clear) and D(x) = diag(x) for x ∈ R n . Note that x = D(x)e and D(D(x)y) = D(x)D(y) = D(y)D(x). When we take the square root of a vector this should be interpreted componentwise (see, e.g., (3.8)). We will also assume that A ≥ 0 and that it has no zero rows.

We first treat the case when A is symmetric. Then the 1-norm version of Algorithm 1 loops the following steps (where A (0) = A and E (0) = I): (k) and we can come up with a more compact form of the algorithm as follows. Let x k be the vector such that D(x k ) = E (k) . Then

(3.7) R = D(A (k) e) 1/2 , E (k+1) = R -1 E (k) , A (k+1) = R -1 A (k) R -1 . Note that A (k) = E (k) AE
E (k+1) = R -1 E (k) = D(E (k) AE (k) e) -1/2 E (k) = D(D(x k )Ax k ) -1/2 D(x k ) = D(x k ) -1/2 D(Ax k ) -1/2 D(x k ) = D(x k ) 1/2 D(Ax k ) -1/2 .
In other words, we can carry out the iteration simply working with the vectors x k . Namely, let x 0 = e (say) and form (3.8)

x k+1 = x k Ax k ,
where the division is performed componentwise. If A is unsymmetric, then we can replace it with

0 A A T 0
and extract the row and column scalings from the top and bottom half of x k .

The reason for introducing the compact form is to allow us to exploit some standard theory for fixed point iteration. In the 1-norm, the balancing problem can be viewed as an attempt to find the positive solution of the nonlinear equation

(3.9) D(x)Ax = e.
This can be rewritten in many ways. In particular, noting that D(x)e = D(Ax) -1 e, we can write

x = D(x)e = D(x) 2 e = D(x)D(Ax) -1 e = D(Ax) -1 x.
In other words (3.8) is a fixed point iteration for solving (3.9). Note that Sinkhorn-Knopp can also be framed in terms of a fixed point iteration for solving (3.9), although some care needs to be taken to ensure that the correct iterates are extracted [START_REF] Knight | The Sinkhorn-Knopp algorithm: Convergence and applications[END_REF]. The consequence of this discussion is that we can prove a result on the asymptotic convergence rate of the algorithm by bounding the norm of the Jacobian of f (x) = D(Ax) -1 x in the environs of a fixed point. We restrict x so that it lies in the positive cone R n + to ensure that all our diagonal matrices are invertible, and we only take the square roots of positive numbers. Not only is f (x) continuous, but it is twice differentiable. Differentiating f (x) term by term, we find that its Jacobian can be written as

(3.10) J(x) = - 1 2 D(x) 1/2 D(Ax) -3/2 A + 1 2 D(Ax) -1/2 D(x) -1/2 .
To establish the asymptotic rate of convergence, we want to bound J(x) at the fixed point x * . Substituting the identity D(Ax * ) = D(x * ) -1 into (3.10) gives

J(x * ) = - 1 2 D(x * ) 1/2 D(x * ) 3/2 A + 1 2 D(Ax * ) 1/2 D(x * ) -1/2 = 1 2 I -D(x * ) 2 A = 1 2 D(x * ) (I -D(x * )AD(x * )) D(x * ) -1 = 1 2 D(x * )(I -P)D(x * ) -1 ,
where P is symmetric and stochastic. In a neighborhood of x * we can bound the asymptotic rate of convergence of the iteration by the largest eigenvalue of (I -P)/2.

Since the spectrum of P lies in [-1, 1], the bound is (1λ min (P))/2. We can couple our observations with the convergence results from section 3 to make a strong statement about the performance of the algorithm.

Theorem 3.3. Let A ≥ 0 be symmetric and be fully indecomposable. Then for any vector x 0 > 0, the iteration (3.8) will converge to the unique positive vector x * such that D(x * )AD(x * ) = P is doubly stochastic and the convergence is asymptotically linear at the rate (1λ min (P))/2 < 1.

Proof. We already know that (3.8) converges from Theorem 3.1. The uniqueness of x * is well known (see, for example, Lemma 4.1 of Knight [START_REF] Knight | The Sinkhorn-Knopp algorithm: Convergence and applications[END_REF]).

Since A is fully indecomposable, so is P. Such a doubly stochastic matrix is primitive [START_REF] Berman | Nonnegative Matrices in the Mathematical Sciences[END_REF], and hence it has a single simple eigenvalue of modulus 1. Therefore, λ min (P) > -1 and the asymptotic rate of convergence, established above, is bounded below one.

Note that Theorem 3.1 establishes the convergence of (3.8) for matrices with total support. If a matrix has total support but is not fully indecomposable, then it can be permuted symmetrically into a block diagonal matrix where each block is irreducible [START_REF] Knight | A fast algorithm for matrix balancing[END_REF]Lemma 3.1]. Now if an irreducible block is also fully indecomposable, then we can apply Theorem 3.3 to it. Suppose that a block is irreducible but not fully indecomposable. It must be permutable to the form (3.1). Since no zero block in an irreducible matrix can be permuted symmetrically to lie in the bottom lefthand corner, the lack of indecomposability manifests itself in a diagonal block. That is, there is a symmetric permutation of the block to

0 B B T Z
where B is fully indecomposable. Since A has total support, Z = 0. For a matrix with this structure, we can still give the rate of convergence: it is covered by a general convergence result for (3.8) on unsymmetric matrices, which we establish next. instead: using (3.7) on this matrix is exactly the same as running Algorithm 1 on A. At a fixed point, the Jacobian matrix can be written as

Rate of convergence for unsymmetric matrices. Much of our analysis carries through if A is unsymmetric and we work with

(3.11) J(x * ) = 1 2 D(x * ) I - 0 P P T 0 D(x * ) -1 ,
where P is doubly stochastic but may be unsymmetric. There is a snag, though. Notice that

0 P P T 0 e -e = -e e ,
from which we conclude that -1 is an eigenvalue of J(x * ), and we can no longer apply the contraction mapping theorem. Nor is the fixed point unique: if P = DAE, then αD and E/α give another solution.

A similar problem arises in trying to apply the contraction mapping theorem in the analysis of the Sinkhorn-Knopp algorithm [START_REF] Knight | The Sinkhorn-Knopp algorithm: Convergence and applications[END_REF]. Here, as there, this difficulty is overcome once we note that we do not really care which fixed point we end up at, just that we get closer and closer to the set of fixed points. But we know this happens: we have a global convergence result (Theorem 3.1).

Theorem 3.4. Let

A = 0 B B T 0 ,
where B ∈ R n×n is nonnegative and fully indecomposable, and let Q be the doubly stochastic matrix one gets from diagonally balancing B.

The iteration (3.8) converges linearly for all initial vectors x 0 > 0 with asymptotic rate of convergence ρ = (1+σ 2 (Q))/2, where σ 2 (Q) is the second largest singular value of Q.

As with Theorem 3.3, we can extend Theorem 3.4 to cover unsymmetric matrices which are not fully indecomposable but have total support by working with a block diagonal permutation.

Proof of Theorem 3.4. Convergence of the algorithm for any starting vector is guaranteed by Theorem 3.1. Let us assume that the limit point is x * = [r T c T ] T . Note that the set of fixed points of (3.8) can be written as

S = βr 1 β c , β ∈ R + .
Let > 0 and choose K so that for k > K, x kx * < , where we choose (for reasons we will explain)

x * = x T D(x * ) -2 x 2n .
We will show that (3.12) min

s∈S x k+1 -s * ≤ ρ min s∈S x k -s * + τ k ,
where τ k = O( 2), and hence we can infer the desired convergence rate for small enough values of . Note that x * may not be the nearest element to x k in S: the nearest point is

s k = αr 1 α c
, where α = 1 + α and α = O( ). Now

x k+1 = f (x k ) = f (x * + p) = x * + J(x * )p + q,
where p * < and q * = O( 2 ). From (3.11) we can deduce that the eigenvalues of J(x * ) are

λ i = 1 2 (1 + σ i (Q)), λ i+n = 1 2 (1 -σ n+1-i (Q)), i = 1, 2, . . . , n,
where σ i (Q) is the ith singular value of Q. Since B is fully indecomposable, then so is Q, and it is primitive. Hence

1 = σ 1 (Q) > σ 2 (Q) ≥ • • • ≥ σ n (Q) ≥ 0.
Let v i be the eigenvector of J(x * ) associated with λ i such that v i * = 1 and write

p = 2n i=1 μ i v i .
From (3.11) we know that J(x * ) is similar to a symmetric matrix. Our choice of norm is motivated by our need to control the size of J(x * )p for all k > K. Since the v i form an orthonormal set with respect to the inner product that induces our norm,

|μ i | ≤ p * < for all i. Noting that λ 1 = 1, v 1 = [r T -c T ]
T , λ 2n = 0, and v 2n = x * , we have

x k+1 = x * + μ 1 v 1 + 2n-1 i=2 λ i μ i v i + q = (1 + μ 1 )r (1 -μ 1 )c + 2n-1 i=2 λ i μ i v i + q = γr c/γ + 0 ηc + 2n-1 i=2 λ i μ i v i + q,
where γ = (1 + μ 1 ) = 1 + O( ) and η = -μ 2 1 /(1 + μ 1 ).

In order to establish (3.12), we need to relate the size of the sum in our expression for x k+1 to min s∈S x ks * = x ks k * . Since

x k -s k = p + x * -s k = p + (1 -α)r (1 -1/α)c = 2n i=1 μ i v i + 1 -α 2 2α v 1 - (α -1) 2 2α v 2n , we have 2n-1 i=2 λ i μ i v i * ≤ λ 2 2n-1 i=2 μ i v i 2 * = λ 2 2n-1 i=2 μ i v i * ≤ λ 2 x k -s k * ,
and letting

τ k = 0 ηc + q * = O( 2 ),
we are done.

Our rate of convergence result allows a direct comparison with other scaling algorithms. In particular, it is shown in [START_REF] Knight | The Sinkhorn-Knopp algorithm: Convergence and applications[END_REF] that the Sinkhorn-Knopp algorithm also converges linearly with an asymptotic rate of convergence of σ 2 (Q). Note that since σ 2 (Q) < 1, we have ρ > σ 2 (Q). Hence the bound on the asymptotic rate of convergence (and empirical evidence suggests that it is generally a sharp bound) is necessarily larger than the corresponding bound for the Sinkhorn-Knopp algorithm.

However, for symmetric matrices the comparison is not so clear cut. The asymptotic rates are λ 2 (P) for the Sinkhorn-Knopp algorithm versus (1λ min (P))/2 for Algorithm 1. For sparse matrices we can expect both of these factors to be close to one (especially if the matrix is close to being reducible). But note that when |A| is symmetric positive definite, the rate of convergence of Algorithm 1 is guaranteed to be less than 1/2, since P is symmetric positive definite, too.

Often, though, we need only a modest degree of scaling to induce beneficial effects in linear solvers. Here, the fact that Algorithm 1 retains symmetry is a distinct advantage: even though a symmetric matrix balanced by the Sinkhorn-Knopp algorithm will (in the limit) be symmetric, the intermediate scalings lose this property.

A number of balancing algorithms [START_REF] Knight | A fast algorithm for matrix balancing[END_REF][START_REF] Livne | Scaling by binormalization[END_REF] have been developed based on Newton iterations for solving (3.9). While these retain symmetry in the intermediate scalings and can achieve quadratic convergence asymptotically, they can behave erratically when the iterates are far from being truly balanced, a phenomenon we have not seen with Algorithm 1. Typically, we have found that this erratic behavior occurs1 when the linear systems in the Newton step are singular (but consistent) and may be due to unwanted elements of the null-space being included in the early iterates. In future work, it may be worth investigating whether Newton-based methods can be employed in the applications we discuss in our numerical experiments. One limitation of these methods is that they are not designed for ∞-norm scaling.

Numerical experiments.

We have implemented the proposed algorithm in MATLAB (version R2009a), and compared it against our implementation of Bunch's and Sinkhorn-Knopp algorithms. We have experimented with a set of matrices obtained from the University of Florida Sparse Matrix Collection (UFL) available at http://www.cise.ufl.edu/research/sparse/matrices/. The matrices in our data set satisfy the following properties: real, 1000 ≤ n ≤ 121000, 2n ≤ nnz ≤ 1790000, without explicit zeros, fully indecomposable, not a matrix with nonzeros only in the set {-1, 1}. There were a total of 224 such matrices at the time of experimentation. In seven out of those 224 matrices, MATLAB's condest function failed to compute an estimated condition number in a reasonable amount of time, on three matrices condest gave inf as the estimated condition number, and one matrix was reported to be rank-deficient at the UFL collection. Therefore, we give results on 213 matrices, among which 64 are unsymmetric, 46 are symmetric positive definite (SPD), and 103 are symmetric but not positive definite.

We present two different sets of experiments. In the first one, we compare different scaling algorithms. In the second one, we investigate the merits of using the proposed scaling algorithm in solving linear systems with a direct method. In terms of convergence, there is little difference in theory which p-norm we choose to scale by in Algorithm 1 (apart from the ∞-norm). But to investigate the potential of the algorithm in practice we have tested both the 1-and 2-norm versions alongside the ∞-norm. For comparison, we have also tested the 1-and 2-norm versions of the Sinkhorn-Knopp algorithm as well as Bunch's ∞-norm algorithm. For an unsymmetric matrix A, we apply Bunch's algorithm to 0 A T A 0

which is equivalent to applying the Sinkhorn-Knopp algorithm using the ∞-norm (SK-inf). For symmetric matrices (both general and SPD), we did not include results with SK-inf, as this alternative loses the symmetry. For SPD matrices, we also use the reciprocal of the square root of diagonal entries as the scaling matrices (on both sides). We refer to this algorithm as rdiag.

Condition numbers and iterations.

We use the performance profiles [START_REF] Dolan | Benchmarking optimization software with performance profiles[END_REF] to assess the effect of the scaling algorithms on the condition number. In order to generate these profiles, we compute the condition numbers of a given matrix and its scaled versions (with different scaling algorithms). This gives us a set of condition estimates {χ 1 , . . . , χ p } corresponding to p different scaling algorithms (the condition number of the original matrix can be thought of as resulting from a scaling algorithm) for a given matrix. We then compute {χ 1 /χ, . . . , χ p /χ} where χ is the smallest value in the set. We take the logarithms of these ratios to avoid large numbers, and designate the resulting number as the performance of an algorithm on the given matrix. Then, we measure the percentage of the number of matrices in which the performance of an algorithm is less than or equal to a given number τ . Plotting these measurements gives the performance profiles shown in Figure 1 (the three subfigures correspond to the three classes of matrices). In other words, for a given value of τ , a plot in Figure 1 shows the fraction of the number of matrices for which the corresponding scaling algorithm obtains a condition estimate within e τ of the smallest (among all the condition estimates). A high fraction indicates that a given algorithm regularly keeps the estimated condition number under control and hence indicates that algorithm's effectiveness (particularly if the fraction is high for small values of τ ).

In these experiments, all the iterative scaling methods are run with an error tolerance of 1.0e-4 and with the number of iterations limited to the matrix dimension. All condition numbers are computed using MATLAB's condest function. For clarity, the profiles for Sinkhorn-Knopp algorithm in 1-and 2-norms are not shown. Generally, they agreed with the proposed algorithm in the corresponding norms to at least four significant digits-exceptions to this rule tend to coincide with cases where either the proposed algorithm or Sinkhorn-Knopp had not converged in n iterations. In all of the subfigures of Figure 1, the profiles of the scaling algorithms (with the single exception of rdiag) are higher than the profile representing no scaling (shown with A), for any range of τ . Therefore, one concludes that the scaling algorithms improve the estimated condition number in an overwhelming majority of cases. For unsymmetric matrices, the algorithms that work with the ∞-norm (inf, SK-inf, and B) outperform (or, have higher profiles than) the 1-and 2-norm scalings until τ > 3.2 after which the 2-norm scaling takes over. For general symmetric matrices, a clear pattern is seen in Figure 1b. This time, the 2-norm scaling has the highest performance profile, marginally beating the 1-norm scaling algorithm which in turn appears superior to the proposed algorithm in the ∞-norm over all shown ranges of τ . Bunch's algorithm is ranked most poorly. For the SPD matrices, all algorithms (except rdiag) are almost equal if one considers estimated condition numbers that are within e of the best, as the fraction of the test cases that are below τ = e are above 0.90.

We give further details in Table 1 which shows that all of the scaling algorithms (except rdiag) improve the condition number in most of the instances (below the row corresponding to 1.0e+0). In fact, rdiag increases the condition number of most The numbers of matrices in which the estimated condition number of the scaled matrix is larger than (the first half of the table) and smaller than (the second half ) the condition number (estimated) of the original matrix by differing amounts, for the unsymmetric, symmetric, and SPD matrices in the data set. Table 2 The statistics on the numbers of iterations of Sinkhorn-Knopp algorithm in 1-and 2-norms (SK-1 and SK-2) and the proposed algorithm in 1-, 2-, and∞-norms (A-1, A-2, andA-∞) to reach an error tolerance ε = 1.0e-4 with a limit of n in the number of iterations.

Unsymmetric (64) Symmetric (103) SPD (46) Ratio SK-inf B A-inf A-1 A-2 B A-inf A-1 A-2 √ d B A-inf A-1 A-2 ≥1.0e+3 0 0 0 3 0 0 0 0 0 21 0 0 0 0 ≥1.0e+2 0 0 0 3 1 1 0 1 0 25 0 0 0 0 ≥1.
Stat. n SK-1 SK-2 A-1 A-2
A-∞ its its/n its its/n its its/n its its/n its its/n unsymmetric ( 64 As seen from these numbers, using the proposed algorithm in 1-and 2-norms for unsymmetric matrices can cause ill-conditioning (problem ids at the UFL collection are 287, 289, 816, 1200, and 1433). In matrices which are close to being decomposable, scaling algorithms "want" to zero out certain elements. This cannot be done exactly and the resultant lop-sided scalings cause the condition number to go through the roof (up by a factor of 10 34 in extremis). We do not advocate scaling as a panacea but, particularly in the symmetric case (when, if there is a full diagonal, we can expect Sinkhorn-Knopp type algorithms to rapidly converge) the scaling is valuable. And note that the negative phenomenon of wildly worse conditioning is not present in ∞norm scalings for unsymmetric matrices as the pressure to remove entries disappears.

We present some statistical indicators (the minimum, the maximum, and the geometric mean) of the number of iterations with the same setting as before for the proposed algorithm in 1-, 2-, and ∞-norms and Sinkhorn-Knopp algorithm in 1-and 2-norms in Table 2. In this table, n is the size of the matrices, "its" is the number of iterations, and "its/n" is the ratio of these two quantities. For the latter column, we compute the ratio for each matrix and then take the minimum, the maximum, and the geometric mean. As seen from the numbers, the proposed algorithm in ∞-norm converges rapidly, displaying a geometric mean of iterations below ten. The 1-norm and 2-norm scaling algorithms have a large number of iterations (see geometric means) for unsymmetric matrices. Those numbers for the Sinkhorn-Knopp algorithm are even larger. The 2-norm variant of both algorithms have larger geometric means, and hence slower convergence than their 1-norm variants. For the symmetric and symmetric positive definite matrices, the number of iterations of the proposed algorithm in the 1-and 2-norms are better than the Sinkhorn-Knopp variants.

4.2.

Uses of the methods in a direct solver. We investigate the use of the proposed scaling algorithm in the solution of linear systems with direct methods. We do not perform experiments with SPD matrices. We use the solver MUMPS [START_REF] Amestoy | A fully asynchronous multifrontal solver using distributed dynamic scheduling[END_REF][START_REF] Amestoy | Hybrid scheduling for the parallel solution of linear systems[END_REF] version 4.10 in our tests with its default parameter settings except for the following. For the unsymmetric matrices, if there are zeros in the diagonal, we find a maximum traversal using dmperm of MATLAB and permute the matrix to put that transversal into the diagonal before calling MUMPS and set the related parameter ICNTL(6) to 0 in MUMPS. For the general symmetric matrices, we set the parameter ICNTL [START_REF] Bunch | Equilibration of symmetric matrices in the max-norm[END_REF] to 2 to use the compressed ordering strategies [START_REF] Duff | Strategies for scaling and pivoting for sparse symmetric indefinite problems[END_REF] along with a columnwise permutation option ICNTL(6)=1, equivalent in essence to dmperm of MATLAB. Before the numerical factorization, MUMPS allocates some memory with amount equal to the estimated one plus X% more where X is specified by the memory relaxation parameter ICNTL [START_REF] Curtis | On the automatic scaling of matrices for Gaussian elimination[END_REF]. MUMPS returns an error message if the actual memory requirements that arise during numerical factorization pass the allocated one. Extra allocation of memory is necessary when pivoting is used to ensure numerical stability. A scaling strategy is deemed to be better if it increases the number of successful runs (by reducing pivoting and hence its effects on memory and work). No other warnings or error messages were emitted by MUMPS for the matrices in the data set.

As the number of iterations of the proposed algorithm can be high (see Table 2), we recommend its use with a fixed and moderate number of iterations rather than to full convergence. In order to suggest a solid scaling approach, we have investigated different strategies specified below by three integers [i 1 , i 2 , i 3 ]. Each strategy proceeds in three phases: i 1 steps of ∞-norm scaling; then i 2 steps of 1-or 2-norm scaling; and finally i 3 steps of ∞-norm scaling. It is understood that each phase continues using the scaling factors found so far, and if convergence occurs within a phase, then the next phase is started. We use the ∞-norm scaling in the first and third phase as a smoother to reduce each nonzero entry below one. The 1-and 2-norm scaling iterations are two lead to different scaled matrices; therefore, we do not use one after another. The results are presented in Table 3. In this table, the scaling strategies are represented as [i 1 , i

(1) 2 , i 3 ] or [i 1 , i (2) 2 , i 3 ] where i (1) 2 and i (2)
2 designate the number of 1-norm or 2-norm iterations performed in phase 2 of the scaling strategy. We also include results with Bunch's algorithm, Sinkhorn-Knopp algorithm with ∞-norm (very close to that of Bunch's), and Sinkhorn-Knopp algorithm with ten iterations of 1-norm scaling. We give results with different levels of memory relaxation parameter. We note that MUMPS can run to completion for all matrices in our data set with one of the scaling strategies with ICNTL(14)= 40. Therefore, any error message returned by MUMPS signifies numerical difficulties.

We comment more on symmetric matrices, as the proposed algorithm is particularly motivated for this case. As seen in the right half of Table 3, the proposed algorithm requires only a few iterations to reduce the number of unsuccessful runs from 18 (with ICNTL(14)=10) and 9 (with ICNTL(14)=80). In general, increasing the number of iterations improves the success rate. The strategies with the 1-norm iterations are more effective than those with the 2-norm iterations (compare, for ex-Table 3 The effects of different scaling strategies on the direct solver MUMPS version 4.10. A strategy is said to be unsuccessful (usuc.), if MUMPS returns an error message with the settings described in the text, where usuc(X) gives the number of erroneous cases with ICNTL(14) = X. The geometric mean of the ratio of the condition number estimate of the scaled matrix to the original one is given in the column cnd(DAE) cnd(A) .

Unsymmetric matrices

General symmetric matrices Strategy usuc. ( 4.00e-02 [0, 3 (1) , 0] 2 0 3.05e-02 4 0 1.48e-02 [1, 3 (1) , 0] 2 0 3.13e-02 4 0 1.39e-02 [1, 3 (1) , 1] 2 0 3.25e-02 4 0 1.68e-02 [1, 3 (1) , 3] 2 0 3.63e-02 5 0 2.08e-02 [0, 3 (2) , 0] 3 0 3.87e-02 7 3 2.18e-02 [1, 3 (2) , 0] 2 0 3.57e-02 ample, the strategies with [•, 3 (1) , •] and [•, 3 (2) , •]). Most of the time, the condition number estimates are better with the scaling strategies having 1-norm iterations (we note that MATLAB's estimator works in the 1-norm, so this may also have some effect). Bunch's algorithm is also very effective for symmetric matrices in reducing the unsuccessful runs with ICNTL(14) = 80; but not that successful with ICNTL(14) = 10. For the unsymmetric matrices, the strategies with the 1-norm iterations are better than those with the 2-norm iterations in reducing the condition estimates and also marginally better in reducing the unsuccessful runs. Bunch's and SK∞ algorithm fail for almost the same number of matrices as with no-scaling (even one more with ICNTL(14) = 10). We conclude that the scaling algorithms (proposed and existing ones) are useful in reducing the need for pivoting so that memory requirements (and work) during numerical factorization are kept close to a priori estimates. In particular, strategies [0, 3 (1) , 0] or [1, 3 (1) , 0] are viable for symmetric and unsymmetric matrices, both in sequential and parallel computing environments. Equivalent SK iterations are also viable for unsymmetric matrices. For symmetric matrices, Bunch's algorithm has merit in sequential computing environments with relatively loose memory limits for the direct solver, but would be harder to parallelize.

7 3 1.97e-02 [1, 3 (2) , 1] 2 0 3.90e-02 7 2 2.27e-02 [1, 3 (2) , 3] 2 0 3.93e-02 9 2 2.61e-02 [0, 10 (1) , 0] 2 0 3.12e-02 3 0 1.71e-02 [1, 10 (1) , 0] 2 0 3.11e-02 3 0 1.69e-02 [1, 10 (1) , 1] 2 0 3.33e-02 3 0 1.88e-02 [0, 10 (2) , 0] 3 0 3.54e-02 5 1 1.57e-02 [1, 10 (2) , 0] 3 0 3.31e-02 4 1 1.63e-02 [1, 10 (2) , 1] 3 0 3.43e-02 4 1 1.74e-02 [1, 100 (1) , 0] 2 0 1.40e-03 3 0 1.08e-02 [1, 100 (2) , 0] 2 0 3.46e-02 3

Conclusion.

We presented an iterative algorithm which scales the ∞-norm of the rows and columns of a matrix to 1. The important features of the proposed algorithm are the following: it preserves symmetry; it is permutation independent;

and it has fast linear convergence with an asymptotic rate of 1/2. We discussed the extension to the 1-norm in detail. Again, the algorithm preserves symmetry and is permutation independent. From the various theorems collecting the convergence analysis results in the different norms, it can be seen that the assumptions are much more restrictive and complicated in their combinatorial aspects in the 1-norm as opposed to the ∞-norm. In the 1-norm, we have established that convergence depends on the nonzero structure of the matrix in much the same way as other Diagonal Product Increasing algorithms, such as Sinkhorn-Knopp; but that for symmetric matrices our algorithm can accelerate this convergence. This is confirmed in experimental results.

The rates of convergence in the particular case of our algorithm, fixed and fast in the ∞-norm, problem-dependent, and potentially much slower in the 1-norm, also illustrate the strong difference between the two. But whatever norm we choose, we have numerical evidence for the algorithm's potential to reduce the condition number of a matrix. Our numerical experiments demonstrate that the proposed algorithm can be helpful for a direct solver in the sense that it reduces the need for numerical pivoting. In particular, we have shown that one step of ∞-norm scaling followed by a few steps of 1-or 2-norm scaling is often good enough. We have also experimentally demonstrated that the proposed algorithm usually converges faster than Sinkhorn-Knopp scaling method for symmetric matrices.

Appendix A. General convergence results in the 1-norm. We extend partly the general convergence theorem of Parlett and Landis [START_REF] Parlett | Methods for scaling to double stochastic form[END_REF] regarding the convergence of Diagonal Product Increasing (DPI) algorithms. We show that some of their hypotheses can be weakened while keeping the same results, and that under some stronger hypothesis, their general convergence result can be strengthened, too. Our demonstrations here very closely follow the discussions by Parlett and Landis [24, Theorem 1, pp. 63-68], and improves only partially the already very general scope of their theoretical analysis. We have just included for completeness all the material about our extensions, since they cover directly the specific case of our Algorithm 1 in the 1-norm. We hope that the demonstrations here might be useful to prove the convergence of scaling alternatives.

For the notation, consider an iterative scaling algorithm producing the iterates (3.2) from section 3.2; that is, (k) , where

A (k) = a (k) ij = D (k) AE
D (k) = diag d (k) 1 , . . . , d (k) n and E (k) = diag e (k) 1 , . . . , e (k) n
, and where D (0) = E (0) = I. Let r (k) i for i = 1, . . . , n, and c (k) j for j = 1, . . . , n, denote the 1-norm of rows and columns, respectively. Assuming A ≥ 0, the row and column sums are given as r

(k) i = n j=1 a (k) ij and c (k) j = n i=1 a (k) ij .
We first restate the theorem of Parlett and Landis [START_REF] Parlett | Methods for scaling to double stochastic form[END_REF] characterizing the convergence of DPI algorithms in a slightly different form.

Theorem A.1. Suppose that a given scaling algorithm produces a sequence (3.2) of diagonal equivalents for a nonnegative matrix A that satisfies the following properties:

(P1) the sequence (s k ) k≥1 of the product of both the row and column scaling factors

s k = n i=1 d (k) i e (k) i , k = 1, 2, . . . ,
is monotonically increasing, (P2) the 1-norm of rows and columns, respectively, converge to 1:

lim k→∞ r (k) i = 1 and lim k→∞ c (k) j = 1, 1 ≤ i, j ≤ n.
Then, if A has support, A = k→∞ A (k) exists and is doubly stochastic. Additionally, if A has total support, then this limit is diagonally equivalent to A.

The property (P1) is exactly the same as in the hypothesis of Parlett and Landis. Property (P2) differs slightly in that it replaces and weakens partially the two other hypothesis made by Parlett and Landis, but is sufficient, combined with (P1), to raise the above conclusions (which are the same as in [START_REF] Parlett | Methods for scaling to double stochastic form[END_REF]Theorem 1]).

We first establish the following intermediate result.

Lemma A.2. Using the same notation as in Theorem A.1 above, property (P1), with the added assumption that all the elements a (k) ij are bounded above by some constant independent of k ≥ 1, as well as the fact that A has support, imply the following:

1. lim k→∞ s k = ξ > 0 exists, 2. there exists a strictly positive constant γ such that

d (k) i e (k) j ≥ γ
for all k ≥ 1 and for each index pair (i, j) such that a ij is strictly positive and can be permuted into a positive diagonal. Proof. Since A has support, there exists a permutation σ such that a

i,σ(i) > 0, 1 ≤ i ≤ n. Let a = min 1≤i≤n (a i,σ(i) ). Then, for all k ≥ 1, n i=1 d (k) i e (k) σ(i) a ≤ n i=1 d (k) i e (k) σ(i) a i,σ(i) = n i=1 a (k) i,σ(i) ≤ β
for some constant β > 0 (the last inequality resulting from the fact that all the elements are bounded above). Then, by the arithmetic-geometric mean inequality,

∀k ≥ 1, s k = n i=1 d (k) i e (k) σ(i) ≤ 1 n n i=1 d (k) i e (k) σ(i) n ≤ β n a n ,
and, combined with (P1), the monotonically increasing sequence (s k ) k≥1 is thus bounded. Consequently, lim k→∞ s k = ξ > 0 exists and is finite. The demonstration of the second point goes around the same ideas. For any a ij > 0 that can be permuted into a positive diagonal, there exists a permutation σ such that j = σ(i) and a l,σ(l) > 0, for 1 ≤ l ≤ n. Then, by (P1),

∀k ≥ 1, d (k) i e (k) j n =1, =i d (k) e (k) σ( ) = s k ≥ s 1 and ∀k ≥ 1, d (k) i e (k) j ≥ s 1 n =1, =i d (k) e (k) σ( ) -1 . Let a = min 1≤i,j≤n (a ij > 0). Then, for all k ≥ 1, n =1, =i d (k) e (k) σ( ) a ≤ n =1, =i d (k) e (k) σ( ) a ,σ( ) = n =1, =i a (k) ,σ( ) ≤ β,
and, by the arithmetic-geometric mean inequality, we can conclude that ∀k ≥ 1, d

(k) i e (k) j ≥ s 1 (n-1)a β n-1 = γ > 0 .
We are now ready to prove Theorem A.1. Proof of Theorem A.1. As a direct consequence from (P2), the sequence of matrices (A (k) ) k≥1 is bounded in the finite dimensional space of real, n × n matrices. Consider any convergent subsequence (A (q) ) q of (A (k) ) k≥1 , and define

A = (â ij ) = lim q→∞ A (q) . From (P2), lim q→∞ r (q) i = 1 and lim q→∞ c (q) j = 1, 1 ≤ i, j ≤ n,
which implies that A is doubly stochastic. Then, since the set of n × n doubly stochastic matrices is the convex hull of the set of n × n permutation matrices (see [START_REF] Birkhoff | Tres observaciones sobre el algebra lineal[END_REF]), A must thus have total support. Therefore, âij = lim q→∞ a (q) ij = 0 whenever a ij > 0 cannot be permuted onto a positive diagonal. Consider any entry a ij > 0 in A that can be permuted onto a positive diagonal, and let

μ ij = âij a ij = lim q→∞ d (q) i e (q) j .
From point 2 in Lemma A.2, we know that μ ij ≥ γ > 0.

Then, applying Lemma 2 of [START_REF] Parlett | Methods for scaling to double stochastic form[END_REF] (which is itself paraphrased from [31, p. 345]), we know that there exist positive sequences (x (q) i ) q and (y (q) i ) q, both with positive limits such that

d (q) i e (q) j = x (q) i y (q) j
for all âij > 0 in A, and for all q ≥ 1.

Then, taking (q) and Y = lim q→∞ Y (q) , we have

X (q) = diag x (q) 1 , . . . , x (q) n , Y (q) = diag y (q) 1 , . . . , y (q) n , X = lim q→∞ X
A = lim q→∞ D (q) SE (q) = lim q→∞ X (q) SY (q) = XSY,
showing that A is diagonally equivalent to S, where S is the submatrix with total support extracted from A where all entries a ij > 0 in A that cannot be permuted into a positive diagonal have been set to zero, the others remaining the same. Finally, if we consider any other convergent subsequence in (A (k) ) k≥1 , for the same reasons as above its limit will also be doubly stochastic and diagonally equivalent to S. Since doubly stochastic equivalents are unique (see [START_REF] Sinkhorn | Problems involving diagonal products in nonnegative matrices[END_REF]), the two limits must be the same. Therefore, we can conclude that lim k→∞ A (k) exists and is doubly stochastic. If, additionally, A has total support, then A = S and is directly diagonally equivalent to the doubly stochastic matrix A, which completes the proof of Theorem A.1.

The previous two properties (P1) and (P2) are insufficient to prove that the two sequences of scaling matrices D (k) and E (k) actually converge. The reason why this additional result holds for the proposed Algorithm 1 is that, not only is the sequence (s k ) k≥1 monotonically increasing, but also both sequences (

n i=1 d (k) i ) k≥1 and ( n i=1 e (k)
i ) k≥1 are independently monotonically increasing. We now extend the conclusions in Theorem A.1 by incorporating this property (denoted by (P1 * ) in section 3.2 and repeated below). Note that (P1 * ) implies property (P1).

Theorem A. (k) and E = lim k→∞ E (k) do exist and A = DAE is doubly stochastic.

Proof. Since (P1 * ) implies (P1), and since (P2) remains the same, we know from Theorem A.1 that if A has total support, A = lim k→∞ A (k) exists and is diagonally equivalent to A. To prove that both D = lim k→∞ D (k) and E = lim k→∞ E (k) exist and A = DAE, we exploit the assumption that A is fully indecomposable. In this case, A is diagonally equivalent to the doubly stochastic limit A = lim k→∞ A (k) , and the diagonal matrices which take place in this equivalence are unique up to a scalar factor [START_REF] Sinkhorn | Concerning nonnegative matrices and doubly stochastic matrices[END_REF].

Let us suppose now that one of the sequences (d

(k)
i ) k≥1 is unbounded for some i. In such a case, there exist a subsequence (d (q) i ) q such that lim q→∞ d (q) i = +∞.

As the matrix A is fully indecomposable, for any index j, 1 ≤ j ≤ n, there exist a chain of positive entries in which the row and column indexes alternately change a ij1 , a i1j1 , a i1j2 , a i2j2 , a i2j3 , . . . , a ip-1jp-1 , a ip-1jp , a ipjp , a ipj "connecting" row i to column j (see, for example, [START_REF] Brualdi | Combinatorial Matrix Theory[END_REF]Theorem 4.2.7]). Consequently, because of the conclusions raised at point 2 in Lemma A.2, we know that the denominator in the fraction Since this can be done for any j, we can conclude that the subsequence ( n j=1 e (q) j ) q (which is strictly positive) goes to zero as (d (q) i ) -n n j=1 μ ij . This last conclusion contradicts with property (P1 * ) above, stating that both sequences ( (k) i k≥1 are monotonically increasing, it is clear that α must be equal to 1 and the two limits must be equal. Therefore, we can conclude that D = lim k→∞ D (k) and E = lim k→∞ E (k) exist, which completes the proof.

The matrices with total support can be permuted into a direct sum of fully indecomposable matrices [START_REF] Perfect | The distribution of positive elements in doubly-stochastic matrices[END_REF]Theorem 1]. Now, if the algorithm generates iterates (3.2) that are insensitive to permutations of A, and if applying it to a block diagonal matrix is equivalent to working directly which each diagonal block separately, we can easily extend the conclusions of Theorem A.3 to the case of matrices with total support by collecting the conclusions from Theorem A.3 independently on each fully indecomposable submatrix in such a direct sum. This is presented as the following corollary.

Corollary A.4. Suppose that a given scaling algorithm produces a sequence of diagonal equivalents, as in (3.2), for a nonnegative matrix A that satisfy the two properties (P1 * ) and (P2), given in Theorem A.3 above, and if, additionally, 1. the algorithm is permutation insensitive, in the sense that, under any row or column permutation of the original matrix A, the scaling elements remain the same and are just permuted along the diagonals of D (k) and E (k) with respect to the corresponding row or column permutations applied to A; 2. the iterates also remain the same when applying the algorithm to a block diagonal matrix as when collecting results obtained by applying the algorithm on each diagonal submatrix separately. Then, if A has total support, both limits D = lim k→∞ D (k) and E = lim k→∞ E (k) do exist and A = DAE is doubly stochastic.

The two extra hypothesis of Corollary A.4 are naturally satisfied by algorithms that generate iterates by relying only on the computation of the one norm of both rows and columns at each iteration (as is the case of the proposed and Sinkhorn-Knopp algorithms). We recall that Parlett and Landis use the fact that their algorithm incorporates a normalization step, with

μ k = 1 n n i=1 r (k+1) i = 1
enforcing the mean of all row sums to be set to 1 at each iteration. This was actually incorporated as property (P3) in the hypothesis of their main convergence theorem.

As we have seen above, this is not mandatory to reach the same conclusions as in Theorem A.1. Additionally, this normalization step may, however, violate the second hypothesis in Corollary A.4, because the averaging of row sums may not be the same when computed globally or blockwise in the case of a block-diagonal matrix. This may not prevent anyway the algorithms proposed in [START_REF] Parlett | Methods for scaling to double stochastic form[END_REF] (they incorporate this normalization step) to produce still a convergent sequence of scaling factors. We have not found how to characterize in a simple way some common generic properties, that might enable to reach the same conclusions in the case of matrices with total support only, and extend the scope of Corollary A.4 to address these particular "normalized " scaling algorithms.

2 .

 2 The algorithm. Consider a general m × n real matrix A. For i = 1, . . . , m, let r i = a T i * ∈ R n×1 denote the ith row of A as a vector, and for j = 1, . . . , n, let c j = a * j ∈ R n×1 denote the jth column of A. Let R and C denote the m × m and n × n diagonal matrices given by (2.1) R = diag r i ∞ i=1,...,m and C = diag c j ∞ j=1,...,n

∞

  , respectively, each one of them being greater than or equal to |a ij | itself.For any subsequent iteration (k ≥ 1), consider the ∞-norm of any row r let indices and p satisfy the equalities |a

i

  (P1 * ) both sequences of rows and column scaling factors products , k = 1, 2, . . . are monotonically increasing; (P2) the 1-norm of rows and columns, respectively, converge to 1, i.e.,

  ij by construction. Proof of Theorem 3.1. By (3.5) and the fact that the iterates in Algorithm 1 satisfy

Fig. 1 .

 1 Fig.1. Performance profiles for the condition number estimates of unsymmetric matrices 1a, general symmetric matrices 1b, and SPD matrices 1c. Horizontal axis in each subfigure corresponds to τ . A marks the condition number estimate of the original matrix; B marks that of Bunch's algorithm; inf, 1, and 2 mark that of the proposed algorithm with ∞-, 1-, and 2-norms, SK-inf marks that of Sinkhorn-Knopp with ∞-norm and rdiag marks that of the reciprocal square root of the diagonal scaling for SPD matrices.

  d

  from zero and. From the demonstration of Theorem A.1 above this fraction has a strictly positive limitlim k→∞ d (k) i e (k) j = μ ij1 μ i1j2 μ i2j3 . . . μ ip-2jp-1 μ ip-1jp μ ipj μ i1j1 μ i2j2 . . . μ ip-1jp-1 μ ipjp = μ ij > 0.

j

  ) k≥1 are monotonically increasing. The same could be done with one of the (e (k) j ) k≥1 instead, and we can finally conclude that each of the sequences (d(k) i ) k≥1and (e (k) i ) k≥1 are bounded. Now, since the two sequences (D (k) ) k≥1 and (E (k) ) k≥1 are bounded in the finite dimensional space of real n × n diagonal matrices, they have convergent subsequences. Let us consider two convergent subsequences:D (p) , E (p) -→ p→+∞ (D 1 , E 1 )andD (q) , E (q) -→ q→+∞ (D 2 , E 2 ) .Obviously, D 1 AE 1 = A = D 2 AE 2 , and thus, because of the uniqueness (shown in[START_REF] Sinkhorn | Concerning nonnegative matrices and doubly stochastic matrices[END_REF]), there exists α > 0 such that D 1 = αD 2 and E 1 = (1/α)E 2 . Then, as mentioned above, since both sequences

Table 1

 1

  being the same as in Theorem A.1 above. Then, if A is fully indecomposable, both limits D = lim k→∞ D

	3. Suppose that a given scaling algorithm produces a sequence (3.2)
	of diagonal equivalents for a nonnegative matrix A, each satisfying the following prop-
	erties:
	(P1 * ) the two products of row and column scaling factors
	n	n
	d (k) i	and	e	(k) i , k = 1, 2, . . . ,
	i=1	i=1
	are monotonically increasing;
	(P2)

Seen for some matrices from the GHS indef group from the UFL collection.

Acknowledgments and availability. We thank Patrick R. Amestoy, Iain S. Duff, Nick Gould, and Mario Arioli for commenting on an earlier draft. We also wish to thank the referees for their detailed reports that helped us improve the manuscript, and Stanley Eisenstat for his very constructive editorial and technical comments.

The MATLAB implementations of all scaling algorithms discussed in this paper and an MPI-based parallel implementation of the proposed scaling algorithm are available at http://perso.ens-lyon.fr/∼bucar/codes.html. The MPI-based parallel version is incorporated in MUMPS [2,4] (since version 4.8) for scaling distributed input matrices. The sequential Fortran implementation of the proposed scaling algorithm is available as MC77 in the HSL Mathematical Software Library.

This work was supported in part by "Agence Nationale de la Recherche" through SOLSTICE project (ANR-06-CIS6-010). http://www.siam.org/journals/simax/35-3/82575.html