Konstantin E Avrachenkov
email: k.avrachenkov@sophia.inria.fr

Jean B Lasserre
email: lasserre@laas.fr

Analytic Perturbation of Generalized Inverses

We investigate the analytic perturbation of generalized inverses. Firstly we analyze the analytic perturbation of the Drazin generalized inverse (also known as reduced resolvent in operator theory). Our approach is based on spectral theory of linear operators as well as on a new notion of group reduced resolvent. It allows us to treat regular and singular perturbations in a unified framework. We provide an algorithm for computing the coefficients of the Laurent series of the perturbed Drazin generalized inverse. In particular, the regular part coefficients can be efficiently calculated by recursive formulae. Finally we apply the obtained results to the perturbation analysis of the Moore-Penrose generalized inverse in the real domain.

Introduction

We study generalized inverses of analytically perturbed matrices:

A(z) = A 0 + zA 1 + z 2 A 2 + ... . (1)
First, we provide the perturbation analysis for the Drazin generalized inverse. In such a case we assume that the matrices A k , k = 0, 1, ..., are square of dimension n with complex entries. Furthermore, since we are interested in the perturbation analysis of the generalized inverse, we assume that the null space of the perturbed matrix A(z) is non-trivial.

Here we distinguish between regular and singular perturbations. The perturbation is said to be regular if it does not change the dimension of null space. Otherwise, the perturbation is said to be singular. The regular perturbation is also called rank-preserving perturbation. One of the main advantages of the complex analytic approach employed in the present work is that it allows us to treat both regular and singular perturbations in a unified framework.

If the coefficient matrices A k , k = 0, 1, ..., are real and we restrict ourselves to real z, the perturbation analysis of the Drazin generalized inverse can be applied to the perturbation analysis of the Moore-Penrose generalized inverse.

The main goals of the present work is to prove the existence of the Laurent series expansion for the perturbed Drazin generalized inverse

A # (z) = +∞ j=-s z j H j , (2)
and to provide a method for efficient computation of coefficients H j , j = -s, -s + 1, ... . We derive recursive formulae for the (matrix) coefficients of the regular part of the Laurent series expansion [START_REF] Avrachenkov | Inversion of analytic matrix functions that are singularly at the origin[END_REF]. Besides their theoretical interest, the recursive formulae are particularly 1 useful when one needs to compute a significant number of terms in the Laurent series expansion. These formulae require knowledge of the singular part coefficients and the latter are obtained via a reduction technique based on the complex analytic approach. In particular, this reduction technique uses a new notion of group reduced resolvent. The order of the pole of ( 2) is obtained as a by-product.

Last but not least, the limit matrix in the Taylor series expansion of the 0-eigenprojection matrix, has a simple exression in terms of (a) the original unperturbed 0-eigenprojection, (b) the perturbation matrices and (c), the singular part of the Laurent series [START_REF] Avrachenkov | Inversion of analytic matrix functions that are singularly at the origin[END_REF]. This provides some insight on how the perturbed 0-eigenvectors relate to the original 0-eigenvectors.

Let us also mention an important issue. In the course of the procedure that we propose, one has to determine the multiplicity of zero eigenvalues and verify whether they are semisimple. Of course, in the general situation neither task is easy in the presence of rounding errors. Nevertheless, we note that in many particular applications this issue can be effectively resolved. For instance, for the perturbation of the Moore-Penrose generalized inverse the semi-simplicity assumption is satisfied automatically, as we transform the problem into an equivalent problem with symmetric matrices.

To our best knowledge, Stewart [START_REF] Stewart | On the perturbation of pseudo-inverses, projections and linear least squares problems[END_REF] was the first to carry out the perturbation analysis of the Moore-Penrose generalized inverse. Even though he has mostly dealt with the regular perturbation case, he has pointed out that the error bounds grow infinitely as the perturbation vanishes in the singular case. He has provided error bounds and a first order approximation for the general componentwise perturbation of the Moore-Penrose generalized inverse. In fact, up to the present most of the work on perturbation analysis of generalized inverses is concentrated on regular componentwise perturbation. An interested reader can find results on regular componentwise perturbation of generalized inverses in [START_REF] Castro-González | Additive perturbations results for the Drazin inverse[END_REF][START_REF] Castro-González | On the Drazin inverse of the sum of two operators and its application to operator matrices[END_REF][START_REF] Wang | Generalized Inverses: Theory and Computations[END_REF] and references therein. The analytic perturbation is more restrictive than the general componentwise perturbation. However, it allows us to perform a deeper analysis of the singular case. The case of singular perturbations is interesting for it provides information on the lack of continuity of solutions to perturbed linear systems. The work on singular perturbation of the generalized inverse is very limited. Deif [START_REF] Deif | Sensitivity analysis in linear systems[END_REF] proposed an algorithm for computing the Laurent series for the generalized inverse of an analytically perturbed matrix. However, to implement his algorithm one has to obtain power series expansions for all the quantities in the singular value decomposition, namely for all the singular eigenvectors and eigenvalues. This can be very costly, especially for matrices of high dimension. In contrast, we propose an algorithm, which works in some reduced subspaces of small dimension.

We note that if matrix valued function A(z) is invertible in some punctured neibourhood around z = 0, then its generalized Drazin inverse is just the ordinary inverse. Therefore, our results are applicable to the classical problem of the inversion of nearly singular operators [START_REF] Avrachenkov | Inversion of analytic matrix functions that are singularly at the origin[END_REF][START_REF] Bart | Meromorphic operator valued functions[END_REF][START_REF] Gohberg | On the local theory of regular analytic matrix functions[END_REF][START_REF] Korolyuk | Mathematical foundations of the state lumping of large systems[END_REF][START_REF] Lancaster | Inversion of lambda-matrices and application to the theory of linear vibrations[END_REF][START_REF] Langenhop | The Laurent expansion for a nearly singular matrix[END_REF][START_REF] Pervozvanski | Theory of Suboptimal Decisions[END_REF][START_REF] Schweitzer | The Laurent expansion of pencils that are singular at the origin[END_REF][START_REF] Vainberg | Theory of branching of solutions of non-linear equations[END_REF][START_REF] Vishik | The solution of some perturbation problems in the case of matrices and self-adjoint and non-self-adjoint differential equations[END_REF], a problem with many important practical applications. In [START_REF] Avrachenkov | Perturbation of null spaces with applications to the eigenvalue problem and generalized inverses[END_REF] the perturbation analysis of the group inverse has been provided. Since the group inverse is a particular case of the Drazin inverse when the space can be decomposed into a direct sum of the null space and the range space, the results of the present work can be applied to the perturbation analysis of the group inverses. However, since the Drazin inverse is a more general notion than the group inverse, the results of [START_REF] Avrachenkov | Perturbation of null spaces with applications to the eigenvalue problem and generalized inverses[END_REF] cannot be applied to the perturbation analysis of the Drazin inverse.

The results of the present paper can be immediately applied to singular perturbations of Markov chains [START_REF] Avrachenkov | The fundamental matrix of singularly perturbed Markov chains[END_REF][START_REF] Delebecque | A reduction process for perturbed Markov chain[END_REF][START_REF] Haviv | Series expansions for stochastic matrices[END_REF][START_REF] Schweitzer | Perturbation series expansions for nearly completely-decomposable Markov chains[END_REF]. For example, the deviation matrix of a Markov chain is just (with minus sign) the reduced resolvent of the Markov chain generator. To demonstrate this application and verify the validity of our theoretical results, a simple illustrative example is provided at the end of the paper. Another application is concerned with perturbations of spectral inverses arising in the chemical networks, as treated in the work of Bohl and Lancaster [START_REF] Bohl | Perturbation of spectral inverses applied to a boundary layer phenomenon arising in chemical networks[END_REF]. Yet another interesting potential application of our results is the perturbation analysis of least square problems [START_REF] Deif | Sensitivity analysis in linear systems[END_REF][START_REF] Stewart | On the perturbation of pseudo-inverses, projections and linear least squares problems[END_REF].

Some preliminary results of the present work have appeared in the research report [START_REF] Avrachenkov | Perturbation analysis of reduced resolvents and generalized inverses[END_REF]. However, there are very important differences between the contents of the research report [START_REF] Avrachenkov | Perturbation analysis of reduced resolvents and generalized inverses[END_REF] and the present article. In the research report the proofs were carried out under the assumption that coefficient matrices of the analytic perturbation are real (even though the perturbation parameter z is complex). Here we have extended the results, whenever possible, to the case of complex coefficient matrices. Also the presentation of the proofs has been significantly improved. There is a new, very important, observation, that in general the Moore-Penrose generalized matrix of the complex perturbation is not analytic. This was absent in the research report.

The paper is organized as follows. In the next Section 2, we provide preliminaries on the complex analytic approach and spectral theory. In Section 3, we prove the existence of a Laurent series expansion for the Drazin generalized inverse of a perturbed matrix. In Section 4, we derive a recursive formula for the coefficients of the regular part of the Laurent series. For clarity of exposition, the proof of the main result is postponed to the Appendix. In Section 5, we propose a reduction process that permits to compute the coefficients of the singular part. In Section 6 we apply the obtained results to the perturbation analysis of the Moore-Penrose and group generalized inverses. Finally, in Section 7, we derive a simple expression for the limiting 0-eigenprojection matrix as the perturbation vanishes.

Preliminaries on complex analytic approach and spectral theory

Let us recall some facts from complex analysis and spectral theory. The book of Kato [START_REF] Kato | Perturbation theory for linear operators[END_REF] is an excellent source of the material on the subject. Any matrix A ∈ C n×n possesses the following spectral representation

A = p i=0 (λ i P i + D i ), (3)
where p + 1 is the number of distinct eigenvalues of A, P i is the eigenprojection and D i is the nilpotent operator corresponding to the eigenvalue λ i . By convention, λ 0 is the zero eigenvalue of A, that is, λ 0 = 0. In the case when there is no zero eigenvalue, the eigenvalues are enumerated from i = 1. The resolvent is another very important object in spectral theory.

Definition 1

The following operator-valued function of the complex parameter ζ is called the resolvent of the operator

A ∈ C n×n R(ζ) = (A -ζI) -1 .
The resolvent satisfies the resolvent identity:

R(ζ 1 ) -R(ζ 2 ) = (ζ 1 -ζ 2 )R(ζ 1 )R(ζ 2 ), (4)
for all ζ 1 , ζ 2 ∈ C. The resolvent has singularities at the points ζ = λ k , where λ k are the eigenvalues of A. In a neighbourhood of each singular point λ k the resolvent can be expanded as a Laurent series

R(ζ) = - m k -1 n=1 1 (ζ -λ k ) n+1 D n k - 1 ζ -λ k P k + ∞ n=0 (ζ -λ k ) n S n+1 k , (5)
where S k is the reduced resolvent corresponding to the eigenvalue λ k with geometric multiplicity m k . In fact, S k is the Drazin generalized inverse of (A -λ k I). And, in particular, we have S 0 = A # . The Drazin generalized inverse has the following basic properties

AA # = I -P 0 , (6)
P 0 A # = 0. (7)
The above equations show that A # is the "inverse" of A in the complementary subspace to the generalized null space of A. Here by generalized null space we mean a subspace which is spanned by all eigenvectors and generalized (Jordan) eigenvectors corresponding to the zero eigenvalue. Note that P 0 is a projection onto this generalized null space. Moreover, if the underlying space admits a decomposition into the direct sum of the null space and the range of the operator A (recall from [START_REF] Campbell | Generalizes Inverses of Linear Transformation[END_REF] that this is a necessary and sufficient condition for the existence of the group inverse), then the Drazin inverse and the group inverse coincide, and the following Laurent expansion holds

R(ζ) = - 1 ζ P 0 + ∞ n=0 ζ n (A # ) n+1 . ( 8
)
Since the Drazin generalized inverse is the constant term in the Laurent series ( 5) at z = λ 0 , it can be calculated via the following Cauchy integral formula

A # = 1 2πi Γ 0 1 ζ R(ζ) dζ, (9)
where Γ 0 is a closed positively-oriented contour in the complex plane, enclosing 0 but no other eigenvalue of A [START_REF] Avrachenkov | Analytic perturbation theory and its applications[END_REF]. The above formula will play a crucial role in the sequel. The Drazin inverse also has a simple expression in terms of eigenprojections, eigenvalues and nilpotent operators of the original operator A. Namely,

A # = p i=1    1 λ i P i + m i -1 j=1 (-1) j 1 λ j+1 i D j i    . ( 10
)
We emphasize that the above sum is taken over all indecies corresponding to non-zero eigenvalues. This expression again demonstrates that the Drazin generalized inverse is the inverse operator in the complementary subspace to the generalized null space. Moreover, this expression exactly represents the inverse operator A -1 whenever A has no zero eigenvalue.

Existence of the Laurent series expansion

In this section we prove the existence of a Laurent series expansion (2) for the Drazin generalized inverse A # (z) of the analytically perturbed matrix A(z). First let us consider the resolvent R(ζ, z) := (A(z) -ζI) -1 of the perturbed A(z). One can expand R(ζ, z) in a power series with respect to the complex variable z near z = z 0 , as follows (see e.g. [START_REF] Kato | Perturbation theory for linear operators[END_REF])

R(ζ, z) = R(ζ, z 0 ) + ∞ n=1 (z -z 0 ) n R (n) (ζ, z 0 ), (11)
where

R (n) (ζ, z 0 ) := ν 1 +...+νp=n (-1) p R(ζ, z 0 )A ν 1 R(ζ, z 0 )A ν 2 • • • R(ζ, z 0 )A νp R(ζ, z 0 ),
where A ν k are the coefficients of A(z) and ν k ≥ 1. The above expansion is called the second Neumann series for the resolvent. It is uniformly convergent for z sufficiently close to z 0 and ζ ∈ D, where D is a compact subset of the complex plane which does not contain the eigenvalues of A(z 0 ) [START_REF] Kato | Perturbation theory for linear operators[END_REF].

Theorem 1 Let A(z) be the analytic perturbation of the matrix A 0 given by (1). Then, the Drazin genralized inverse A # (z) of the perturbed operator A(z) can be expanded as a Laurent series (2).

Proof: We first show that there exists a domain 0 < |z| < z max such that A # (z) can be expanded in Taylor series at any point z 0 in this domain. For a fixed, arbitrary z > 0, (9) becomes

A # (z) = 1 2πi Γ 0 (z) 1 ζ R(ζ, z)dζ, (12)
where Γ 0 (z) is a closed positively-oriented curve enclosing the origin but no other eigenvalue of A(z).

With z max less than the modulus of any non-zero eigenvalue of A 0 , expand the perturbed resolvent in the power series [START_REF] Castro-González | On the Drazin inverse of the sum of two operators and its application to operator matrices[END_REF], around the point z 0 (with 0 < |z 0 | < z max ). Then, the substitution of that series in the integral formula [START_REF] Coderh | Hierarchical aggregation of linear systems with multiple time scales[END_REF], yields

A # (z) = 1 2πi Γ 0 (z 0 ) 1 ζ R(ζ, z 0 ) + ∞ n=1 (z -z 0 ) n R (n) (ζ, z 0 ) dζ.
Since the power series for R(ζ, z) is uniformly convergent for z sufficiently close to z 0 , we can integrate the above series term by term,

A # (z) = 1 2πi Γ 0 (z 0 ) 1 ζ R(ζ, z 0 ) dζ + ∞ n=1 (z -z 0 ) n 1 2πi Γ 0 (z 0 ) 1 ζ [R (n) (ζ, z 0 )] dζ = A # (z 0 ) + ∞ n=1 (z -z 0 ) n H n (z 0 ), (13)
where the coefficients are defined by

H n (z 0 ) := 1 2πi Γ 0 (z 0 ) 1 ζ [R (n) (ζ, z 0 )] dζ.
The convergence of power series [START_REF] Deif | Sensitivity analysis in linear systems[END_REF] in some non-empty domain 0 < |z| < z max can be shown by using the bounds for the contour integrals (for very similar development see for example [19, Ch.2, Sec.3]). From the power series (13), we can see that A # (z) is holomorphic in the domain 0 < |z| < z max . Consequently, by Laurent's theorem (see, e.g., [23, v. II, p. 7]), we conclude that A # (z) possesses a Laurent series expansion at z = 0 (with radius of convergence z max ), i.e.,

A # (z) = +∞ n=-∞ z n H n . (14)
We next show that the pole at z = 0 can be at most of finite order. Consider the spectral representation [START_REF] Castro-González | Additive perturbations results for the Drazin inverse[END_REF] for the reduced resolvent of the perturbed operator A(z)

A # (z) = p i=1    1 λ i (z) P i (z) + m i -1 j=1 (-1) j 1 λ j+1 i (z) D j i (z)    .
From Kato [START_REF] Kato | Perturbation theory for linear operators[END_REF], we know that the perturbed eigenvalues λ i (z) are bounded in |z| ≤ z max and they have at most algebraic singularities. Furthermore, the eigenprojections P i (z) and nilpotents D i (z) can also have only algebraic singularities and poles of finite order. Therefore, none of the functions λ i (z), P i (z) and D i (z) can have an essential singularity. This latter fact implies that their finite sums, products or divisions as in A # (z), do not have an essential singularity as well and, consequently, the order of pole in ( 14) is finite. This completes the proof. ✷

4 Recursive formula for the regular part coefficients

Here we derive recursive formulae for the coefficients of the regular part of the Laurent series (2). We use an analytic technique based on Cauchy contour integrals and resolvent-like identities.

First, observe that the structure of the perturbed Drazin inverse (A 0 + zA 1 + z 2 A 2 + ...) # is similar to the structure of the classical resolvent (A 0 -zI) -1 . Moreover, A # (z) becomes precisely the resolvent if A 1 = -I and A k = 0 for k ≥ 2. Therefore, one can expect that these two mathematical objects have some similar features. It turns out that the Drazin inverse of an analytically perturbed matrix A(z) satisfies an identity similar to the resolvent identity (4).

Lemma 1

The reduced resolvent A # (z) of the analytically perturbed operator A(z) = ∞ k=0 z k A k satisfies the resolvent-like identity:

A # (z 1 ) -A # (z 2 ) = ∞ k=1 (z k 2 -z k 1 )A # (z 1 )A k A # (z 2 ) + A # (z 1 )P 0 (z 2 ) -P 0 (z 1 )A # (z 2 ), ( 15
)
where P 0 (z) is the eigenprojection matrix corresponding to the zero eigenvalue.

Proof: Consider the following expression

A(z 2 ) -A(z 1 ) = ∞ k=1 (z k 2 -z k 1 )A k .
Premultiplying by A # (z 1 ) and postmultiplying by A # (z 2 ), yields

A # (z 1 )A(z 2 )A # (z 2 ) -A # (z 1 )A(z 1 )A # (z 2 ) = ∞ k=1 (z k 2 -z k 1 )A # (z 1 )A k A # (z 2 ).
Then, using (6), we get

A # (z 1 )[I -P 0 (z 2 )] -[I -P 0 (z 1 )]A # (z 2 ) = ∞ k=1 (z k 2 -z k 1 )A # (z 1 )A k A # (z 2 ).
Equivalently,

A # (z 1 ) -A # (z 2 ) = ∞ k=1 (z k 2 -z k 1 )A # (z 1 )A k A # (z 2 ) + A # (z 1 )P 0 (z 2 ) -P 0 (z 1 )A # (z 2 ),
which is the desired identity [START_REF] Gohberg | On the local theory of regular analytic matrix functions[END_REF]. ✷ In the next theorem, we obtain a general relation between the coefficients of the Laurent series (2).

Theorem 2 Let H k , k = -s, -s + 1, ... be the coefficients of the Laurent series (2) and P 0 (z) = ∞ k=0 z k P 0k be a power series for the eigenprojection corresponding to the zero eigenvalue of the perturbed operator. Then the coefficients H k , k = -s, -s + 1, ... satisfy the following relation

∞ k=1 k-1 i=0 H n-i A k H m+i-k+1 = -(η n + η m -1)H n+m+1 - 0, m < 0, 1 2πi Γ 1 z -n-1 1 A # (z 1 )[P 0m+1 + z 1 P 0m+2 + ...]dz 1 , m ≥ 0 - 0, n < 0 1 2πi Γ 2 z -m-1 2 [P 0n+1 + z 2 P 0n+2 + ...]A # (z 2 )dz 2 , n ≥ 0, ( 16
)
where

η m := 1, m ≥ 0, 0, m < 0.
For the sake of clarity of presentation, the detailed proof is postponed to the Appendix. Now the recursive formula for the coefficients of the regular part of the Laurent series (2) becomes a corollary of the above general result.

Corollary 1 Suppose that the coefficients H k , k = -s, ..., -1, 0 and P 0k , k = 0, 1, ... are given. Then, the coefficients of the regular part of the Laurent expansion (2) can be computed by the following recursive formula:

H m+1 = - m+s i=0   s j=0 H -j A i+j+1   H m-i - m i=1 P 0m+1-i H i (17)
-(P 0m+1 H 0 + ... + P 0m+1+s H -s ) -(H -s P 0m+1+s + ... + H 0 P 0m+1 )
for m = 0, 1, ... .

Proof:

Let us take n = 0, m > 0 and then simplify the last two terms in [START_REF] Hassin | Mean passage times and nearly uncoupled Markov chains[END_REF], so that 1 2πi

Γ 1 z -n-1 1 A # (z 1 )[P 0m+1 + z 1 P 0m+2 + ...]dz 1 = 1 2πi Γ 1 1 z 1 [ 1 z s 1 H -s + ... + 1 z 1 H -1 + H 0 + ...][P 0m+1 + z 1 P 0m+2 + ...]dz 1 = 1 2πi Γ 1 1 z 1 1 z m+1 1 1 z s 1 [H -s P 0m+1+s + ... + H 0 P 0m+1 ]dz 1 = H -s P 0m+1+s + ... + H 0 P 0m+1 (18)
The last term can be dealt with in a similar fashion.

1 2πi Γ 2 z -m-1 2 [P 0n+1 + z 2 P 0n+2 + ...]A # (z 2 )dz 2 = 1 2πi Γ 2 z -m-1 2 [P 01 + z 2 P 02 + ...]A # (z 2 )dz 2 = 1 2πi Γ 2 1 z 2 1 Z m 2 [P 01 + z 2 P 02 + ...][ 1 z s 1 H -s + ... + 1 z 1 H -1 + H 0 + ...]dz 2 = 1 2πi Γ 2 1 z 2 [P 01 H m + P 02 H m-1 + ... + P 0m H 1 + P 0m+1 H 0 + ... + P 0m+1+s H -s ]dz 2 = m i=1 P 0m+1-i H i + (H -s P 0m+1+s + ... + H 0 P 0m+1 ) (19)
Substituting ( 18) and ( 19) into [START_REF] Hassin | Mean passage times and nearly uncoupled Markov chains[END_REF] with n = 0 and m > 0, we obtain

∞ k=1 k-1 i=0 H n-i A k H m+i-k+1 = -H m+1 -(H -s P 0m+1+s + ... + H 0 P 0m+1 ) -(P 0m+1 H 0 + ... + P 0m+1+s H -s ) - m i=1 P 0m+1-i H i .
Rearranging terms in the above expression, we get

H m+1 = - m+s i=0   s j=0 H -j A i+j+1   H m-i - m i=1 P 0m+1-i H i -(P 0m+1 H 0 + ... + P 0m+1+s H -s ) -(H -s P 0m+1+s + ... + H 0 P 0m+1 )
which is the recursive formula [START_REF] Haviv | Series expansions for stochastic matrices[END_REF]. ✷ If the perturbed operator A(z) is invertible for 0 < |z| < z max , then the inverse A -1 (z) can be expanded as a Laurent series

A -1 (z) = 1 z -s H -s + ... + 1 z H -1 + H 0 + zH 1 + ... (20)
and the formula [START_REF] Haviv | Series expansions for stochastic matrices[END_REF] becomes

H m+1 = - m+s i=0   s j=0 H -j A i+j+1   H m-i , m = 0, 1, ...
Furthermore, if the perturbed operator is invertible and the perturbation is linear A(z) = A 0 + zA 1 we retrieve the recursive formula from [START_REF] Langenhop | The Laurent expansion for a nearly singular matrix[END_REF] H m+1 = (-H 0 A 1 )H m , m = 0, 1, ... .

Reduction process

We have seen that the regular terms in the Laurent series expansion (2) of A # (z) can be computed recursively by [START_REF] Haviv | Series expansions for stochastic matrices[END_REF]. However, to apply [START_REF] Haviv | Series expansions for stochastic matrices[END_REF], one first needs to compute the terms H -s , . . . , H -1 , H 0 , that is, the terms of the singular part.

The complex analytic approach allows us to treat the cases regular and singular perturbations in a unified framework. In fact, let us first obtain some results on the regular perturbation which will be useful in the reduction process for the singular perturbation.

Regular case. Let us apply analytic function techniques to express the power series for the Drazin generalized inverse of the perturbed operator in the case of regular perturbation, that is, when the dimension of the null space of the matrix does not change if the perturbation parameter deviates from zero. In other words, there is no splitting of the zero eigenvalue of the perturbed matrix at z = 0. The latter implies that the expansion ( 13) is valid in some neighbourhood of z 0 = 0, and any contour Γ 0 := Γ 0 (0) chosen so that it does not enclose eigenvalues other than zero. Namely, the expansion (13) takes the form

A # (z) = A # 0 + ∞ n=1 z n A # n , (21)
where

A # 0 = A # (0), A # n = 1 2πi Γ 0 1 ζ [R (n) (ζ)] dζ and R (n) (ζ) = ν 1 +...+νp=n (-1) p R(ζ)A ν 1 R(ζ)A ν 2 • • • R(ζ)A νp R(ζ).
It turns out that it is possible to express the coefficients A # n in terms of

• the unperturbed Drazin generalized inverse A # (0),

• the eigenprojection P 0 corresponding to the zero eigenvalue of A 0 ,

• the perturbation matrices A n , n = 1, 2, . . ..

The next theorem gives the precise statement.

Theorem 3 Suppose that the operator A 0 is perturbed analytically as in (1) and assume that the zero eigenvalue of A 0 is semi-simple and the perturbation is regular. Then, the matrices A # n , n = 1, 2, . . . in the expansion (21) are given by the following formula

A # n = n p=1 (-1) p ν 1 + ... + ν p = n µ 1 + ... + µ p+1 = p + 1 ν j ≥ 1, µ j ≥ 0 S µ 1 A ν 1 S µ 2 ...A νp S µ p+1 , (22)
where S 0 := -P 0 and S k := (A # 0 ) k , k = 1, 2, ....

Proof: Since Γ 0 encloses only the zero eigenvalue, we have by ( 21)

A # n = 1 2πi Γ 0 1 ζ R (n) (ζ)dζ = ν 1 +...νp=n (-1) p 1 2πi Γ 0 1 ζ R(ζ)A ν 1 R(ζ)...A νp R(ζ)dζ. = ν 1 +...νp=n (-1) p Res ζ=0 { 1 ζ R(ζ)A ν 1 R(ζ)...A νp R(ζ)}
In order to compute the above residue, we replace R(ζ) by its Laurent series (8) in the expression

1 ζ R(ζ)A ν 1 R(ζ)...A νp R(ζ),
and collect the terms with 1/ζ, that is, the terms

σ 1 +...+σ p+1 =0 S σ 1 +1 A ν 1 S σ 2 +1 ...A νp S σp+1 .
Next, we change indicies µ k := σ k + 1, k = 1, ..., p + 1, and rewrite the above sum as:

µ 1 +...+µ p+1 =p+1 S µ 1 A ν 1 S µ 2 ...A νp S µp ,
which yields [START_REF] Langenhop | The Laurent expansion for a nearly singular matrix[END_REF]. ✷

Remark 1 Of course, formula [START_REF] Langenhop | The Laurent expansion for a nearly singular matrix[END_REF] is computationally demanding due to the combinatorial explosion. However, only few terms will be computed by this formula (see the arguments developed below).

Singular case. We now show that by using a reduction process, we can transform the original singular problem into a regular one. Our reduction process can be viewed as complimentary to the existing reduction process based on spectral theory (see, e.g., the book of Kato [START_REF] Kato | Perturbation theory for linear operators[END_REF]) which is applied to the eigenvalue problem. To the best of our knowledge, applying the reduction technique to analytical perturbations of generalized inverses is new.

To develop the reduction technique in the context of the generalized inverses, we need to introduce a new notion of group reduced resolvent. A definition based on spectral representation is as follows: Definition 2 Let A : C n → C n be a linear operator with the spectral representation (3). Then, the group reduced resolvent A #Λ relative to the group of eigenvalues Λ := {λ i } k i=0 is defined as follows:

A #Λ def = p i=k+1    1 λ i P i + m i -1 j=1 (-1) j 1 λ j+1 i D j i    .
We note that the Drazin generalized inverse (see formula [START_REF] Castro-González | Additive perturbations results for the Drazin inverse[END_REF]) is a particular case of the group reduced resolvent. In this case, the group of eigenvalues consists only of the zero eigenvalue.

From our definition, the properties of a group reduced resolvent follow easily. In particular, in the next theorem, we will obtain an alternative analytic expression of the group reduced resolvent that will play a crucial role in perturbation analysis.

Theorem 4 Let A be a linear operator with representation (3). Then, the group reduced resolvent relative to the eigenvalues Λ = {λ i } k i=0 is given by

A #Λ = 1 2πi Γ 1 ζ (A -ζI) -1 dζ, ( 23
)
where Γ is a contour in the complex plane which encloses the set of eigenvalues {λ i } k i=0 , but none of the other eigenvalues {λ i } p i=k+1 .

Proof: It is a well known fact [START_REF] Kato | Perturbation theory for linear operators[END_REF] that the resolvent can be represented by

- p i=0   1 ζ -λ i P i + m i j=1 1 (ζ -λ i ) j+1 D j i   .
Substituting the above expression into the integral of ( 23) yields

1 2πi Γ 1 ζ (A -ζI) -1 dζ = - 1 2πi Γ p i=0   1 ζ(ζ -λ i ) P i + m i j=1 1 ζ(ζ -λ i ) j+1 D j i   dζ.
Using the fact that for every positive integer l

Res ζ=0 1 ζ(ζ -λ) l = 1 (-λ) l and Res ζ=λ 1 ζ(ζ -λ) l = - 1 (-λ) l , we obtain 1 2πi Γ 1 ζ (A -ζI) -1 dζ = k i=0 Res ζ=λ i    - p i=0   1 ζ(ζ -λ i ) P i + m i j=1 1 ζ(ζ -λ i ) j+1 D j i      = p i=k+1    1 λ i P i + m i -1 j=1 (-1) j 1 λ j+1 i D j i    .
According to Definition 2, the latter expression is equal to the group reduced resolvent, so that the proof is complete. ✷

(λ i P i + D i )[I -P ] # = p i=k+1 (λ i P i + D i ) # .
Using formula [START_REF] Castro-González | Additive perturbations results for the Drazin inverse[END_REF], we obtain

(A[I -P ]) # = p i=k+1    1 λ i P i + m i -1 j=1 (-1) j 1 λ j+1 i D j i    .
The latter is equal to the group reduced resolvent A #Λ by Definition 2. ✷ Now equipped with this new notion of group reduced resolvent, we go back to perturbation analysis. The group of the perturbed eigenvalues λ i (z) such that λ i (z) → 0 as z → 0 is called the 0-group. We denote the 0-group of eigenvalues by Ω. The eigenvalues of the 0-group split from zero when the perturbation parameter differs from zero. Since the eigenvalues of the perturbed operator are algebraic functions of the perturbation parameter [START_REF] Kato | Perturbation theory for linear operators[END_REF], each eigenvalue of the 0-group (other than 0) can be written as

λ i (z) = z ν λ iν + o(z ν ), ( 24
)
with λ iν = 0 and ν is a positive rational number. The reduction technique is essentially based on the semi-simplicity assumption of reduced operators [START_REF] Kato | Perturbation theory for linear operators[END_REF], which will be introduced below. Under that assumption, the power ν in [START_REF] Pervozvanski | Theory of Suboptimal Decisions[END_REF] must be an integer. The latter implies that we can partition the 0-group into subsets that we call z l -groups. Namely, we say that the eigenvalue λ i (z) belongs to the z l -group if λ i (z) = z l λ il + o(z l ), with λ il = 0. We denote z l -group by Λ l . Consider now the spectral representation of the perturbed reduced resolvent.

A # (z) = p i=1    1 λ i (z) P i (z) + m i -1 j=1 (-1) j 1 λ j+1 i (z) D j i (z)    ,
where {λ i (z)} k i=1 is the 0-group. From the above formula one can see that in this case, the Laurent expansion for the reduced resolvent A # (z) will possess terms with negative powers of z. Moreover, it turns out that under our assumptions, the z k -group eigenvalues contribute to the terms of the Laurent expansion for A # (z) with negative powers -k, -k + 1, ..., -1 as well as to the regular part of the Laurent expansion.

The basic idea is to first treat the part of the perturbed operator corresponding to the eigenvalues which do not tend to zero as z → 0. Then we treat subsequently the parts of the perturbed operator corresponding to the eigenvalues which belong to the z-group, z 2 -group and so on.

What is helpful is that to treat the part of A(z) corresponding to z k+1 -group we have to perform the same algorithm as for the part of the perturbed operator corresponding to the z k -group. These steps constitute the (finite) reduction process. Now we implement the general idea we have just briefly outlined above. Consider a fixed contour Γ 0 that encloses only the zero eigenvalue of the unperturbed operator A 0 . Note that by continuity of eigenvalues the 0-group of eigenvalues of the perturbed operator A(z) lies inside Γ 0 for z sufficiently small. Therefore, we may define the group reduced resolvent relative to the 0-group of eigenvalues as follows:

A #Ω (z) = 1 2πi Γ 0 1 ζ R(ζ, z)dζ = 1 2πi Γ 0 1 ζ (A(z) -ζI) -1 dζ.
Since A #Ω (z) is an analytic function in some neighbourhood of the origin, it can be expanded as a power series

A #Ω (z) = A #Ω 0 + ∞ i=1 z i A #Ω i . (25)
Note that A #Ω 0 = (A 0 ) # and from Theorem 3 it follows that the other coefficients A #Ω i , i = 1, 2, ... can be calculated by the formula [START_REF] Langenhop | The Laurent expansion for a nearly singular matrix[END_REF]. We would like to emphasize that in general the group reduced resolvent A #Ω (z) is different from the reduced resolvent A # (z). However, we note that A #Ω (z) does coincide with A # (z) in the case of regular perturbations.

Another operator that is used extensively in the reduction process is the group projection

P (z) = 1 2πi Γ 0 R(ζ, z) dζ,
which describes the subspace corresponding to the eigenvalues which split from zero. The group projection is an analytic function in some small neighbourhood of the origin (see e.g. [START_REF] Kato | Perturbation theory for linear operators[END_REF]). Next, as in the classical reduction process [START_REF] Baumgartel | Analytic perturbation theory for matrices and operators[END_REF][START_REF] Kato | Perturbation theory for linear operators[END_REF], we define the restriction B(z) of the operator A(z) to the subspace determined by the group projection P (z), that is,

B(z) := 1 z A(z)P (z) = 1 2πiz Γ 0 ζR(ζ, z) dζ, Assumption S2
The zero eigenvalue of B 0 is semi-simple.

We would like to emphasize that the subsequent reduction steps are totally identical with the first one. At each reduction step, we make the assumption Sk that the analogue of B 0 at step k has a semi-simple 0-eigenvalue. The final result is stated in the next theorem.

Theorem 5 Let Assumptions Sk hold. Then, the reduction process terminates after a finite number of steps, say s, and the perturbed Drazin inverse A # (z) has the following expression:

A # (z) = A #Ω (z) + 1 z B #Ω (z) + 1 z 2 C #Ω (z) + ... + 1 z s Z # (z). ( 29
)
Proof: Consider the first reduction step. Since R(P (z)) and R(I -P (z)) represent a direct decomposition of C n and the subspace R(P (z)) is invariant under the operator A(z), we can write

A # (z) = (A(z)[I -P (z)] + A(z)P (z)) # = (A(z)[I -P (z)]) # + (A(z)P (z)) # = A #Ω (z) + z -1 (z -1 A(z)P (z)) # ,
where Lemma 2 was used to get the first term of the right-hand-side. In view of Assumption S1, the operator B(z) = z -1 A(z)P (z) is analytic in z and hence, one can apply the next reduction step. Similarly, Assumptions Sk, k = 1, 2, . . . guarantee that the reduction process can be carried out. Since the splitting of the zero eigenvalue has to terminate after a finite number of steps [START_REF] Kato | Perturbation theory for linear operators[END_REF], we conclude that the reduction process has to terminate after a finite number of steps as well. Indeed, we successively eliminate the eigenvalues of the z-group, the z 2 -group, etc ... Let λ i (z) = z s λ is + ... be the last eigenvalue which splits from zero. Then the corresponding reduced operator Z(z) is regularly perturbed and the associated reduced resolvent Z # (z) has the power series defined by Theorem 3. This completes the proof. ✷ Summarizing, to obtain the Laurent series for A # (z), there are two cases to distinguish. First, if one needs only few regular terms of A # (z), then it suffices to replace A #Ω (z), B #Ω (z), ... in (29) by their respective power series [START_REF] Schweitzer | Perturbation series expansions for nearly completely-decomposable Markov chains[END_REF] computed during the reduction process. Note that only few terms of the power series A #Ω (z), B #Ω (z), ... are needed. Otherwise, if one wishes to compute a significant number of regular terms, then compute only H -s , . . . , H -1 , H 0 as above (in which case, again, only a few terms of A #Ω (z), B #Ω (z), ... are needed) and then use the recursive formula [START_REF] Haviv | Series expansions for stochastic matrices[END_REF]. Of course, one needs first to compute the power series expansion of the eigenprojection P 0 (z), which can be obtained by a number of methods [START_REF] Avrachenkov | Analytic perturbation theory and its applications[END_REF][START_REF] Baumgartel | Analytic perturbation theory for matrices and operators[END_REF][START_REF] Kato | Perturbation theory for linear operators[END_REF].

Remark 2 If the operator A(z) has an inverse for z = 0, then the above algorithm can be used to calculate its Laurent expansion. Hence, the inversion problem A -1 (z) is a particlular case of the presented above complex analytic approach.

Example 1 As was mentioned in the introduction, the perturbation analysis of the reduced resolvent can be applied directly to the theory of singularly perturbed Markov chains [START_REF] Avrachenkov | The fundamental matrix of singularly perturbed Markov chains[END_REF][START_REF] Delebecque | A reduction process for perturbed Markov chain[END_REF][START_REF] Haviv | Series expansions for stochastic matrices[END_REF][START_REF] Schweitzer | Perturbation series expansions for nearly completely-decomposable Markov chains[END_REF]. Namely, the reduced resolvent of the generator of a Markov chain is just with minus sign the deviation matrix of this chain. The deviation matrix plays a crucial role in the Markov chain theory. For example, it is used to obtain mean first-passage times. Taking into account the above remark, we consider an example of a perturbed Markov chain. Let us consider the following perturbed operator. In this instance, the Laurent expansion for A # (z) has a simple pole. Using the method of [START_REF] Hassin | Mean passage times and nearly uncoupled Markov chains[END_REF] for the determination of the singularity order of the perturbed Markov chains, one can check that

A(z) = A 0 + zA 1 =   0 0 0 0 0.5 -0.5 0 -0.5 0.5   + z   2 -1 -1 -3 1 2 -4 3
A # (z) = 1 z H -1 + H 0 + zH 1 + ...
By applying the reduction process, we compute the singular coefficient H -1 and the first regular coefficient H 0 . Since the zero eigenvalues of the reduced operators are always semi-simple in the case of perturbed Markov chains [START_REF] Delebecque | A reduction process for perturbed Markov chain[END_REF], we conclude from Theorem 5 that

H -1 = B # 0 , and
H 0 = A # 0 + B #Ω 1 .
To compute B # 0 and B #Ω 1 , we need to calculate the first two terms of the expansion for the reduced operator B(z). In particular,

B 0 = P (0)A 1 P (0) =   2 -1 -1 -3.5 1.75 1.75 -3.5 1.75 1.75   B 1 = -(A # 0 A 1 P (0)A 1 P (0) + P (0)A 1 A # 0 A 1 P (0) + P (0)A 1 P (0)A 1 A # 0 ) = 1 8   0 4 -4 - 24
5 19 20 -17 -3   calculated with the help of [START_REF] Schweitzer | The Laurent expansion of pencils that are singular at the origin[END_REF]. Next, we calculate the eigenprojection corresponding to the zero eigenvalue of the operator B 0 , that is,

Q(0) = 1 22   14
B #Ω 1 = Q(0)B 1 (B # 0 ) 2 -B # 0 B 1 B # 0 + (B # 0 ) 2 B 1 Q(0) = 1 2662   - 16
52 -36 -236 41 195 248 -201 -47   . Note that the above H -1 coincides with the one given in [START_REF] Schweitzer | Perturbation series expansions for nearly completely-decomposable Markov chains[END_REF].

Thus, we finally obtain

H -1 = B # 0 = 1 5.5 2   1 -1.75 -1.75   2 -1 -1 = 1 121   8 -4 -4 -14 7
If we have in hand the expansion for the ergodic projection, we can use the recursive formula [START_REF] Haviv | Series expansions for stochastic matrices[END_REF] to compute the regular coefficients. Let us compute by the recursive formula the coefficient H 1 for our example. First, applying the reduction process for the eigenproblem [START_REF] Avrachenkov | Analytic perturbation theory and its applications[END_REF][START_REF] Baumgartel | Analytic perturbation theory for matrices and operators[END_REF][START_REF] Kato | Perturbation theory for linear operators[END_REF], one can compute the coefficients for the expansion of the ergodic projection associated with z and z 2 . 6 Perturbation of the Moore-Penrose generalized inverse and the group inverse Note that the Laurent series for the perturbed group inverse does not always exist. Indeed, the existence of the group inverse of the unperturbed operator does not imply the existence of the group inverse of the perturbed operator.

P 01 = 1 121   2 -

Example 2 Consider

A(z) =   0 z 0 0 0 0 0 0 1   .
The space C n can be decomposed in a direct sum of the null space and range of A(0) but no such decomposition exists if z = 0. Thus, the unperturbed operator A(0) has a group inverse and the perturbed operator does not.

The following is a general sufficient condition for the existence of the Laurent series for the perturbed group inverse.

Theorem 6 Let the group inverse A g (z) of the analytically perturbed matrix A(z) exist in some non-empty (possibly punctured) neighbourhood of z = 0. Then the group inverse A g (z) can be expanded as a Laurent series around z = 0 with a non-zero radius of convergence.

Proof: Whenever the group inverse exists, it coincides with the Drazin generalized inverse. Therefore, the existence of a Laurent series follows from Theorem 1. ✷ As one can see from the following example, even though the Moore-Penrose generalized inverse always exists, it might be not analytic function of the perturbation parameter. A # (z 1 ) -A # (z 2 ) z 2 -z 1 -A # (z 1 )P 0 (z 2 ) z 2 -z 1 + P 0 (z 1 )A # (z 2 ) z 2 -z 1 dz 2 dz 1 .

Thus, we obtain

∞ k=1 k-1 i=0
H n-i A k H m+i-k+1 = I 1 -I 2 + I 3 , with:

I 1 := 1 2πi 2 Γ 1 Γ 2 z -n-1 1 z -m-1 2 A # (z 1 ) -A # (z 2 ) z 2 -z 1 dz 2 dz 1 , I 2 := 1 2πi 2 Γ 1 Γ 2 z -n-1 1 z -m-1 2
A # (z 1 )P 0 (z 2 ) z 2 -z 1 dz 2 dz 1 ,

I 3 := 1 2πi 2 Γ 1 Γ 2 z -n-1 1 z -m-1 2 P 0 (z 1 )A # (z 2 ) z 2 -z 1 dz 2 dz 1 .
Let us separately calculate the integrals I 1 , I 2 and I 3 . The integral I 1 can be written as

I 1 = 1 2πi 2 Γ 1 Γ 2 z -n-1 1 z -m-1 2 z 2 -z 1 A # (z 1 )dz 2 dz 1 - 1 2πi 2 Γ 1 Γ 2 z -n-1 1 z -m-1 2 z 2 -z 1 A # (z 2 )dz 2 dz 1 = 1 2πi 2 Γ 1   Γ 2 z -m-1 2 z 2 -z 1 dz 2   z -n-1 1 A # (z 1 )dz 1 - 1 2πi 2 Γ 2   Γ 1 z -n-1 1 z 2 -z 1 dz 1   z -m-1 2 A # (z 2 )dz 2
In the last equality we used Fubini Theorem to change the order of integration. Using the auxiliary integrals (36) and (37), we obtain

I 1 = 1 2πi Γ 1 (-η m z -m-1 1 )z -n-1 1 A # (z 1 )dz 1 - 1 2πi Γ 2 (-(1 -η n )z -n-1 2 )z -m-1 2 A # (z 2 )dz 2 = - η n + η m -1 2πi Γ 1 z -n-m-2 1 A # (z 1 )dz 1 = -(η n + η m -1)H n+m+1 ,
where the second integral can be taken over Γ 1 by the principle of deformation of contours. We calculate the second integral I 2 as follows: where, in the above expressions, the auxiliary integral (38) and the property A # (z)P 0 (z) = 0 has been used. Now, we calculate the last integral I 3 with the help of the auxiliary integral (39).

I 2 = 1 2πi Γ 1 z -n-1 1 A # (z 1 ) 1 2πi Γ 2 z -m-1 2 P 0 (z 2 ) z 2 -z 1 dz 2 dz 1 =      1 2πi Γ 1 0z -n-1 1 dz 1 , m < 0, -1 2πi Γ 1 z -n-
I 3 = 1 2πi Γ 2 1 2πi Γ 1 z -n-1 1 P 0 (z 1 ) z 2 -z 1 dz 1 z -m-1 2 A # (z 2 )dz 2 =      -1 2πi Γ 2
z -n-m-2 2 P 0 (z 2 )A # (z 2 )dz 2 , n < 0, Finally, summing up the three integrals I 1 , I 2 and I 3 , we obtain the relation [START_REF] Hassin | Mean passage times and nearly uncoupled Markov chains[END_REF]. ✷

-1 2πi Γ 2 z -m-

Lemma 2

 2 Let P = k i=0 P i be the projection corresponding to the group of eigenvalues Λ = {λ i } k i=0 , then A #Λ = (A[I -P ]) # . Proof: Since P i P j = δ ij P i and D i P j = δ ij D i , we have P i [I -P ] = 0, D i [I -P ] = 0 for i = 0, ..., k and P i [I -P ] = P i , D i [I -P ] = D i for i = k + 1, ..., p. Then, (3) and above yields (A[I -P ]) # = p i=0

1 

 1 Note that -A(z) is the generator of a Markov chain. The zero eigenprojection and the reduced resolvent of the unperturbed matrix A 0 is given by

  Now using formula[START_REF] Langenhop | The Laurent expansion for a nearly singular matrix[END_REF] from Theorem 3, we obtain

1 ,

 1 which is not analytic since it depends on z.

1 1A 1 [ 1 z 1 A 1 A 1 ) -z m+1 1 P 0m+1 -z m+2 1 P 1 A # (z 1 )[z m+1 1 P 0m+1 + z m+2 1 P 1 A

 1111111111111 # (z 1 )z -m-1 P 00 + z 1 P 01 + ...+ z m 1 P 0m ]dz 1 , m ≥ 0 -n-m-2 # (z 1 )[P 00 + z 1 P 01 + ... + z m 1 P 0m ]dz 1 , m ≥ # (z 1 )[P 0 (z 0m+2 -...]dz 1 , m ≥ 0m+2 + ...]dz 1 , m ≥ # (z 1 )[P 0m+1 + z 1 P 0m+2 + ...]dz 1 , m ≥ 0

  = -(H 0 A 1 )H 0 -(P 01 H 0 + P 02 H -1 ) -(H -1 P 02 + H 0 P 01 )

		12 10 2 -12 10 2 -12 10	  P 02 =	1 1331	 	32 -192 160 32 -192 160 32 -192 160	 
	Then, according to formula (17), we have
	H 1 =	1 14641	 	-368 1744 -10816 1856 -1488 9072 -2128 12416 -10288	  .

1 2

 2 [P 0n+1 + z 2 P 0n+2 + z 2 2 P 0n+3 + ...]A # (z 2 )dz 2 , n ≥ 0, [P 0n+1 + z 2 P 0n+2 + z 2 2 P 0n+3 + ...]A # (z 2 )dz 2 , n ≥ 0.

		0,		n < 0,
	=	-1 2πi	Γ 2	z -m-1 2

where Γ 0 is some fixed contour enclosing only the zero eigenvalue of the unperturbed operator A 0 . For the operator B(z) to be analytic at zero, we need the following assumption.

Assumption S1

The zero eigenvalue of the operator A 0 is semi-simple, that is the nilpotent operator D 0 corresponding to λ 0 = 0 is equal to zero.

Note that this assumption is not too restrictive. For example, in the case of a self-adjoint perturbation operator, the zero eigenvalue of A 0 is semi-simple. This is in particular the case when one studies the Moore-Penrose generalized inverse of an analytically perturbed matrix since it reduces to study a symmetric perturbation of the Drazin inverse (see Section 6). Whenever Assumption S1 is satisfied, the operator B(z) can be expressed as a power series [START_REF] Kato | Perturbation theory for linear operators[END_REF]

with B 0 = P (0)A 1 P (0), and

where S 0 := -P (0) and

Since the operator B(z) is analytic in some neighbourhood of the origin, we can again construct the expansion for its group reduced resolvent

The coefficients B #Ω i , i = 1, 2, ... are calculated by the formula given in Theorem 3. This is the first reduction step. To continue, we must distinguish between two cases:

(i) If the splitting of the zero eigenvalue terminates, and consequently B(z) is a regular perturbation of B 0 , then B #Ω (z) = B # (z) and the Drazin inverse of the perturbed operator A(z) is given by

By substituting the series expansions [START_REF] Schweitzer | Perturbation series expansions for nearly completely-decomposable Markov chains[END_REF] and [START_REF] Stewart | On the perturbation of pseudo-inverses, projections and linear least squares problems[END_REF] for A #Ω (z) and B # (z) into [START_REF] Vainberg | Theory of branching of solutions of non-linear equations[END_REF], we obtain the Laurent series expansion for A # (z), which has a simple pole at zero. (ii) If the zero eigenvalue splits further, the expression

represents only the group reduced resolvent relative to the eigenvalues constituting the 0-group but not the z-group, and we have to continue the reduction process. In fact, we now consider B(z) as a singular perturbation of B 0 , and repeat the procedure with B(z). The 0-group of eigenvalues of B 0 contains all the z k -groups of A(0) (with k ≥ 2), but not the z-group. Specifically, we construct the next-step reduced operator

where Q(z) is the eigenprojection corresponding to the 0-group of the eigenvalues of B(z).

Again, to ensure that C(z) is an analytic function of z, we assume However, if we restrict A k , k = 0, 1, ... to the matrices with real entries and if z belongs to some interval of the real line, we can state the following existence result.

Theorem 7 Let A † (ε) be the Moore-Penrose generalized inverse of the analytically perturbed matrix

where A k ∈ IR n×m , ε ∈ IR and the series converges for 0 < |ε| < ε max . Then, A † (ε) possesses a Laurent series expansion

in some non-empty punctured vicinity around ε = 0.

Proof: Applying the formula

for the perturbed operator A(ε) yields

Note that the group inverse of a symmetric matrix always exists. Hence, by Theorem 6, (A T (ε)A(ε)) g has a Laurent series expansion and so the Moore-Penrose generalized inverse A † (ε) does. ✷ We would like to emphasize that according to (31), computing the Laurent series of the perturbed Moore-Penrose generalized inverse A † (ε) reduces to computing the Laurent series of a group inverse. Moreover, A T (ε)A(ε) is a symmetric perturbation, that is each term of its power series has a symmetric matrix coefficient. This guarantees that the reduction process defined in the previous section and restricted to the real line is indeed applicable in this case.

7 Asymptotics for P 0 (z)

We know that the eigenprojection P 0 (z) of the perturbed operator corresponding to the identically zero eigenvalue is analytic in some (punctured) neighbourhood of z = 0 (see [START_REF] Avrachenkov | Analytic perturbation theory and its applications[END_REF] and [START_REF] Kato | Perturbation theory for linear operators[END_REF]), that is,

for z sufficiently small but different from zero. For regular perturbations, P 00 is just P 0 (0), and the group projection coincides with the eigenprojection. This is not the case for singular perturbations.

Therefore, an interesting question is how P 00 in (32) relates to the original matrix P 0 (0) in the general case and how the power series (32) can be computed. The answers to these questions are provided below.

Proposition 1 The coefficients of the power series (32) for the perturbed eigenprojection are given by

Proof:

The above formula is obtained by equating the terms with the same power of z in the identity [START_REF] Bart | Meromorphic operator valued functions[END_REF] for the perturbed operators. ✷ Corollary 2 As z → 0, the limit eigenprojection matrix P 00 satisfies

where P 0 (0) is the 0-eigenprojection of the unperturbed matrix A 0 .

Proof: When substituting the Laurent series expansion ( 14) into ( 6)-( 7), and equating the terms of same power, one obtains in particular,

In addition, from A(z)P 0 (z) = 0, we immediately obtain:

A 0 P 00 = 0 so that P 00 = P 0 (0)V,

for some matrix V . Moreover, as P 0 (0) 2 = P 0 (0), we also have P 0 (0)P 00 = P 0 (0) 2 V = P 0 (0)V = P 00 . Therefore, pre-multiplying both sides of (34) by P 0 (0), and using P 0 (0)A 0 = 0, one obtains (33), the desired result. ✷ Hence, (33) relates in a simple manner, the limit matrix P 00 to the original 0-group P 0 (0), in terms of the perturbation matrices A k , k = 1, ..., s, the original matrix P 0 (0), and the coefficients H -k , k = 1, ..., s of the singular part of A(z) # . This shows how the perturbed 0-eigenvectors compare to the unperturbed ones, for small z. Observe that in the case of a linear (or first-order) perturbation, only the singular term H -1 is involved. Finally, the regular case is obtained as a particular case, since then, H -k , k = 1, . . . , s vanish so that P 00 = P 0 (0).

Appendix: The proof of Theorem 2

To prove Theorem 2 we use the Cauchy contour integration and the residue technique. First we present some auxiliary results.

Lemma 3 Let Γ 1 and Γ 2 be two closed positively-oriented contours in the complex plane around zero and let z 1 ∈ Γ 1 , z 2 ∈ Γ 2 . Furthermore, assume that the contour Γ 2 lies inside the contour Γ 1 . Then the following formulae hold:

with

Proof: The proof of formulae (36),(37) is given in the book of Kato [START_REF] Kato | Perturbation theory for linear operators[END_REF] and in the book of Korolyuk and Turbin [START_REF] Korolyuk | Mathematical foundations of the state lumping of large systems[END_REF].

Let us establish the auxiliary integral (38). If m < 0, then the function

is analytic inside the area enclosed by the contour Γ 2 and hence the auxiliary integral (38) is equal to zero by the Cauchy Integral Theorem. To deal with the case m ≥ 0, we first expand the function

Then, according to the Residue Theorem, we have

Thus, we have calculated the integral (38). The same method is applied to calculate the auxiliary integral (39). ✷ The proof of Theorem 2: Each coefficient of the Laurent series (2) can be represented by the contour integral formula

where Γ is a closed positively-oriented contour in the complex plane, which encloses zero but no other eigenvalues of A 0 . Let us substitute (40) into the following expression

As in Lemma 3, we assume without loss of generality that the contour Γ 2 lies inside the contour Γ 1 . Then, we can rewrite the above expressions as double integrals

Using the resolvent-like identity [START_REF] Gohberg | On the local theory of regular analytic matrix functions[END_REF], we obtain


Emmanuelle Anceaume
email: emmanuelle.anceaume@irisa.fr

Yann Busnel
email: yann.busnel@univ-nantes.fr

Sketch ⋆-metric: Comparing Data Streams via Sketching

In this paper, we consider the problem of estimating the distance between any two large data streams in small-space constraint. This problem is of utmost importance in data intensive monitoring applications where input streams are generated rapidly. These streams need to be processed on the fly and accurately to quickly determine any deviance from nominal behavior. We present a new metric, the Sketch ⋆-metric, which allows to define a distance between updatable summaries (or sketches) of large data streams. An important feature of the Sketch ⋆-metric is that, given a measure on the entire initial data streams, the Sketch ⋆-metric preserves the axioms of the latter measure on the sketch. Extensive experiments conducted on both synthetic traces and real data sets allow us to validate the robustness and accuracy of the Sketch ⋆-metric.

I. INTRODUCTION

The main objective of this paper is to propose a novel metric that reflects the relationships between any two discrete probability distributions in the context of massive data streams. Specifically, this metric, designated as Sketch ⋆metric in the following, allows us to efficiently estimate a broad class of distances measures between any two large data streams by computing these distances only using compact synopses or sketches of the streams. The Sketch ⋆metric is distribution-free and makes no assumption about the underlying data volume. It is thus capable of comparing any two data streams, identifying their correlation if any, and more generally, it allows us to acquire a deep understanding of the structure of the input streams. Formalization of this metric is the first contribution of this paper.

The interest of estimating distances between any two data streams is important in data intensive applications. Many different domains are concerned by such analyses including machine learning, data mining, databases, information retrieval, and network monitoring. In all these applications, it is necessary to quickly and precisely process a huge amount of data. For instance, in IP network management, the analysis of input streams will allow us to rapidly detect the presence of anomalies or intrusions when changes in the communication patterns occur [START_REF] Lakhina | Mining anomalies using traffic feature distributions[END_REF]. Actually, the problem of detecting changes or outliers in a data stream is similar to the problem of identifying patterns that do not conform to the expected behavior, which has been an active area of research for many decades. To accurately analyze streams of data, a panel of information-theoretic measures and distances have been proposed to answer the specificities of the analyses. Among them, the most commonly used are the Kullback-Leibler (KL) divergence [START_REF] Kullback | On information and sufficiency[END_REF], or more generically, the fdivergences, introduced by Csiszar, Morimoto and Ali & Silvey [START_REF] Ali | General Class of Coefficients of Divergence of One Distribution from Another[END_REF], [START_REF] Csiszár | Information Measures: A Critical Survey[END_REF], [START_REF] Morimoto | Markov processes and the h-theorem[END_REF], the Jensen-Shannon divergence and the Battacharyya distance [START_REF] Bhattacharyya | On a measure of divergence between two statistical populations defined by their probability distributions[END_REF]. Unfortunately, computing information theoretic measures of distances in the data stream model is challenging essentially because one needs to process a huge amount of data sequentially, on the fly, and by using very little storage with respect to the size of the stream. In addition the analysis must be robust over time to detect any sudden change in the observed streams (which might be the manifestation of routers deny of service attack or worm propagation). We tackle this issue by presenting an approximation algorithm that constructs a sketch of the stream from which the Sketch ⋆-metric is computed. This algorithm is a one-pass algorithm. It uses very basic computations, little storage space (i.e., O(t(log n+k log m))

where k and t are precision parameters, and m and n are respectively the size of the input stream and the number of items in the stream), and does not need any information on the structure of the input stream. This constitutes the second contribution of the paper.

Finally, the robustness of our approach is validated with a detailed experimentation study based on both synthetic traces that range from stable streams to highly skewed ones, and real data sets.

The paper is organized as follows. First, Section II reviews the related work on classical generalized metrics. Section III describes the data stream model. Section IV presents the necessary background that makes the paper self-contained. Section V formalizes the Sketch ⋆-metric. Section VI presents the algorithm that approximates the Sketch ⋆-metric in one pass and Section VII presents extensive experiments of our algorithm. Finally, we conclude in Section VIII.

II. RELATED WORK

Work on data stream analysis mainly focuses on efficient methods (data-structures and algorithms) to answer different kind of queries over massive data streams. Mostly, these methods consist in deriving statistic estimators over the data stream, in creating summary representations of streams (to build histograms, wavelets, and quantiles), and in comparing data streams. Regarding the construction of estimators, a seminal work is due to Alon et al. [START_REF] Alon | The space complexity of approximating the frequency moments[END_REF]. The authors have proposed estimators of the frequency moments F k of a stream, which are important statistical tools that allow to quantify specificities of a data stream. Subsequently, a lot of attention has been paid to the strongly related notion of the entropy of a stream, and all notions based on entropy [START_REF] Cover | Elements of information theory[END_REF]. These notions are essentially related to the quantification of the amount of randomness of a stream (e.g, [START_REF] Anceaume | An information divergence estimation over data streams[END_REF], [START_REF] Chakrabarti | Estimating entropy and entropy norm on data streams[END_REF], [START_REF] Guha | Streaming and sublinear approximation of entropy and information distances[END_REF], [START_REF] Lall | Data streaming algorithms for estimating entropy of network traffic[END_REF]). The construction of synopses or sketches of the data stream have been proposed for different applications (e.g, [START_REF] Charikar | Finding frequent items in data streams[END_REF], [START_REF] Cormode | Sketching probabilistic data streams[END_REF]). Actually in [START_REF] Guha | Sketching information divergences[END_REF], the authors propose a characterization of the information divergences that are not sketchable. They have proven that any distance that has not "norm-like" properties is not sketchable.

Distance and divergence measures are key measures in statistical inference and data processing problems [START_REF] Basseville | On entropies, divergences, and mean values[END_REF]. There exists two largely used broad classes of measures, namely the f -divergences and the Bregman divergences that are very important to quantify the amount of information that separates two distributions. Our goal in this paper is to formalize a metric that allows to efficiently and accurately estimate a broad class of distances measures between any two large data streams by computing these distances only on compact synopses or sketches of streams.

III. DATA STREAM MODEL

We consider a system in which a node P receives a very large data stream σ = a 1 , a 2 , . . . , a m of data items that arrive sequentially. In the following, we describe a single instance of P , but clearly multiple instances of P may coexist in a system (e.g., in case P represents a router, or a base station in a sensor network). Each data item a i of the stream σ is drawn from the universe Ω = {1, 2, . . . , n} where n should be very large. Data items can be repeated multiple times in the stream. In the following, we suppose that the length m of the stream is not known. Items in the stream arrive regularly and quickly, and due to memory constraints, need to be processed sequentially and in an online manner. Therefore, node P can locally store only a small fraction of the items and perform simple operations on them. The algorithms we consider in this work are characterized by the fact that they can approximate some function on σ with a very limited amount of memory. We refer the reader to [START_REF] Muthukrishnan | Data Streams: Algorithms and Applications[END_REF] for a detailed description of data streaming models and algorithms.

IV. INFORMATION DIVERGENCE OF DATA STREAMS

We first present notations and background that make this paper self-contained.

A. Preliminaries

• A natural approach to study a data stream σ is to model it as an empirical data distribution over the universe

B. Metrics and divergences 1) Metric definitions:

The classical definition of a metric is based on a set of four axioms.

Definition 1 (Metric) Given a set X, a metric is a function d : X × X → R such that, for any x, y, z ∈ X, we have: Non-negativity: d(x, y) ≥ 0

(1)

Identity of indiscernibles: d(x, y) = 0 ⇔ x = y (2) Symmetry: d(x, y) = d(y, x) (3)
Triangle inequality:

d(x, y) ≤ d(x, z) + d(z, y) (4)
In the context of information divergence, usual distance functions are not precisely metric. Indeed, most of divergence functions do not verify the 4 axioms, but only a subset of them. For instance, a pseudometric is a function that verifies the axioms of a metric with the exception of the identity of indiscernible, while a premetric is a pseudometric that relax both the symmetry and the triangle inequality axioms.

Two classes of generalized metrics, usually denoted as divergences, that allow to measure the separation of distributions have been proposed, namely the class of f -divergences and the class of Bregman divergences. Note that in the following by abusing the notation, we denote "|Ω|-point distribution" by "Ω-point distribution".

2) f -divergence: The class of f -divergences provides a set of relations that is used to measure the "distance" between two distributions p and q. Mostly used in the context of statistics and probability theory, a f -divergence D f is a premetric that guarantees monotonicity and convexity.

Definition 2 (f -divergence) Let p and q be two Ω-point distributions. Given a convex function f : (0, ∞) → R such that f (1) = 0, the f -divergence of q from p is D f (p||q) = i∈Ω q i f p i q i ,
where by convention 0f ( 0 0 ) = 0, af ( 0 a ) = a lim u→0 f (u), and 0f ( a 0 ) = a lim u→∞ f (u)/u if these limits exist.

Property 3 (Monotonicity) Given κ an arbitrary transition probability that respectively transforms two Ω-point distributions p and q into p κ and q κ , we have:

D f (p||q) ≥ D f (p κ ||q κ ).
Property 4 (Convexity) Let p 1 , p 2 , q 1 and q 2 be four Ωpoint distributions. Given any λ ∈ [0, 1], we have:

D f (λp 1 + (1 -λ)p 2 ||λq 1 + (1 -λ)q 2 ) ≤ λD f (p 1 ||q 1 ) + (1 -λ)D f (p 2 ||q 2 ).
3) Bregman divergence: Initially proposed in [START_REF] Bregman | The relaxation method of finding the common point of convex sets and its application to the solution of problems in convex programming[END_REF], the Bregman divergences are a generalization of the notion of distance between points. This class of generalized metrics always satisfies the non-negativity and identity of indecernibles. However they do not always satisfy the triangle inequality and their symmetry depends on the choice of the differentiable convex function F . Specifically, Definition 5 (Bregman divergence (BD)) Given a continuously-differentiable and strictly convex function F defined on a closed convex set C, the Bregman divergence of p from q is

B F (p||q) = F (p) -F (q) -∇F (q), (p -q) .
where the operator •, • denotes the inner product, and ∇F (q) is the gradient of F at q.

In the context of data stream, it is possible to reformulate this definition as follows. Specifically, Definition 6 (Decomposable BD) Let p and q be any two Ω-point distributions. Given a strictly convex function F : (0, 1] → R, the Bregman divergence of q from p is defined as

B F (p||q) = i∈Ω (F (p i ) -F (q i ) -(p i -q i )F ′ (q i )) .
The Bregman divergence verifies non-negativity and convexity properties in its first argument, but not necessarily in the second argument. Another interesting property is given by thinking of the Bregman divergence as an operator of the function F . Property 7 (Linearity) Let F 1 and F 2 be any two strictly convex and differentiable functions. Given any λ ∈ [0, 1], we have that

B F1+λF2 (p||q) = B F1 (p||q) + λB F2 (p||q).

4) Classical metrics:

Based on these definitions, we present several commonly used metrics in Ω-point distribution context. These specific metrics are used in the evaluation part presented in Section VII.

Kullback-Leibler divergence:

The Kullback-Leibler (KL) divergence [START_REF] Kullback | On information and sufficiency[END_REF], also called the relative entropy, is a robust metric for measuring the statistical difference between two data streams. The KL divergence owns the special feature that it is both a f -divergence and a Bregman one (with f (t) = F (t) = t log t).

Given p and q two Ω-point distributions, the Kullback-Leibler divergence is defined as

D KL (p||q) = i∈Ω p i log p i q i . ( 5
)
Jensen-Shannon divergence: The Jensen-Shannon divergence (JS) is a symmetrized version of the Kullback-Leibler divergence. Also known as information radius (IRad) or total divergence to the average, it is defined as

D JS (p||q) = 1 2 D KL (p||ℓ) + 1 2 D KL (q||ℓ), (6)
where ℓ = 1 2 (p + q). Note that the square root of this divergence is a metric.

Bhattacharyya distance: The Bhattacharyya distance is derived from his proposed measure of similarity between two multinomial distributions, also known as the Bhattacharya coefficient (BC) [START_REF] Bhattacharyya | On a measure of divergence between two statistical populations defined by their probability distributions[END_REF]. It is a semimetric as it does not verify the triangle inequality. It is defined as

D B (p||q) = -log(BC(p, q)) where BC(p, q) = i∈Ω √ p i q i .

V. SKETCH ⋆-METRIC

We now present a method to sketch two input data streams σ 1 and σ 2 , and to compute any generalized metric φ between these sketches such that this computation preserves all the properties of φ computed on σ 1 and σ 2 .

Definition 8 (Sketch ⋆-metric) Let p and q be any two Ωpoint distributions. Given a precision parameter k, and any generalized metric φ on the set of all Ω-point distributions, there exists a Sketch ⋆-metric φ k defined as follows

φ k (p||q) = max ρ∈P k (Ω) φ( p ρ || q ρ ),
with ∀a ∈ ρ, p ρ (a) = i∈a p i and where P k (Ω) is the set of all partitions of Ω into exactly k nonempty and mutually exclusive cells.

Remark 9 Note that for k > n, it does not exist a partition of Ω into k nonempty parts. By convention, we consider that φ k (p||q) = φ(p||q) in this specific context.

In this section, we focus on the preservation of axioms and properties of a generalized metric φ by the corresponding Sketch ⋆-metric φ k .

A. Axioms preserving

Theorem 10 Given any generalized metric φ then, for any k ∈ N, the corresponding Sketch ⋆-metric φ k preserves all the axioms of φ.

Proof: The proof is directly derived from Lemmata 16-19 in the companion paper [START_REF] Anceaume | Sketch ⋆-metric: Comparing Data Streams via Sketching[END_REF]. The first three ones say that using sets operations and sum we get that (i) from non-negative numbers it is impossible to generate negative numbers, (ii) 0 always remains 0, and (iii) it is impossible to generate asymmetry. Finally, the triangle inequality is preserved as there exists ρ ∈ P k (Ω) a k-cell partition such that φ( p ρ || q ρ ) = max ρ∈P k (Ω) φ( p ρ || q ρ ). Due to space constraints, the interested readers are invited to look through [START_REF] Anceaume | Sketch ⋆-metric: Comparing Data Streams via Sketching[END_REF] for detailed proofs of these four technical lemmata.

B. Properties preserving

Theorem 11 Given a f -divergence φ then, for any k ∈ N, the corresponding Sketch ⋆-metric φ k is also a f -divergence.

Proof: From Theorem 10, φ k preserves the axioms of the generalized metric. Thus, φ k and φ are in the same equivalence class. Moreover, from Lemma 13, φ k verifies the monotonicity property. Thus, as the f -divergence is the only class of decomposable information monotonic divergences (cf. [START_REF] Csiszár | Information Measures: A Critical Survey[END_REF]), φ k is also a f -divergence. Theorem 12 Given a Bregman divergence φ then, for any k ∈ N, the corresponding Sketch ⋆-metric φ k is also a Bregman divergence.

Proof: From Theorem 10, φ k preserves the axioms of the generalized metric. Thus, φ k and φ are in the same equivalence class. Moreover, the Bregman divergence is characterized by the property of transitivity (cf. [START_REF] Csiszár | Why least squares and maximum entropy? an axiomatic approach to inference for linear inverse problems[END_REF]) defined as follows. Given p, q and r three Ω-point distributions such that q = Π(L|r) and p ∈ L, with Π is a selection rule according to the definition of Csiszár in [START_REF] Csiszár | Why least squares and maximum entropy? an axiomatic approach to inference for linear inverse problems[END_REF] and L is a subset of the Ω-point distributions, we have the Generalized Pythagorean Theorem:

φ(p||q) + φ(q||r) = φ(p||r).
Moreover the authors in [START_REF] Amari | Information geometry of divergence functions[END_REF] show that the set S n of all discrete probability distributions over n elements ({x 1 , . . . , x n }) is a Riemannian manifold, and it owns another different dually flat affine structure. They also show that these dual structures give rise to the generalized Pythagorean theorem. This is verified for the coordinates in S n and for the dual coordinates [START_REF] Amari | Information geometry of divergence functions[END_REF]. Combining these results with the projection theorem [START_REF] Csiszár | Why least squares and maximum entropy? an axiomatic approach to inference for linear inverse problems[END_REF], [START_REF] Amari | Information geometry of divergence functions[END_REF], we obtain that

φ k (p||r) = max ρ∈P k (n) φ( p ρ || r ρ ) = max ρ∈P k (n) (φ( p ρ || q ρ ) + φ( q ρ || r ρ )) = max ρ∈P k (n) φ( p ρ || q ρ ) + max ρ∈P k (n) φ( q ρ || r ρ ) = φ k (p||q) + φ k (q||r)
Finally, by the characterization of Bregman divergence through transitivity [START_REF] Csiszár | Why least squares and maximum entropy? an axiomatic approach to inference for linear inverse problems[END_REF], and reinforced with Lemma 15 statement, φ k is also a Bregman divergence.

In the following, we show that the Sketch ⋆-metric preserves the properties of divergences.

Lemma 13 (Monotonicity) Given any generalized metric φ verifying the Monotonicity property then, for any k ∈ N, the corresponding Sketch ⋆-metric φ k preserves the Monotonicity property.

Proof: Let p and q be any two Ω-point distributions. Given c < n, consider a partition µ ∈ P c (Ω). As φ is monotonic, we have φ(p||q) ≥ φ( p µ || q µ ) [START_REF] Amari | α-Divergence Is Unique, Belonging to Both f -Divergence and Bregman Divergence Classes[END_REF]. We split the proof into two cases:

Case [START_REF] Ali | General Class of Coefficients of Divergence of One Distribution from Another[END_REF]. Suppose that c ≥ k. Computing φ k ( p µ || q µ ) amounts in considering only the k-cell partitions ρ ∈ P k (Ω) that verify ∀b ∈ µ, ∃a ∈ ρ : b ⊆ a.

These partitions form a subset of P k (Ω). The maximal value of φ( p ρ || q ρ ) over this subset cannot be greater than the maximal value over the whole P k (Ω). Thus we have

φ k (p||q) = max ρ∈P k (Ω) φ( p ρ || q ρ ) ≥ φ k ( p µ || q µ ).
Case [START_REF] Alon | The space complexity of approximating the frequency moments[END_REF]. Suppose now that c < k. By definition, we have

φ k ( p µ || q µ ) = φ( p µ || q µ ). Consider ρ ′ ∈ P k (Ω) such that ∀a ∈ ρ ′ , ∃b ∈ µ, a ⊆ b.
It then exists a transition probability that respectively transforms p ρ ′ and q ρ ′ into p µ and q µ . As φ is monotonic, we have

φ k (p||q) = max ρ∈P k (Ω) φ( p ρ || q ρ ) ≥ φ( p ρ ′ || q ρ ′ ) ≥ φ( p µ || q µ ) = φ k ( p µ || q µ ).
Finally for any value of c, φ k guarantees the monotonicity property. This concludes the proof. Lemma 14 (Convexity) Given any generalized metric φ verifying the Convexity property then, for any k ∈ N, the corresponding Sketch ⋆-metric φ k preserves the Convexity property.

Proof: Let p 1 , p 2 , q 1 and q 2 be any four Ω-point distributions. Given any λ ∈ [0, 1], we have:

φ k (λp 1 + (1 -λ)p 2 ||λq 1 + (1 -λ)q 2 ) = max ρ∈P k (Ω) φ λ p 1ρ + (1 -λ) p 2ρ ||λ q 1ρ + (1 -λ) q 2ρ Let ρ ∈ P k (Ω) such that φ λ p 1ρ + (1 -λ) p 2ρ ||λ q 1ρ + (1 -λ) q 2ρ = max ρ∈P k (Ω) φ λ p 1ρ + (1 -λ) p 2ρ ||λ q 1ρ + (1 -λ) q 2ρ .
Input: Two input streams σ 1 and σ 2 ; the distance φ, k and t settings; Output: The distance φ between σ 1 and σ 2 1 Choose t functions h : [n] → [k], each from a 2-universal hash function family;

2 C σ1 [1...t][1...k] ← 0, C σ2 [1...t][1...k] ← 0; 3 for a j ∈ σ 1 do 4 v ← a j ; 5 for i = 1 to t do 6 C σ1 [i][h i (v)] ← C σ1 [i][h i (v)] + 1; 7 for a j ∈ σ 2 do 8 w ← a j ; 9 for i = 1 to t do 10 C σ2 [i][h i (w)] ← C σ2 [i][h i (w)] + 1; 11 On query φk (σ 1 ||σ 2 ) return φ = max 1≤i≤t φ(C σ1 [i][-], C σ2 [i][-]); Figure 1. Sketch ⋆-metric algorithm
As φ verifies the Convexity property, we have:

φ k (λp 1 + (1 -λ)p 2 ||λq 1 + (1 -λ)q 2 ) = φ λ p 1ρ + (1 -λ) p 2ρ ||λ q 1ρ + (1 -λ) q 2ρ ≤ λφ( p 1ρ || q 1ρ ) + (1 -λ)φ( p 2ρ || q 2ρ ) ≤ λ max ρ∈P k (Ω) φ( p 1ρ || q 1ρ ) + (1 -λ) max ρ∈P k (Ω) φ( p 2ρ || q 2ρ ) = λ φ k (p 1 ||q 1 ) + (1 -λ) φ k (p 2 ||q 2 )
that concludes the proof.

Lemma 15 (Linearity) The Sketch ⋆-metric definition preserves the Linearity property.

Proof: For space limitation reasons, proof is presented in the companion paper [START_REF] Anceaume | Sketch ⋆-metric: Comparing Data Streams via Sketching[END_REF].

To summarize, we have shown that the Sketch ⋆-metric preserves all the axioms of a metric as well as the properties of f -divergences and Bregman divergences. We now show how to efficiently implement such a metric.

VI. APPROXIMATION ALGORITHM

In this section, we propose an algorithm that computes the Sketch ⋆-metric in one pass on the stream. By definition of the metric (cf. Definition 8), we need to generate all the possible k-cell partitions. The number of these partitions follows the Stirling numbers of the second kind, which is equal to

S(n, k) = 1 k! k j=0 (-1) k-j k j j n ,
where n is the size of the items universe. Therefore, S(n, k) grows exponentially with n. As the generating function of S(n, k) is equivalent to x n , it is unreasonable in term of space complexity. We show in the following that generating t = ⌈log(1/δ)⌉ random kcell partitions, where δ is the probability of error of our randomized algorithm, is sufficient to guarantee good overall performance of our metric.

Our algorithm is inspired from the Count-Min Sketch algorithm proposed by Cormode and Muthukrishnan [START_REF] Cormode | An improved data stream summary: the count-min sketch and its applications[END_REF]. Specifically, the Count-Min algorithm is an (ε, δ)-approximation algorithm that solves the frequency-estimation problem. For any items in the input stream σ, the algorithm outputs an estimation fv of the frequency of item v such that P{| fv -

f v | > ε(m -f v )} < δ
, where m represent the size of the input stream and ε, δ > 0 are given as parameters of the algorithm. The estimation is computed by maintaining a two-dimensional array C of t × k counters, and by using t 2-universal hash functions h i (1 ≤ i ≤ t), where k = 2/ε and t = ⌈log(1/δ)⌉. Each time an item v is read from the input stream, this causes one counter of each line to be incremented, i.e.,

C[i][h i (v)] is incremented by one for each i ∈ [1..t].
To compute the Sketch ⋆-metric of two streams σ 1 and σ 2 , two sketches σ 1 and σ 2 of these streams are constructed according to the above description. Note that there is no particular assumption on the length of both streams σ 1 and σ 2 . That is their respective length is finite but unknown. By construction of the 2-universal hash functions h i (1 ≤ i ≤ t), the i th line of C σ1 and C σ2 corresponds to the same partition ρ i of the Ω-point empirical distributions of both σ 1 and σ 2 . Thus when a query is issued to compute the given distance φ between these two streams, the maximal value over all the t partitions ρ i of the distance φ between σ 1ρ i and σ 2ρ i is returned, i.e., the distance φ applied to the i th lines of C σ1 and C σ2 for 1 ≤ i ≤ t. Figure 1 presents the pseudo-code of our algorithm.

Lemma 16

Given parameters k and t, Algorithm 1 gives an approximation of the Sketch ⋆-metric, using O (t(log n + k log m)) bits of space.

Proof: The matrices C σi , for any i ∈ {1, 2}, are composed of t × k counters, which uses O (log m). On the other hand, with a suitable choice of hash family, we can store the hash functions above in O(t log n) space.

VII. PERFORMANCE EVALUATION

A. Settings of the experiments

We have implemented our Sketch ⋆-metric and have conducted a series of experiments on different types of streams and for different parameters settings. We have fed our algorithm with both real-world data sets and synthetic traces. Real data give a realistic representation of some existing systems, while the latter ones allow to capture phenomenon which may be difficult to obtain from real-world traces, and thus allow to check the robustness of our metric. We have varied all the significant parameters of our algorithm, that is, the maximal number of distinct data items n in each stream, the number of cells k of each generated partition, and the number of generated partitions t. Synthetic traces of streams have been generated from 7 distributions showing very different shapes, that is the Uniform distribution, the Zipfian or power law one with parameter α = 1, 2, 4, the Poisson distribution with parameter λ = N/2, the Binomial and the Negative Binomial (or Pascal) ones. For each parameters setting, we have conducted and averaged 100 trials of the same experiment, leading to a total of more than 300, 000 experiments for the evaluation of our metric. Real data have been downloaded from the repository of Internet network traffic [START_REF]the Internet Traffic Archive[END_REF]. We have used five large traces among the available ones. Two of them represent two weeks logs of HTTP requests to the Internet service provider ClarkNet WWW server -ClarkNet is a full Internet access provider for the Metro Baltimore-Washington DC area -the other two ones contain two months of HTTP requests to the NASA Kennedy Space Center WWW server, and the last one represents seven months of HTTP requests to the WWW server of the University of Saskatchewan, Canada. Table I presents some statistics of these data traces. Note that all these benchmarks share a Zipfian behavior, with a lower α for the University of Saskatchwan.

B. Main lessons drawn from the experiments

In this section, we evaluate the accuracy of the Sketch ⋆-metric by comparing φ k (p||q) with φ k (p||q), for φ ∈ {Kullback-Leiber, Jensen-Shannon, Bhattacharyya}, and for p and q generated from the 7 distributions and the 5 real data sets. Distances computed from the sketches of the stream are referred to as Sketch in the legend of the graphs, while the ones computed from the full streams are mentioned as Ref.

Due to space constraints, only a subset of the results are presented in the paper. The interested reader is invited to read [START_REF] Anceaume | Sketch ⋆-metric: Comparing Data Streams via Sketching[END_REF] for the complete evaluation.

Figure 2 shows the accuracy of our metric as a function of the different input streams and the different generalized metrics applied on these streams. The first noticeable remark is that Sketch ⋆-metric behaves perfectly well when the two compared streams follow the same distribution, whatever the generalized metric φ used. This can be observed from both synthetic traces (cf. Figure 2(a) with both p and q following the Pascal distribution, and Figure 2(b) with both p and q uniformly distributed), and real data sets (cf. Figures 2(c) and 2(d) with the NASA (July and August) and ClarkNet (August and September) traces). This tendency is further observed when the distributions of input streams are close to each other (e.g., Zipf-α = 2, 4 and Pascal distributions, or Uniform and Zipf-α = 1). This makes the Sketch ⋆-metric a very good candidate as a parametric method for making distribution parameters inference. Another interesting result is shown when the two input streams exhibit a totally different shape. Specifically, let us consider Figures 2(a) and 2(b). Sketching the Uniform distribution leads to k-cell partitions whose value is well distributed, that is, for a given partition φ, all the k cell values have with high probability the same value. Now, when sketching the Pascal distribution, the repartition of the data items in the cells of any given partitions is such that a few number of data items (those with high frequency) populate a very few number of cells. However, the values of these cells is very large compared to the other cells, which are populated by a large number of data items whose frequency is small. Thus, the contribution of data items exhibiting a small frequency and sharing the cells of highly frequent items is biased compared to the contribution of the other items. Thus although the input streams show a totally different shape, the accuracy of φ k is only slightly lowered in these scenarios which makes it a very powerful tool to compare any two different data streams. The same observation holds with real data sets. When the shapes of the input streams are different (which is the case for Saskatchewan with respect to the 4 other input streams), the accuracy of the Sketch ⋆-metric decreases a little bit but in a very small proportion. Notice that the scales on the y-axis differ significantly in Figures 2(a We have also observed in [START_REF] Anceaume | Sketch ⋆-metric: Comparing Data Streams via Sketching[END_REF] the strong impact of the non-symmetry of the Kullback-Leibler divergence on the computation of the distance (computed on full streams or on sketches) with a clear influence when the input streams follow a Pascal and Zipf-α = 1 distributions.

Figure 3 summarizes the good properties of φ k by illustrating how, for any generalized metric φ, and for any variations in the shape of the two input distributions, φ k remains close to φ. Recall that increasing values of the r parameter of the Negative Binomial distribution makes the shape of the distribution flatter, while maintaining the same mean value.

Figure 4 presents the impact of the number of cells per generated partition on the accuracy of the ⋆-metric on both synthetic traces and real data. It clearly shows that by increasing k the number of data items per cell in the generated partition shrinks and thus the absolute error on the computation of the distance decreases. The same feature appears when the number n of distinct data items in the stream increases. Indeed, when n increases (for a given k), the number data items per cell augments and thus the precision of our metric decreases. This gives rise to a shift of the inflection point, as illustrated in Figure 4(b) as data sets have almost twenty to forty times more distinct data items than the synthetic ones. As aforementioned, the input streams exhibit very different shapes which explain the strong impact of k. Note also that k has the same influence on the Sketch ⋆-metric for all the generalized distances φ.

Finally, it is interesting to note that the number t of generated partitions has a slight influence on the accuracy of our metric The reason comes from the use of 2-universal hash functions, which guarantee for each of them and with high probability that data items are uniformly distributed over the cells of any partition. As a consequence, augmenting the number of such hash functions has a weak influence on the accuracy of the metric.

VIII. CONCLUSION AND OPEN ISSUES

In this paper, we have introduced a new metric, the Sketch ⋆-metric, that allows to compute any generalized metric φ on the summaries of two large input streams. We have presented a simple and efficient algorithm to sketch streams and compute this metric, and we have shown that it behaves pretty well whatever the considered input streams. We are convinced of the undisputable interest of such a metric in various domains including machine learning, data mining, databases, information retrieval and network monitoring.

Regarding future works, we plan to consider a distributed setting, where each site would be in charge of analyzing its own streams and then would propagate its results to the other sites of the system for comparison or merging. An immediate application of such a tool would be to detect massive attacks in a decentralized manner (e.g., by identifying specific connection profiles as with worms propagation, and massive port scan attacks or by detecting sudden variations in the volume of received data).

  )-2(b) and in Figures 2(c)-2(d).

Figure 2 .

 2 Figure 2. Comparison between the Sketch ⋆-metric and the φ metric as a function of the input stream q either generated from a distribution or real traces. For synthetic traces, m = 200, 000 and n = 4, 000. Parameters of the count-min sketch data structure are k = 200 and t = 4. All the histograms share the same legend, but for readability reasons, this legend is only indicated on histogram 2(a).

  distance Ref -Kullback-Leibler divergence Sketch -Kullback-Leibler divergence Ref -Jensen-Shannon divergence Sketch -Jensen-Shannon divergence

Figure 3 .

 3 Figure 3. Comparison between the Sketch ⋆-metric and the φ metric as a function of the parameters of the Negative Binomial distribution NB(r, n/(2r + n)), where distribution p follows a Uniform distribution and q follows the Negative Binomial distribution NB(r, n/(2r + n)).

  distance Ref -Kullback-Leibler divergence Sketch -Kullback-Leibler divergence Ref -Jensen-Shannon divergence Sketch -Jensen-Shannon divergence (a) Synthetic traces -Distribution p follows a Uniform distribution and q follows a Negative Binomial NB(3; 0.99) distance Sketch -Bhattacharyya distance Ref -Kullback-Leibler divergence Sketch -Kullback-Leibler divergence Ref -Jensen-Shannon divergence Sketch -Jensen-Shannon divergence (b) Real datasets -The input stream p is the ClarkNet (August) trace and q is the Saskatchewan one

Figure 4 .

 4 Figure 4. Comparison between the Sketch ⋆-metric and the φ metric as a function of the number of cells k per partition (the number of partitions t of the count-min sketch data structure is set to 4). For synthetic traces, m = 200, 000 and n = 4, 000.

  Ω, given by (p 1 , p 2 , . . . , p n ) with p i = x i /m, and x i = |{j : a j = i}| representing the number of times data item i appears in σ. Note that x i , the number of times item i appears in a stream, is commonly called the frequency of i. We have m = i∈Ω x i .

	• 2-universal Hash Functions A collection H of hash functions h : {1, . . . , M } → {0, . . . , M ′ } is said to be 2-universal if for every h ∈ H and for every two M ′ , different items i, j ∈ [M ], P{h(i) = h(j)} ≤ 1 which is exactly the probability of collision obtained if
	the hash function assigned truly random values to any
	i ∈ [M ], where notation [M ] means {1, . . . , M }.

Table I STATISTICS

 I OF THE FIVE REAL DATA TRACES.

	Data trace	# items (m)	# distinct items (n)	max. freq.
	NASA (July)	1,891,715	81,983	17,572
	NASA (August)	1,569,898	75,058	6,530
	ClarkNet (August)	1,654,929	90,516	6,075
	ClarkNet (September)	1,673,794	94,787	7,239
	Saskatchewan	2,408,625	162,523	52,695