Paul-Marie Samson
email: paul-marie.samson@univ-mlv.fr

Keywords: concentration inequalities, transportation inequalities, infimum-convolution operator, empirical processes, bin packing problem

Infimum-convolution description of concentration properties of product probability measures, with applications.

Description de propriétés de concentration des probabilités produit par opérateurs d'infimum-convolution, et applications.

1 Introduction to the main theorems.

The concentration of measure phenomenon on product spaces has been widely investigated by M. Talagrand in [START_REF] Talagrand | Concentration of measure and isoperimetric inequalities in product spaces[END_REF], [START_REF] Talagrand | New concentration inequalities in product spaces[END_REF] and [START_REF] Talagrand | A new look at independence[END_REF]. In these papers, a large variety of powerful dimension-free concentration inequalities are obtained by induction over the dimension of the product space. Consider a probability space (E 1 , E 1 , µ 1 ) and its product E = E n 1 , µ = µ ⊗n 1 . The basic idea of concentration is that if the measure of a subset A ⊂ E is not too small (µ(A) ≥ 1/2) then most of the points x ∈ E are "close" to A. Talagrand expressed it in the following form: for a given "distance" f (x, A) between x and A one has

e Kf (x,A) dµ(x) ≤ 1 µ(A) , ( 1
)
where K is a non-negative constant. If µ(A) is not too small, using Chebyshev inequality, it follows that the set of elements x ∈ E which are "far" from A is of small measure, since for t ≥ 0,

µ ({x ∈ E, f (x, A) ≥ t}) ≤ 1 µ(A) e -Kt .
The main issue here is to define an interesting notion of "closeness". The results of the present paper are connected with the distances associated to the socalled convex hull approach from [START_REF] Talagrand | Concentration of measure and isoperimetric inequalities in product spaces[END_REF] and [START_REF] Talagrand | New concentration inequalities in product spaces[END_REF]. We obtain refinements of Talagrand's results by extending the so-called infimum-convolution description of concentration introduced by Maurey [Mau] (see also [B-G-L]). One of the main motivation for these investigations is to provide new optimal deviation bounds for suprema of sums of random variables (see section 3).

Our approach has some of its origins in the so-called "convex property (τ )" of [Mau] (see also [Sam]). It is a variant of Maurey's "property (τ )". This property was studied by several authors in connection with concentration properties of log-concave measures, as the Gaussian and the exponential measure (see [Mau], [B-G-L], [B-G], [Sch]). The "convex property (τ )" is a dimension-free functional inequality which is valid for every product measure µ := µ 1 ⊗ • • • ⊗ µ n , when each µ i is a probability measure on a normed vector space (F i , • ) supported by a set of diameter less than one. Maurey's result states: for every convex measurable function f on the product space F := n i=1 F i , e Qf dµ e -f dµ ≤ 1.

Here, Q is the infimum-convolution operator associated to the quadratic cost function C:

Qf (x) := inf y∈F [f (y) + C(x -y)] , x ∈ F, (2)
with C(z) := 1 4 n i=1 z i 2 , z = (z 1 , . . . , z n ) ∈ F. We define Q α , a first variant of the operator Q, which is suitable for some abstract probabilistic situation where metric structure is not a priori provided. Thus, we extend the "convex property (τ )" from product Banach spaces to any product probability spaces (E, E) = n i=1 (E i , E i ). For further measurability considerations, we assume that every singleton {x i }, x i ∈ E i , belongs to the σfield E i . Let µ = µ 1 ⊗ • • • ⊗ µ n be a product probability measure on (E, E). We establish that for every parameter α ∈ (0, 1), and for every bounded measurable function f on E one has

e αQαf dµ 1 α e -(1-α)f dµ 1 1-α ≤ 1.
(3)

This first result is a consequence of a transportation inequality presented further in the introduction (see Theorem 1.2). We observe that (3) will still holds when α equal 0 or 1: when α goes to 1, (3) provides

log e Q1f dµ ≤ f dµ, (4)
and when α goes to 0, it yields log e -f dµ ≤ Q 0 f dµ.

(5)

In order to define Q α , we introduce some notations. For a given measurable space (F, F), let P(F ) be the set of probability measures on F and T (F ) be the set of all transition probabilities from (F, F) to (F, F). For every ζ (1) ∈ P(F ) and p ∈ T (F ), we define ζ (1) p ∈ P(F × F ) and ζ (1) p ∈ P(F ) as follows: for

B ∈ F ⊗ F, ζ (1) p(B) := 1 1 B (x, y)p(x, dy)ζ (1) (dx),
and for A 2 ∈ F, 1) .

ζ (1) p(A 2 ) := p(•, A 2 )dζ (
ζ (1) 2) . We define a cost function c α , α ∈ (0, 1):

c α ( ) := α(1 -) log(1 -) -(1 -α ) log(1 -α ) α(1 -α) , for 0 ≤ ≤ 1,
and c α ( ) := +∞ if > 1. The family of function c α , 0 < α < 1, dominates the usual quadratic cost function: c 1 ( ) ≥ c α ( ) ≥ c 0 ( ) ≥ 2 /2 for every ≥ 0. For every y = (y 1 , . . . , y n ) ∈ E, let y i denote the vector (y 1 , . . . , y i ), 1 ≤ i ≤ n. For every x ∈ E, we define

Q α f (x) := inf p∈T (E)
f (y)p(x, dy)

+ n i=1 c α 1 1 xi =yi p i (x i , dy i |y i-1 ) p(x, dy) , (6)
where the infimum is taken over all transition probabilities p in T (E) that can be written as

p(x, dy) = p 1 (x 1 , dy 1 ) p 2 (x 2 , dy 2 |y 1 ) • • • p n (x n , dy n |y n-1 ), (7)
where p i (•, •|y i-1 ) ∈ T (E i ) and

p i (x i , dy i |y i-1 ) := λ i (y i-1 )δ xi (dy i ) + (1 -λ i (y i-1 ))ν i (dy i ), (8)
with λ i being a measurable function on

i-1

j=1 (E j , E j ) satisfying 0 ≤ λ i (y i-1 ) ≤ 1, and ν i being a probability measure on E i absolutely continuous with respect to µ i (µ i << ν i ), for which ν i ({x i }) = 0. Actually p i (x i , •|y i-1 ) is a convex combination of ν i and the Dirac measure at point x i . Choosing p i (x i , dy i |y i-1 ) = δ xi (dy i ) for all 1 ≤ i ≤ n in the definition (6), we observe that Q α f ≤ f.

For a better understanding, let us consider the one-dimensional definition of

Q α : for n = 1, Q α f (x) := inf p f (y)p(x, dy) + c α 1 1 x =y p(x, dy) , with p(x, dy) := λδ x + (1 -λ)ν, 0 ≤ λ ≤ 1. ( 9
)
Comparing this definition with the usual definition (2), we see that f (y) has been replaced by f (y)p(x, dy), and the cost C(x -y) by c α 1 1 x =y p(x, dy) , the cost to pay to move from the initial position x. This cost only depends on the probability 1 -λ to move from x and is independent of the way to move from x given by ν.

A second result of this paper is a functional inequality of the form (3) with a new operator R α for which the cost term depends on both (1 -λ) and the measure ν. More precisely, in dimension one, we set for every

x ∈ E, R α f (x) := inf p f (y)p(x, dy) + d α 1 1 x =y dp(x, y) dµ(y) dµ(y) (10)
and in dimension n we set

R α f (x) := inf p∈T (E)
f (y)p(x, dy)

+ n i=1 d α 1 1 xi =yi dp i (x i , y i |y i-1 ) dµ i (y i ) dµ i (y i )p(x, dy) , (11)
where the infimum is taken over all transition probabilities p in T (E) such that (7) and (8) hold. The function d α is defined by its convex Legendre-transform

d * α : d α ( ) := sup h∈R [h -d * α (h)] , ∈ R, where d * α (h) := αe (1-α)h + (1 -α)e -αh -1 α(1 -α) , h ∈ R. ( 12
)
We notice that d α ( ) is equivalent to 2 /2 at zero. One has

d 0 ( ) := (1 + ) log(1 + ) -, ≥ 0, d 1/2 ( ) := 2 log 2 + 1 + 2 4 -4 1 + 2 4 -1 , ≥ 0, d 1 ( ) := (1 -) log(1 -) + , if 0 ≤ ≤ 1,
and d 1 ( ) := +∞ if > 1.

Theorem 1.1 For every bounded measurable function f , for every parameter 0 < α < 1, one has

e αRαf dµ 1 α e -(1-α)f dµ 1 1-α ≤ 1. ( 13
)
The proof of Theorem 1.1 is obtained by induction over n. A similar method works for (3). However, as already mentioned, we deduce (3) from the transportation type inequality (14) below. Let us recall the definition of the relative entropy of a measure ν with respect to a measure µ with density dν/dµ:

Ent µ dν dµ := dν dµ log dν dµ dµ.
We define a pseudo-distance C α between probability measures: for ζ (1) , ζ (2) ∈ P(E):

C α (ζ (1) , ζ (2) ) := inf p∈T (ζ (1) ,ζ (2) ) n i=1 c α 1 1 xi =yi p i (x i , dy i |y i-1 ) p(x, dy)ζ (1) (dx),
where the infimum runs over all p in T (ζ (1) , ζ (2) ) satisfying ( 7) and ( 8).

Theorem 1.2 For every ζ (1) and ζ (2) in P(E) absolutely continuous with respect to µ, one has

C α (ζ (1) , ζ (2) ) ≤ 1 α Ent µ dζ (1) dµ + 1 1 -α Ent µ dζ (2) dµ . ( 14
)
The connection between (3) and ( 14) is obtained following the lines of [B-G]. We easily check that

Q α f dζ (1) -f dζ (2) ≤ C α (ζ (1) , ζ (2) ).
Then (3) follows by applying Theorem 1.2 and by choosing ζ (1) and ζ (2) with respective densities dζ (1) dµ = e αQαf e αQαf dµ and dζ (2) dµ = e -(1-α)f e -(1-α)f dµ .

Remarks:

1. Theorem 1.2 improves the transportation inequality of Marton for contracting Markov chains (see [Mar]), when we restrict the study to product probability measures. Marton introduces a distance d 2 defined by 2) ). Then, Marton's inequality follows by optimizing over α in (14).

n d 2 (ζ (1) , ζ (2) ) 2 := inf p∈T (ζ (1) ,ζ (2) ) n i=1 1 1 xi =yi p i (x i , dy i |y i-1 ) 2 p(x, dy)ζ (1) (dx). Since c α ( ) ≥ 2 /2, ≥ 0, one has n d 2 (ζ (1) , ζ (2) ) 2 ≤ 2C α (ζ (1) , ζ (
2. A first byproduct of (3) and Theorem 1.1 is Talagrand's concentration inequalities of the form (1). By monotone convergence, we extend (3) to any real measurable function f on E satisfying e -(1-α)f dµ < ∞. Then by adopting the convention (+∞)×0 ≤ 1, (3) is extended to all f taking values in R∪{+∞}, for which e -(1-α)f dµ < ∞ holds. For every measurable set A ∈ E, let φ A (x) := 0 if x ∈ A and φ A (x) := +∞ otherwise. Then, for every x ∈ E, one has

Q α φ A (x) = inf p,p(x,A)=1 n i=1 c α 1 1 xi =yi p i (x i , dy i |y i-1 ) p(x, dy).
The infimum is taken over all p in T (E) which transport the mass on A, that is: p(x, A) = 1. This defines a pseudo-distance C α between x and A:

C α (x, A) := Q α φ A (x)
. Applying (3) with the function φ A , we get Theorem 4.2.4 of [START_REF] Talagrand | Concentration of measure and isoperimetric inequalities in product spaces[END_REF]:

e αCα(x,A) dµ(x) ≤ 1 µ(A) α 1-α . ( 15
)
The same argument with the functional inequality ( 13) gives

e αDα(x,A) dµ(x) ≤ 1 µ(A) α 1-α .
where

D α (x, A) := inf p,p(x,A)=1 n i=1 d α 1 1 xi =yi dp i (x i , y i |y i-1 ) dµ i (y i ) dµ i (y i )p(x, dy).
For the careful reader of [START_REF] Talagrand | New concentration inequalities in product spaces[END_REF], this new result improves Theorem 4.2 of [START_REF] Talagrand | New concentration inequalities in product spaces[END_REF] and Theorem 1 of [START_REF] Panchenko | A note on Talagrand's concentration inequalities[END_REF].

A simple procedure to reach deviation tails for a given function f is the following. If f is sufficiently "regular", Q α f or R α f can be estimated from below, and (3) or (13) provide bounds for the Laplace transform of f . With the operator Q α , the regularity of f is given by the "best" non-negative measurable functions h

i : E → R + ∪ {+∞}, 1 ≤ i ≤ n such that f (x) -f (y) ≤ n i=1 h i (x)1 1 xi =yi for µ-almost every x, y ∈ E. (16)
On the right-hand side, the weighted Hamming distance measures how far f (x) is from f (y). Let c * α be the Legendre-transform of c α . If ( 16) holds, then

Q α f (x) ≥ f (x) - n i=1 c * α (h i (x)) for µ-almost every x ∈ E. (17)
Actually, the "best" functions h i 's minimize the quantity

n i=1 c * α (h i (x)
). Similarly, with R α , the regularity of f is given by the "best"

h i : E × E i → R + ∪ {+∞}, 1 ≤ i ≤ n such that f (x) -f (y) ≤ n i=1 h i (x, y i )1 1 xi =yi for µ-almost every x, y ∈ E. ( 18
)
This implies: for µ-almost every

x in E R α f (x) ≥ f (x) - n i=1 d * α (h i (x, y i ))dµ i (y i ), (19)
and the "best" h i 's minimize the quantity

n i=1 d * α (h i (x, y i ))dµ i (y i ).
At the beginning of section 3, we show how to derive deviation tails from ( 3) and ( 17) when f represents a supremum of sum of non-negative random variables (see Corollary 3.3). This is a generic simple example. Then, as a main application of this paper, the same procedure with ( 13) and ( 19) provides deviation's bounds when the variables of the sums are not necessarily non-negative (see Corollary 3.5 and Corollary 3.6). In this matter, the entropy method (developed by [D-S], [A-M-S]) has been first used by Ledoux [Led1] and then by several authors to improve the results by Talagrand, and to reach optimal constants (see [Mas], [START_REF] Rio | Inégalités de concentration pour les processus empiriques de classes de parties[END_REF], [START_REF] Rio | Une inégalité de Bennett pour les maxima de processus empiriques[END_REF]). The entropy method is a general procedure that yields deviation inequalities from a logarithmic Sobolev type inequality via a differential inequality on Laplace transforms. This is well-known as the Herbst argument. Our approach is an alternative to the entropy method. For the suprema of empirical sums, Theorem 3.1 and Theorem 3.4 provide exponential deviations tails that improve in some way the ones obtained with the entropy method. To be complete, recall that Panchenko [START_REF] Panchenko | Symmetrization approach to concentration inequalities for empirical processes[END_REF] introduced a symmetrization approach that allows to get other concentration results for suprema of empirical processes from (3).

[B-L-M1], [B-L-M2], [B-B-L-M], [Bous1], [Bous2], [Kle], [K-R], [Led2],
In Section 4.1, we recover from (3) the Talagrand's deviation tails (around the mean instead of the median) for the one-dimensional bin packing problem.

2 Dimension-free functional inequalities for product measures.

The first part of this section is devoted to the proof of Theorem 1.2, and the second part to the proof of Theorem 1.1.

A transportation type inequality.

Theorem 1.2 is obtained by tensorization of the one-dimensional case. Let 2) be two probability measures on (E, E) with respective densities a and b with respect to µ.

ζ (1) = ζ (
For n = 1, the proof of ( 14) is based on the construction of an optimal transition probability p

* ∈ T (ζ (1) , ζ (2) ). Let [c] + = max(c, 0), c ∈ R. For every x ∈ E with a(x) = 0, we set p * (x, •) := λ(x)δ x (•) + (1 -λ(x))ν(•), (20)
with λ(x) := inf 1, b a (x) , and ν :=

[b-a]+ ζ (1) -ζ (2) T V µ, where ζ (1) -ζ (2) T V
denotes the total variation distance between ζ (1) and ζ (2) ,

ζ (1) -ζ (2) T V := 1 2 |a -b|dµ = [a -b] + dµ.
If a(x) = 0, p * (x, •) is any probability measure in P(E).

Lemma 2.1 According to the above definitions, one has

C α (ζ (1) , ζ (2) ) = c α 1 1 x =y p * (x, dy) ζ (1) (dx), and
c α 1 1 x =y p * (x, dy) ζ (1) (dx) ≤ 1 α Ent µ dζ (1) dµ + 1 1 -α Ent µ dζ (2) dµ .
Proof of Lemma 2.1.

For any p ∈ T (ζ (1) , ζ (2) ), we consider the measure ζ := ζ (1) p, ζ ∈ P(E × E). Let ∆ denote the diagonal of E × E and ζ ∆ be the restriction of ζ to ∆, ζ ∆ (A) = ζ(∆ ∩ (A × A)) := ζ({(x, x), x ∈ A}) for every A ∈ E. Since ζ ∈ M(ζ (1) , ζ (2) ), we have ζ ∆ = cµ with c ≤ inf(a, b). Consequently, since c α is increasing on R + , we get c α 1 - inf(a, b) a ζ (1) (dx) ≤ c α 1 1 x =y p(x, dy) ζ (1) (dx).
This lower bound is reached for p * . Indeed, setting

ζ * = ζ (1) p * , from the definition (20) one has ζ * ∆ = inf(a, b)µ. The function c α satisfies for every u, v > 0, c α 1 - v u ≤ 1 αu c 0 (1 -u) + 1 (1 -α)u c 0 (1 -v).
Consequently integrating with respect to ζ (1) = aµ and using the identity

c 0 1 - dζ (i) dµ dµ = Ent µ dζ (i) dµ ,
for i = 0 or 1, we get Lemma 2.1. 2

Let us now consider the n-dimensional case:

(E, E) = n i=1 (E i , E i ) and µ := µ 1 ⊗ • • • ⊗ µ n . Let us first introduce some notations. For every 1 ≤ k ≤ n, let ζ (1)k denote the marginal of ζ (1) on (E k , E k ) := k i=1 (E i , E i ). The density of ζ (1)k with respect to µ 1 ⊗ • • • ⊗ µ k is a k (x k ) := a(x k , y k+1 , . . . , y n )µ k+1 (dy k+1 ) . . . µ n (dy n ). One has ζ (1)n = ζ (1) and a n = a. Let ζ (1) k (•|x k-1 ) denote the probability measure on (E k , E k ) with density a k (x k |x k-1 ) := a k (x k ) a k-1 (x k-1 ) , with respect to µ k . One has ζ (1)k = ζ (1) 1 • • • ζ (1) k . The same notations are still available with ζ (2) with its density b: ζ (2)k = ζ (2) 1 • • • ζ (2) k .
For every vectors x, y in E, let us consider a sequence of transition probabilities

p k (•, •|x k-1 , y k-1 ) ∈ T ζ (1) k (•|x k-1 ), ζ (2) k (•|y k-1 ) , 1 ≤ k ≤ n. De- fine ζ k (dx k , dy k |x k-1 , y k-1 ) := ζ (1) k (dx k |x k-1 ) p k (x k , dy k |x k-1 , y k-1 ), one has ζ k (•, •|x k-1 , y k-1 ) ∈ M ζ (1) k (•|x k-1 ), ζ (2) k (•|y k-1 ) , and ζ k is a transition prob- ability from k-1 i=1 (E i , E i ) × k-1 i=1 (E i , E i ) to (E k , E k ) × (E k , E k ). The marginals of ζ k := ζ 1 • • • ζ k are ζ (1)k and ζ (2)k . Setting p k (x k , dy k ) := p 1 (x 1 , dy 1 ) • • • p k (x k , dy k |x k-1 , y k-1 ), we exactly have ζ k = ν (1)k p k and p k ν (1)k = ν (2)k
. By definition, we say that the transition probability p(x, dy) := p n (x n , dy n ) is a well-linked sequence of the transition probabilities

p k , 1 ≤ k ≤ n. One has p ∈ T (ζ (1) , ζ (2) ).
To get Theorem 1.2, we define p * as a well-linked sequence of optimal transition probabilities

p * k (•, •|x k-1 , y k-1 ) ∈ T ζ (1) k (•|x k-1 ), ζ (2)
k (•|y k-1 ) defined as in ( 20). And we show that

n i=1 c α 1 1 xi =yi p * i (x i , dy i |y i-1 , x i-1 ) p * (x, dy)dζ (1) (dx) ≤ 1 α Ent µ dζ (1) dµ + 1 β Ent µ dζ (2)
dµ .

The left hand-side of this inequality is equal to

n i=1 c α 1 1 xi =yi p * i (x i , dy i |x i-1 , y i-1 ) ζ (1) i (dx i |x i-1 )ζ * i-1 (dx i-1 , dy i-1 ).
Lemma 2.1 ensures that

c α 1 1 xi =yi p * i (x i , dy i |x i-1 , y i-1 ) ζ (1) i (dx i |x i-1 ) ≤ 1 α Ent µi dζ (1) ( • |x i-1 ) dµ i + 1 β Ent µi dζ (2) ( • |y i-1 ) dµ i .
The proof of Theorem 1.2 is ended by integrating this inequality with respect to the coupling measure

ζ * i-1 ∈ M ζ (1)i-1 , ζ (2)
i-1 , and then using the classical tensorization property of entropy

Ent µ dζ (1) dµ i = n i=1 Ent µi dζ (1) ( • |x i-1 ) dµ i ζ (1)i-1 (dx i-1 ).
2.2 A dimension-free functional inequality.

Theorem 1.1 is obtained by induction over the dimension n using a simple contraction argument given at the end of this section.

For n = 1, according to (10), one has for every

x in E, R α f (x) = f (x) -sup 0≤λ≤1 sup ν (f (x) -f (y))(1 -λ) dν dµ (y) -d α (1 -λ) dν dµ (y) dµ(y),
where ν << µ and ν({x}) = 0. Since d * α denotes the Legendre transform of d α , it follows

R α f (x) ≥ f (x) -d * α ([f (x) -f (y)] + ) dµ(y).
Therefore, Theorem 1.1 implies the next statement.

Lemma 2.2 For every bounded measurable function f ,

e αf (x)-α d * α ([f (x)-f (y)]+)dµ(y) dµ(x) 1 α e -(1-α)f dµ 1 1-α ≤ 1. (21)
Actually, for n = 1, Theorem 1.1 is equivalent to this Lemma. To show it, we need the following observation.

Lemma 2.3 Let f be a bounded measurable function on E. There exists a measurable function f on E satisfying f (x) ≤ f (x) and

R α f (x) ≤ f (x) -d * α [ f (x) -f (y)] + dµ(y), for every x ∈ E.
To simplify the notations, let • γ , γ ∈ R denote the L γ (µ)-norm. For a given function f , Lemma 2.3 ensures that there exists f for which

e Rαf α ≤ exp α f (x) -α d * α [ f (x) -f (y)] + dµ(y) dµ(x)
µ ({x ∈ E, f (x) ≤ t}) = λ ({u ∈ [0, 1], g(u) ≤ t}) , t ∈ R,
and therefore, ( 21) is equivalent to

1 0 e αg(t)-α t 0 d * α (g(t)-g(s))ds dt 1 α 1 0 e -(1-α)g(ω) dω 1 1-α ≤ 1, ( 22
)
where g is an increasing right-continuous function.

We observe that d * α is the best non-negative convex function with d * α (0) = 0, satisfying such an inequality. Applying it to the test functions g u (t) := a1 1 (u,1] (t) with a ≥ 0 and 0 ≤ u ≤ 1, we get u + (1 -u)e αa-αd * α (a)

1 α u + (1 -u)e -(1-α)a 1 1-α ≤ 1, 0 ≤ u ≤ 1.
When u → 0, this provides a lower bound for d * α (a) which gives exactly its definition (12). The cost function d α is therefore optimal.

By an elementary approximation argument, we only need to prove (22) for every increasing simple function g with a finite number of values. This is obtained by induction over the number of values of the simple function g. Clearly (22) holds for constant functions g since d * α (0) = 0. Then, we apply the next induction step.

Proposition 2.4 Let g be an increasing simple function on [0, 1] that reach its maximum value on (v, 1]. For u ∈ (v, 1] and a > 0 let g u (t) = a1 1 (u,1] (t), t ∈ [0, 1]. Then, g satisfies ( 22) implies g + g u satisfies ( 22).

The proof of this Proposition is given in Appendix.

2

Proof of Lemma 2.3. It suffices to find a transition probability p satisfying ( 9) and a function f , f ≤ f , for which

f (y)p (x, dy) + d α 1 1 x =y dp (x, y) dµ(y) dµ(y) = f (x) -d * α [ f (x) -f (y)] + dµ(y) (23) Since d * α is convex, the function ψ α : u -→ (d * α ) ([u -f (y)] + ) dµ(y), is increasing. By differentiating, one has (d * α ) (h) = e (1-α)h -e -αh , h ≥ 0. If u → +∞, then ψ α (u) → +∞, and if u → -∞ then ψ α (u) → 0.
Therefore, by continuity of ψ α , there exists a real number f 0 such that

(d * α ) ([f 0 -f (y)] + ) dµ(y) = 1.
Let A := {x ∈ E, f (x) ≤ f 0 }. We define f by f (x) := f (x) if x ∈ A, and f (x) := f 0 otherwise. Clearly, one has f ≤ f . For a given x ∈ E, we define λ(x) by

1 -λ(x) := (d * α ) [ f (x) -f (y)] + dµ(y). Since (d * α ) is increasing and (d * α ) (0) = 0, if x ∈ A then 0 ≤ 1 -λ(x) = (d * α ) ([f (x) -f (y)] + ) dµ(y) ≤ 1, and if x ∈ A then 1 -λ(x) = (d * α ) ([f 0 -f (y)] + ) dµ(y) = 1. Define p (x, dy) = λ(x)δ x (dy) + (1 -λ(x))ν x (dy)
, where ν x is the probability measure with density

dν x dµ (y) = (d * α ) [ f (x) -f (y)] + 1 -λ(x) .
One has ν x ({x}) = 0. With these definitions, the left hand-side of (23) is

f (x) - (f (x) -f (y))(1 -λ(x)) dν x dµ (y) -d α (1 -λ(x)) dν x dµ (y) dµ(y) = f (x) - (f (x) -f (y))(d * α ) [ f (x) -f (y)] + -d α (d * α ) [ f (x) -f (y)] + dµ(y). If x ∈ A, this expression is f (x)- ( f (x) -f (y))(d * α ) [ f (x) -f (y)] + -d α (d * α ) [ f (x) -f (y)] + dµ(y), that is f (x) -d * α [ f (x) -f (y)] + dµ(y).
Next, if x ∈ A, the same equality holds: from the definition of f 0 , one has

f (x) - (f (x) -f (y))(d * α ) [ f (x) -f (y)] + dµ(y) = f (x) -f (x) (d * α ) ([f 0 -f (y)] + ) dµ(y) + f (y)(d * α ) [ f (x) -f (y)] + dµ(y) = f (y)(d * α ) [ f (x) -f (y)] + dµ(y) = f (x) -( f (x) -f (y))(d * α ) [ f (x) -f (y)] + dµ(y).
2 Let us now present the contractivity argument that extends Theorem 1.1 to any dimension n. We sketch the induction step from n = 1 to n = 2. Then it suffices to repeat the same argument.

Let µ = µ 1 ⊗ µ 2 and E = E 1 × E 2 .
We want to show that e Rαf α ≤ e f -β , with β = 1 -α. From definition (11), one has for every (

x 1 , x 2 ) ∈ E, R α f (x 1 , x 2 ) = inf p1 inf p2 f (y 1 , y 2 )p 2 (x 2 , dy 2 |y 1 ) + d α 1 1 x2 =y2 dp 2 dµ 2 (x 2 , y 2 |y 1 )dµ 2 (y 2 ) p 1 (x 1 |dy 1 ) + d α 1 1 x1 =y1 dp 1 dµ 1 (x 1 , y 1 ) dµ 1 (y 1 ) .
We easily check that R α f = R (1) R (2) f , with for g :

E i → R, i = 1, 2, R (i) g(x i ) := inf pi g(y i )p i (x i , dy i ) + d α 1 1 xi =yi dp i dµ i (x i , y i ) dµ i (y i ) . Let • γ(i) , denotes the L γ (µ i )-norm, γ ∈ R. One has e Rαf α = e R (1) R (2) f α(1) α(2)
.

Applying Theorem 1.1 with the measure µ 1 and the function

x 1 → R (2) f (x 1 , x 2 ), we get e Rαf α ≤ e R (2) f -β(1) α(2) = e -βR (2) f dµ 1 -1 β -α β (2)
.

For γ < 0, the function h → h γ , is concave on the set of positive measurable functions h (see [H-L-P]). Therefore, by Jensen it follows that

e Rαf α ≤ e -βR (2) f -α β (2) dµ 1 -1 β = e R (2) f α(2) -β(1)
.

Then, we apply Theorem 1.1 again with the measure µ 2 and the function

x 2 → f (x 1 , x 2 ), we get e Rαf α ≤ e f -β(2) -β(1) = e f -β .
3 Deviation inequalities for suprema of sums of independent random variables.

Let F be a countable set and let (X 1,t ) t∈F , . . . , (X n,t ) t∈F be n independent processes. In this part, we state deviations inequalities for the random variable

Z := sup t∈F n i=1 X i,t .
It is enough to consider a finite set F = {1, . . . , N }. The results settled in this section then extend to any countable set by monotone convergence. One has Z = f (X), where X = ((X 1,t ) t∈F , . . . , (X n,t ) t∈F ) and f (x

) := sup 1≤t≤N n i=1 x i,t , for x = (x 1 , . . . , x n ) ∈ R N n with x i = (x i,t ) 1≤t≤N . For a given x in R N n , if τ (x) := inf t ∈ F, f (x) = n i=1 x i,t , then f (x) = n i=1 x i,τ (x) . Since f (y) ≥ n i=1 y i,τ (x) , x, y ∈ R N n , we notice that f (x) -f (y) ≤ n i=1 (x i,τ (x) -y i,τ (x) )1 1 xi =yi .
The condition (18) holds setting h i (x, y i ) = [x i,τ (x) -y i,τ (x) ] + . If all x i,t 's and y i,t 's are non-negative, then the condition ( 16) holds with h i (x) = x i,τ (x) . These observations are the key point to reach deviation results for Z from (3) and Theorem 1.1 using the estimates ( 17) and ( 19).

We first consider the sums of non-negative random variables. Applying ( 4) and ( 5) to λf , λ ≥ 0, and using the estimate (17), we get the following basic results.

Theorem 3.1 Assume X i,t ≥ 0 for every i and t, then for every λ ≥ 0, log

E n i=1 (1 + λX i,τ (X) ) ≤ λE[Z], (24)
and

log E [exp(-λZ)] ≤ n i=1 E exp(-λX i,τ (X) ) -1 . (25)
As first application of this Theorem, we recover the well-known bound for the Laplace transform of Z when the X i,t 's belongs to [0, 1] (see [Mas]).

Corollary 3.2 Assume X i,t ∈ [0, 1] for every i and t, then for every

λ ∈ R log E [exp(λZ)] ≤ E[Z](e λ -1).
Proof . For λ < 0, the above inequality is an obvious consequence of (25). Next, for λ > 0, by convexity of the exponential function, we observe that

log E [exp(λZ)] ≤ log E n i=1 (1 + (e λ -1)X i,τ (X) ) .
Then we apply inequality (24).

2 For any r ∈ (1, 2], define

Σ r := n i=1 X r i,τ (X) ≤ sup t∈F n i=1 X r i,t .
When the X i,t 's are not upper-bounded, Theorem 3.1 induces the following new results.

Corollary 3.3 If E[Σ r ] < ∞ for some r ∈ (1, 2]
, then one has: for any u ≥ 0,

P Z ≥ E[Z] + u + u r Σ r E[Σ r ] -1 ≤ exp - r -1 r u r E[Σ r ] 1/(r-1) , (26)
and

P [Z ≤ E[Z] -u] ≤ exp - r -1 r u r E[Σ r ] 1/(r-1) . (27)
Proof . Since e -x + x -1 ≤ (x r /r) for x ≥ 0 and r ∈ (1, 2], the inequality (25) implies: for any λ ≥ 0,

log E [exp(-λ(Z -E[Z]))] ≤ λ r r E[Σ r ].
The inequality (27) then follows from Chebyshev inequality by optimizing over all λ ≥ 0.

From the inequality (24), using x -log(1 + x) ≤ (x r /r), x ≥ 0, r ∈ (1, 2], we also get: for any λ > 0,

log E exp λZ - λ r r Σ r ≤ λE[Z].
Therefore, by Chebyshev inequality, for any v ≥ 0,

P Z ≥ E[Z] + λ r-1 r Σ r + v λ ≤ exp(-v).
When Σ r is close to its mean E[Σ r ], the optimal choice for λ is given by rv = (r -1)λ r E[Σ r ]. Then, the proof of ( 26) is complete by taking

λ r-1 = u/E[Σ r ]. 2
The non-negativity of the X i,t 's is a strong restriction that can be relaxed with Theorem 1.1. Using the estimate (19), it gives: for any α ∈ (0, 1), β = 1-α,

E e αλZ-αS α,λ 1 α E e -βλZ 1 β ≤ 1, λ ≥ 0, (28)
with

S α,λ := n i=1 E d * α λ[X i,τ (X) -X i,τ (X) ] + X ,
where (X i,t ) t∈F is an independent copy of X i = (X i,t ) t∈F and E[ • |X] is the conditional expectation given X = (X 1 , . . . , X n ). By Chebychev inequality, it follows that for every m ∈ R, v ≥ 0,

P Z ≥ m + v λ + S α,λ λ 1-α P [Z ≤ m] α ≤ e -α(1-α)v , (29)
Consequently, it suffices to control the fluctuation of S 2 α,λ to derive deviation bounds for Z. For r ∈ (1, 2], we define

V r := n i=1 E [X i,τ (X) -X i,τ (X) ] r + X .
By a Taylor expansion, S α,λ is of order λ 2 2 V 2 as λ → 0. V 2 is a variance factor, one has

V 2 ≤ V := sup t∈F n i=1 E [X i,t -X i,t ] 2 + X i,t .
When α or β goes to zero, (28) gives the following main estimates for the Laplace transform of Z.

Theorem 3.4 For any λ ≥ 0, one has

log E [exp(λZ -S 1,λ )] ≤ λE [Z] , (30)
and

log E [exp(-λ(Z -E[Z]))] ≤ E [S 0,λ ] . (31)
For any random variable Y , define the ψ 1 -norm by

Y ψ1 := inf {c > 0 | E [exp(|Y |/c)] ≤ 2} . Corollary 3.5 1. If E[V r ] < ∞ for some r ∈ (1, 2],
then for any u ≥ 0,

P Z ≥ E[Z] + u + u r V r E[V r ] -1 ≤ exp - r -1 r u r E[V r ] 1/(r-1) . (32) 2. If C i := [X i,τ (X) -X i,τ (X) ] + ψ1
< ∞ for every i, then for any u ≥ 0,

P [Z ≤ E[Z] -u] ≤ exp - u 2 4C 2 e -1 + M u , (33)
with C 2 := n i=1 C 2 i and M := max(C i ).

Proof . The inequality (32) is a consequence of (30) using d * 1 (x) ≤ (x r /r), x ≥ 0, r ∈ (1, 2], and then following the proof of ( 26).

Let Y i := [X i,τ (x) -X i,τ (X) ] + , for any λ ≥ 0,

E [d * 0 (λY i )] = +∞ 0 C i λ e Ciλu -1 P [Y i /C i ≥ u] du. If 0 ≤ λ < 1/C i , then one has sup u≥0 e Ciλu -1 e u = C i λ exp (1 -C i λ) log(1 -C i λ) C i λ ≤ e -1 C i λ 1 -C i λ ,
and therefore

E [d * 0 (λY i )] ≤ e -1 C 2 i λ 2 1 -C i λ +∞ 0 e u P [Y i /C i ≥ u] du = e -1 C 2 i λ 2 1 -C i λ .
Consequently, for any 0 ≤ λ < 1/M ,

E [S 0,λ ] ≤ e -1 C 2 λ 2 1 -M λ .
Finally, (33) follows using (31), applying Chebychev inequality and then optimizing over all 0 ≤ λ < 1/M . 2 The main interest of (29) and Corollary 3.5 is that no boundedness condition are needed on the X i,t 's. If they are bounded, the results can be refined as follows.

Corollary 3.6 1. Assume that X i,t ≤ M i,t , and E (M i,t -X i,t ) 2 ≤ 1, for every i and t, then for any u ≥ 0,

P(Z ≥ E[Z] + u) ≤ exp   - u 2 1 + ε u E[V ] log 1 + u E[V ]   , (34) ≤ exp - u 2 2E[V ] + 2u , with ε(u) := d * 1 (log (1 + u)) log (1 + u) .
2. Assume that m i,t ≤ X i,t ≤ M i,t , with M i,t -m i,t = 1 for every i and t, then for any u ≥ 0,

P(Z ≤ E[Z] -u) ≤ exp -E[V 2 ]d 0 u E[V 2 ] , (35)
≤ exp - u 2 2E[V 2 ] + 2 3 u , with d 0 (u) := (1 + u) log(1 + u) -u.

Remarks:

1. These results extend the well-known Bennett's or Bernstein's inequalities (see [Ben], [Bern]). These inequalities apply under the assumption of Theorem 3.6 when F = {t 0 }, Z = n i=1 X i,t0 , and m i,t0 ≤ X i,t0 ≤ M i,t0 , with M i,t0m i,t0 = 1. The Bennett's inequality asserts

P(Z ≥ E[Z] + u) ≤ exp -Var(Z)d 0 u Var(Z) , u ≥ 0,
or equivalently

P(Z ≤ E[Z] -u) ≤ exp -Var(Z)d 0 u Var(Z) , u ≥ 0. If F = {t 0 } then E[V ] = E[V 2 ] = Var(Z)
. Therefore, (35) exactly recover the above Bennett inequality. From the Cramér-Chernoff large deviation theorem, we know that d 0 is an optimal rate function, reached when Z approximates a Poisson distribution (for more details, see [Mas]). In the right-hand side deviation inequality (34), since ε(u) → 0 as u → 0, the rate of deviation is optimal in the moderate deviations bandwidth, and since ε(u) → 1 as u → +∞, it is suboptimal in the large deviations bandwidth.

2. If m i,t ≤ X i,t ≤ M i,t with M i,t -m i,t = 1, then by using the classical symmetrization and contraction argument of Ledoux and Talagrand (see [L-T] Lemma 6.3 and Theorem 4.12, see also Lemma 14 [Mas]), one can show that

E[V ] ≤ σ 2 + 16 E sup t∈F n i=1 (X i,t -E[X i,t ]) , (36)
with σ 2 := sup t∈F n i=1 Var(X i,t ). Therefore E[V ] is often close to the maximal variance σ 2 .

3. If the X i,t 's are moreover centered and

Z = sup t∈F | n i=1 X i,t |, then (36) yields E[V ] ≤ σ 2 + 16 E[Z].
Therefore, the results of Corollary 3.6 are of the same nature as the recent results of Klein and Rio [K-R]. A basic difference is that we do not assume that the X i,t 's are centered. Moreover the proofs to reach deviation bounds from Theorem 1.1 are simple, especially for the left-hand side deviations, as regard to the entropy method used in [K-R].

4. We observe that other left-hand side deviation's bounds for centered X i,t 's in (-∞, 1] were given by Klein [Kle] under additional moment conditions.

Proof of Corollary

3.6. If m i,t ≤ X i,t ≤ M i,t with M i,t -m i,t = 1, then S 2 0,λ ≤ d * 0 (λ)V 2 , λ ≥ 0. Consequently (31) gives log E [exp(-λ(Z -E[Z]))] ≤ d * 0 (λ)E [V 2 ] .
Then (35) follows by Chebyshev inequality and optimizing over all λ ≥ 0. Finally, it is well-known that (35) implies the Bernstein inequality.

For the right-hand side deviations, by Hölder inequality and using (30), we get for any q > 1 and 1/p + 1/q = 1

E [exp (λ(Z -E[Z]))] ≤ E exp p q S 2 1,λq 1 p , λ ≥ 0. ( 37
) Define S 2 1,λ := sup t∈F n i=1 Y i,t (λ) with Y i,t (λ) := E d * 1 λ[X i,t -X i,t ] + X i,t . We observe that the function x → d * 1 ( √ x) is concave, therefore since X i,t ≤ M i,t and E (M i,t -X i,t ) 2 ≤ 1, by Jensen inequality, Y i,t (λ) ≤ d * 1 λ M i,t -X i,t 2 ≤ d * 1 (λ),
where . 2 denotes the L 2 -norm. Consequently, applying Corollary 3.2 with S 2 1,λ , by homogeneity, we obtain

log E exp p q S 2 1,λq ≤ log E exp p q S 2 1,λq ≤ e p q d * 1 (λq) -1 d * 1 (λq) E S 2 1,λq .
Then from (37) and since

S 2 1,λ ≤ λ 2 2 V , we get log E [exp (λ(Z -E[Z]))] ≤ q(q -1)λ 2 e d * 1 (λq) q-1 -1 2d * 1 (λq) E [V ] . (38)
For each λ > 0, let q(λ) denote the largest q ≥ 1 so that H(q) = 0 with

H(q) = d * 1 (λq) -q(q -1)λ.
Such a q(λ) exists since H(1) ≥ 0 and H(q) → -∞ as q → +∞. For any λ ≥ 0, (

) ensures that log E [exp (λ(Z -E[Z]))] ≤ λ 2 e q(λ)λ -1 E[V ], 38
and by Chebyshev inequality

P(Z ≥ E[Z] + u) ≤ exp -uλ + λ 2 e q(λ)λ -1 E[V ] u ≥ 0. ( 39
)
One has λq(λ)(q(λ) -1) = d * 1 (λq(λ)) ≤ q(λ) 2 d * 1 (λ), and therefore

λq(λ) 1 - d * 1 (λ) λ ≤ λ. Since d * 1 (λ) λ
→ 0, as λ → 0, it follows that λq(λ) → 0 as λ → 0, and we easily check that λq(λ) → +∞ as λ → +∞. Consequently, for any u ≥ 0 there exists λ u such that u = e q(λu)λu -1

E[V ]. Hence, q(λ u )λ u = log 1 + u E[V ] , and
H(q(λ u )) = 0 gives q(λ u ) = 1 + ε u E[V ]
. The proof of (34) ends by choosing λ = λ u in (39). Then the Bernstein inequality easily follows. 2

4 The one-dimensional bin packing problem.

The one-dimensional bin-packing problem can be described as follows: given a finite collection of independent random weights X 1 , • • • , X n and a collection of identical bins with capacity C (which exceeds the largest of the weights), what is the minimum number of bins N = N (X) (X = {X i , 1 ≤ i ≤ n}) into which the weights can be placed without exceeding the bin capacity C. To simplify we choose C = 1 so that 0

≤ X i ≤ 1, 1 ≤ i ≤ n. Theorem 4.1 Denote Σ 2 = n i=1 X 2 i .
Then for every u ≥ 0

P(N ≥ E[N ] + 1 + u) ≤ exp - u 2 16E[Σ 2 ] + 3u , (40)
and

P(N ≤ E[N ] -1 -u) ≤ exp - u 2 4E[Σ 2 ] . (41)
Observe that for u > n -1, P(|N -E[N ]| ≥ 1 + u) = 0 and one early result [R-T], using martingales approach, is

P(|N -E[N ]| ≥ u) ≤ 2 exp - 2u 2 n , u ≥ 0.
Therefore, the variance factor E[Σ 2 ] improves the factor n for u ≤ E[Σ 2 ]. Obviously, the sum of the variance of the X i 's is expected instead of E[Σ 2 ]. Up to constants, our result is Talagrand's Theorem 6.1 [START_REF] Talagrand | Concentration of measure and isoperimetric inequalities in product spaces[END_REF], where the mean is replaced by a median of N . Talagrand's inequality with the median is a consequence of (15). Lugosi [Lug] also obtains these inequalities from the entropy method, by first recovering some Talagrand's convex distance inequalities (see [Lug]).

Proof . We first show that for any x = {x 1 , . . . , x n }, y = {y 1 , . . . , y n },

N (x) -N (y) ≤ 1 + 2 n i=1
x i 1 1 xi =yi .

To see this, observe that N (x) ≤ N (y)+N ({x i , x i = y i }), then that N ({x i , x i = y i }) ≤ 1 + 2 n i=1 x i 1 1 xi =yi (where . denotes the integer part), since among 2 n i=1 x i 1 1 xi =yi bins, at least one of them is half-empty. Then, we deduce that for any x = {x 1 , . . . , x n }, λ ≥ 0,

Q α (λN )(x) ≥ λN (x) + λ + n i=1 c * α (2λx i ) ≥ λN (x) + λ + 2λ 2 n i=1 x 2 i .
Following the lines of the proof of ( 27), we get (41). For the proof of (40), inequality (3) for α = 1 combined with Cauchy-Schwarz inequality implies

E [exp (λ(N -E[N ]) -λ)] ≤ E exp Σ 2 1,4λ 1 2 , λ ≥ 0,
where Σ 2 1,λ := n i=1 c * 1 (λX i ). By convexity and since 0 ≤ X i ≤ 1, = (1 -u)φ(x)e -α∆(u,x) ((1 -u)e -βx + A(u)) B(u) + (1 -u)e αa-α∆(u,x) , where φ(x) := -e -βx B(u) + (1 -u)e αa-α∆(u,x) e α∆ (u,x) + (1 -u)e -βx + A(u) 1 -A(u)e βx + B(u)e -αx e αa = -B(u)e -βx e α∆(u,a) + B(u) (1 -u)e -βx + A(u) +A(u)e (α+β)x ue -βx -A(u) Since ∆(u, x) ≥ 0 for every u ≥ 0 and x ≥ 0, it follows: for every x ≥ 0, φ(x) ≤ -B(u) ue -βx -A(u) + A(u)e (α+β)x ue -βx -A(u)

E [exp c * 1 (4λX i )] ≤ 1 + E X 2 i e 4λ -4λ -1 ≤ exp E X 2 i d * 0 (

=

A(u)e (α+β)x -B(u) ue -βx -A(u) .

Observing that for every x ≥ 0, A(u)e (α+β)x -B(u) ≥ A(u) -B(u) ≥ 0, and ue -βx -A(u) ≤ u -A(u) ≤ 0, we conclude that φ(x) ≤ 0 and therefore ψ (x) ≤ 0 for every x ≥ 0. This ends the proof of Proposition 2.4 2

  β = 1 -α, and the proof of Theorem 1.1 for n = 1 is finished since f ≤ f . Proof of Lemma 2.2. Let λ denote the Lebesgue measure on [0, 1]. Let h be the increasing repartition function of f under µ: h(t) := µ ({x ∈ E, f (x) ≤ t}) , t ∈ R. Let g denote the inverse repartition function: g(u) := inf{t, h(t) ≥ u}, u ∈ [0, 1]. Then, one has

e

  To simplify the notations, d * denotes d * α . For every measurable increasing functions g on [0, 1], the left-hand side of (22-β(g(ω)-g(t)) dωα β e -α t 0 d * (g(t)-g(s))ds dt Setting ψ(a) := A(u) + (1 -u)e -βa 1 β B(u) + (1 -u)e αa-α∆(u,a)1 α , this inequality is exactly ψ(a) ≤ ψ(0). Therefore the proof of Proposition 2.4 is complete if ψ is a non-decreasing function.By differentiating, we get for every x ≥ 0, u)e -βx (1 -u)e -βx + A(u) + (1 -u) 1 -A(u)e βx + B(u)e -αx e αa-α∆(u,a) B(u) + (1 -u)e αa-α∆(u,x)

Acknowledgments. The author thanks the referee for useful comments that contribute to extend some applications. The author is also grateful to Michel Ledoux, Mathieu Meyer and Krzysztof Oleszkiewicz for helpful comments.

According to Proposition 2.4, we assume that (22) holds for a given simple increasing function g and we want to show that (22) still holds for g + g u , that is:

) dω e βa + (1 -u)

α β e -α u 0 d * (g(t)-g(s)+a)ds dt.

We set A(u) := u 0 e β(g(1)-g(s)) ds, and B(u) := e -α(g(1)-g(s)) ds. Since g is increasing, and since α and β are non-negative, we observe that

and 22) holds with the function g, one has M (u, 0) + N (u, 0) ≤ 1. Therefore, in order to get Proposition 2.4, it suffices to show that

The monoticity properties of the functions g and d * ensure that for 0 ≤ t ≤ u, Consequently, since β and a are non-negative, one has

.e -α(g(1)-g(t)) e -α u 0 d * (g(1)-g(s))ds dt


Quan-Lin Li

John C S Lui

Yang Wang

⋆ Quan

Lin Li

A Matrix-Analytic Solution for Randomized Load Balancing Models with PH Service Times

In this paper, we provide a matrix-analytic solution for randomized load balancing models (also known as supermarket models) with phase-type (PH) service times. Generalizing the service times to the phase-type distribution makes analysis of the supermarket models more difficult and challenging than that of the exponential service time case which has been extensively discussed in the literature. We describe the supermarket model as a system of differential vector equations, provide a doubly exponential solution to the fixed point of the system of differential vector equations, and analyze the exponential convergence of the current location of the supermarket model to its fixed point.

Introduction

In the past few years, a number of companies (e.g., Amazon, Google and Microsoft) are offering the cloud computing service to enterprises, and many content publishers and application service providers are increasingly using Data Centers to host their services. This emerging computing paradigm allows service providers and enterprises to concentrate on developing and providing the own services/goods without worrying about computing system maintenance or upgrade, and thereby significantly reduce their operating cost. For companies that offer cloud computing service in the data centers, they can take advantage of the variation of computing workloads from these customers and achieve the computational multiplexing gain. One of the important technical challenges that they have to address is how to utilize these computing resources in the data centers efficiently since many of these servers can be used. There is a growing interest to examine simple and robust load balancing strategies to efficiently utilize the computing resource of the server farms.

Randomized load balancing is a simple and effective mechanism to fairly utilize computing resources. It can deliver surprisingly good performance measures such as reducing collisions, waiting times and backlogs. In a supermarket model, each arriving job randomly picks a small subset of servers and examines their instantaneous workload, and the job is routed to the least loaded server. When a job is committed to a server, jockeying is not allowed and each server uses the first-come-first-service (FCFS) discipline to process all jobs, e.g., see Mitzenmacher [START_REF] Mitzenmacher | The power of two choices in randomized load balancing[END_REF]. For the supermarket models, most of recent research applied density dependent jump Markov processes to deal with the simple case with Poisson arrival processes and exponential service times, and illustrated that there exists a fixed point which decreases doubly exponentially. Readers may refer to, such as, a simple supermarket model by [START_REF] Mitzenmacher | The power of two choices in randomized load balancing[END_REF][START_REF] Vvedenskaya | Queueing system with selection of the shortest of two queues: An asymptotic approach[END_REF]; simple variations by [START_REF] Mitzenmacher | Analyses of load stealing models using differential equations[END_REF][START_REF] Mitzenmacher | On the analysis of randomized load balancing schemes[END_REF][START_REF] Mitzenmacher | The power of two random choices: a survey of techniques and results[END_REF][START_REF] Luczak | On the maximum queue length in the supermarket model[END_REF]; load information by [START_REF] Dahlin | Interpreting stale load information[END_REF][START_REF] Mitzenmacher | How useful is old information[END_REF]; fast Jackson network by Martin and Suhov [START_REF] Martin | Fast Jackson networks[END_REF][START_REF] Martin | Point processes in fast Jackson networks[END_REF][START_REF] Suhov | Fast Jackson Networks with Dynamic Routing[END_REF]; and general service times by Bramson, Lu and Prabhakar [START_REF] Bramson | Randomized load balancing with general service time distributions[END_REF]. When the arrival processes or the service times are more general, the available results of the supermarket models are few up to now. The purpose of this paper is to provide a novel approach for studying a supermarket model with PH service times, and show that the fixed point decreases doubly exponentially. Also, the PH approximation of order 2 gains new numerical insights on practical applications of the general service times to supermarket models.

The remainder of this paper is organized as follows. In the next section, we describe the supermarket model with the PH service times as a system of differential vector equations based on density dependent jump Markov processes. In Section 3, we set up a system of nonlinear equations satisfied by the fixed point, provide a doubly exponential solution to the system of nonlinear equations, and compute the expected sojourn time of a tagged arriving customer. In Section 4, we study exponential convergence of the current location of the supermarket model to its fixed point.

Supermarket Model

In this section, we describe a supermarket model with a Poisson arrival process and PH service times as a system of differential vector equations based on density dependent jump Markov processes.

Let us formally describe the supermarket model, which is abstracted as a multi-server multi-queue stochastic system. Customers arrive at a queueing system of n > 1 servers as a Poisson process with arrival rate nλ for λ > 0. The service time of each customer is of phase type with irreducible representation (α, T ) of order m. Each arriving customer chooses d ≥ 1 servers independently and uniformly at random from the n servers, and waits for service at the server which currently contains the fewest number of customers. If there is a tie, servers with the fewest number of customers will be chosen randomly. All customers in every server will be served in the FCFS manner. Please see Figure 1 for an illustration of the supermarket model.

For the supermarket models, the PH distribution allows us to model more realistic systems and understand their performance implication under the randomized load balancing strategy. As indicated in [START_REF] Harchol-Balter | Exploiting process lifetime distributions for dynamic load balancing[END_REF], the process lifetime of many parallel jobs, in particular, jobs to the data centers, tends to be non-exponential. For the PH service time distribution, we use the following irreducible representation: (α, T ) of order m, the row vector α is a probability vector whose jth entry is the probability that a service begins in phase j for 1 ≤ j ≤ m; T is an th entry is denoted by t i,j with t i,i < 0 for 1 ≤ i ≤ m, and t i,j ≥ 0 for 1 ≤ i, j ≤ m and i = j. Let T 0 = -T e 0, where e is a column vector of ones with a suitable dimension in the context. The expected service time is given by 1/µ = -αT -1 e. Unless we state otherwise, we assume that all the random variables defined above are independent, and that the system is operating in the stable region ρ = λ/µ < 1.

We define n (i) k (t) as the number of queues with at least k customers and the PH service time in phase i at time t ≥ 0. Clearly, 0 ≤ n

(i) k (t) ≤ n for k ≥ 1 and 1 ≤ i ≤ m. Let X (0) n (t) = n n = 1,
and k ≥ 1

X (k,i) n (t) = n (i) k (t) n ,
which is the fraction of queues with at least k customers and the service time in phase i at time t ≥ 0. We write

X (k) n (t) = X (k,1) n (t) , X (k,2) n (t) , . . . , X (k,m) n (t) , k ≥ 1,
X n (t) = X (0) n (t) , X (1) n (t) , X (2) n (t) , . . . .

The state of the supermarket model may be described by the vector X n (t) for t ≥ 0. Since the arrival process to the queueing system is Poisson and the service time of each server is of phase type, the stochastic process {X n (t) , t ≥ 0} is a Markov process whose state space is given by

Ω n = { g (0) n , g (1)
n , , g (2) n . . . :

g (0) n = 1, g (k) n ≥ g (k+1) n ≥ 0, and ng (k) n is a vector of nonnegative integers for k ≥ 1}. Let s 0 (n, t) = E X (0) n (t) and k ≥ 1 s (i) k (n, t) = E X (k,i) n (t) .
Clearly, s 0 (n, t) = 1. We write

S k (n, t) = s (1) k (n, t) , s (2)
k (n, t) , . . . , s

(m) k (n, t) , k ≥ 1.
As shown in Martin and Suhov [START_REF] Martin | Fast Jackson networks[END_REF] and Luczak and McDiarmid [4], the Markov process {X n (t) , t ≥ 0} is asymptotically deterministic as n → ∞, this is due to the fact that from a block-structured point of view, the PH distribution can keep many excellent properties of the exponential distribution. Thus lim n→∞ E X (0)

n (t) and lim n→∞ E X (k,i) n
always exist by means of the law of large numbers. Based on this, we write

S 0 (t) = lim n→∞ s 0 (n, t) = 1, for k ≥ 1 s (i) k (t) = lim n→∞ s (i) k (n, t) , S k (t) = s (1)
k (t) , s (2)
k (t) , . . . , s (m) k (t)
and

S (t) = (S 0 (t) , S 1 (t) , S 2 (t) , . . .) . Let X (t) = lim n→∞ X n (t).
Then it is easy to see from the Poisson arrivals and the PH service times that {X (t) , t ≥ 0} is also a Markov process whose state space is given by Ω = g (0) , g (1) , g (2) , . . . :

g (0) = 1, g (k) ≥ g (k+1) ≥ 0 for k ≥ 1 .
If the initial distribution of the Markov process {X n (t) , t ≥ 0} approaches the Dirac delta-measure concentrated at a point g ∈ Ω, then its steady-state distribution is concentrated in the limit on the trajectory S g = {S (t) : t ≥ 0}. This indicates a law of large numbers for the time evolution of the fraction of queues of different lengths. Furthermore, the Markov process {X n (t) , t ≥ 0} converges weakly to the fraction vector S (t) = (S 0 (t) , S 1 (t) , S 2 (t) , . . .), or for a sufficiently small ε > 0,

lim n→∞ P {||X n (t) -S (t) || ≥ ε} = 0,
where ||a|| is the L ∞ -norm of vector a.

In what follows we provide a system of differential vector equations in order to determine fraction vector S (t). To that end, we introduce the Hadamard Product of two matrices A = (a i,j ) and B = (b i,j ) as follows:

A ⊙ B = (a i,j b i,j ) .
Specifically, for k ≥ 2 we have

A ⊙k = A ⊙ A ⊙ • • • ⊙ A k matrix A .
To determine the fraction vector S (t), we need to set up a system of differential vector equations satisfied by S (t). To that end, we provide a heuristic description in terms of a concrete example, and indicate how to derive the differential vector equations.

In the supermarket model with n servers, we determine the expected change in the number of queues with at least k customers over a small time period of length dt. The probability vector that during the time period [0,dt), any arriving customer joins a queue with k -1 customers is given by n λS ⊙d k-1 (n, t) -λS ⊙d k (n, t) dt. Similarly, the probability vector that during the time period [0,dt), a customer leaves a server queued by k customers is given by

n S k (n, t) T + S k+1 (n, t) T 0 α dt. Therefore we can obtain dE [n k (n, t)] =n λS ⊙d k-1 (n, t) -λS ⊙d k (n, t) dt + n S k (n, t) T + S k+1 (n, t) T 0 α dt, which leads to dS k (n, t) dt = λS ⊙d k-1 (n, t) -λS ⊙d k (n, t) + S k (n, t) T + S k+1 (n, t) T 0 α. (1)
Taking n → ∞ in both sides of Equation ( 1), we have

dS k (t) dt = λS ⊙d k-1 (t) -λS ⊙d k (t) + S k (t) T + S k+1 (t) T 0 α. (2)
Using a similar analysis to Equation (2), we obtain a system of differential vector equations for the fraction vector S (t) = (S 0 (t) , S 1 (t) , S 2 (t) , . . .) as follows:

S 0 (t) = 1,

d dt S 0 (t) = -λS d 0 (t) + S 1 (t) T 0 , (3)
d dt S 1 (t) = λαS d 0 (t) -λS ⊙d 1 (t) + S 1 (t) T + S 2 (t) T 0 α, (4)
and for k ≥ 2,

d dt S k (t) = λS ⊙d k-1 (t) -λS ⊙d k (t) + S k (t) T + S k+1 (t) T 0 α. (5)
Note that the above derivations or the ordinary differential equations (ODE) descriptions are the approximation to the large scale stochastic system. One can be made the formal derivation using the "mean field" or the "stochastic density jump Markov process" techniques. While it is in form of ODE, we believe readers can easily understand the ODE of the above equations since the supplementary variable method are extensively used in the study of stochastic systems.

A Matrix-Analytic Solution

In this section, we provide a doubly exponential solution to the fixed point of the system of differential vector equations ( 3), ( 4) and [START_REF] Martin | Point processes in fast Jackson networks[END_REF].

A row vector π = (π 0 , π 1 , π 2 , . . .) is called a fixed point of the fraction vector S (t) if lim t→+∞ S (t) = π. In this case, it is easy to see that

lim t→+∞ d dt S (t) = 0.
Therefore, as t → +∞ the system of differential vector equations ( 3), ( 4) and ( 5) can be simplified as

-λπ d 0 + π 1 T 0 = 0, ( 6
)
λαπ d 0 -λπ ⊙d 1 + π 1 T + π 2 T 0 α = 0, (7)
and for k ≥ 2,

λπ ⊙d k-1 -λπ ⊙d k + π k T + π k+1 T 0 α = 0. ( 8
)
In general, it is more difficult and challenging to express the fixed point of the supermarket models with more general arrival processes or service times. Fortunately, we can derive a closed-form expression for the fixed point π = (π 0 , π 1 , π 2 , ...) for the supermarket model with PH service times by means of a novel matrix-analytic approach.

Noting that S 0 (t) = 1 for all t ≥ 0, it is easy to see that π 0 = 1. It follows from Equation ( 6) that

π 1 T 0 = λ. (9)
To solve Equation ( 9), we denote by ω the stationary probability vector of the irreducible Markov chain T + T 0 α. Obviously, we have

ωT 0 = µ, λ µ ωT 0 = λ. (10)
Thus, we obtain π 1 = λ µ ω = ρ • ω. Based on the fact that π 0 = 1 and

π 1 = ρ • ω, it follows from Equation (7) that λα -λρ d • ω ⊙d + ρ • ωT + π 2 T 0 α = 0, which leads to π 2 T 0 = λρ d ω ⊙d e.
Let θ = ω ⊙d e. Then it is easy to see that θ ∈ (0, 1), and

π 2 T 0 = λθρ d .
Using a similar analysis to Equation (10), we have

π 2 = λθρ d µ ω = θρ d+1 • ω. (11)
Based on π 1 = ρ • ω and π 2 = θρ d+1 • ω, it follows from Equation ( 8) that for k = 2,

λρ d • ω ⊙d -λθ d ρ d 2 +d • ω ⊙d + θρ d+1 • ωT + π 3 T 0 α = 0, which leads to π 3 T 0 = λθ d+1 ρ d 2 +d .
Using a similar analysis on Equation ( 10), we have

π 3 = λθ d+1 ρ d 2 +d µ ω = θ d+1 ρ d 2 +d+1 • ω. (12)
Based on Equations ( 11) and ( 12), we may infer that there is a structured expression

π k = θ d k-2 +d k-3 +•••+d+1 ρ d k-1 +d k-2 +•••+d+1 • ω for k ≥ 1.
To that end, the following theorem states this important result.

Theorem 1 The fixed point π = (π 0 , π 1 , π 2 , . . .) is given by

π 0 = 1, π 1 = ρ • ω and for k ≥ 2, π k = θ d k-2 +d k-3 +•••+1 ρ d k-1 +d k-2 +•••+1 • ω, (13)
or

π k = θ d k-1 -1 d-1 ρ d k -1 d-1 • ω = ρ d k-1 (θρ) d k-1 -1 d-1
• ω.

Proof: By induction, one can easily derive the above result. It is clear that Equation ( 13) is correct for the cases with l = 2, 3 according to Equations [START_REF] Mitzenmacher | The power of two random choices: a survey of techniques and results[END_REF] and [START_REF] Suhov | Fast Jackson Networks with Dynamic Routing[END_REF]. Now, we assume that Equation ( 13) is correct for the cases with l = k. Then it follows from Equation (8) that for l = k + 1, we have

λθ d k-2 +d k-3 +•••+d ρ d k-1 +d k-2 +•••+d • ω ⊙d -λθ d k-1 +d k-2 +•••+d ρ d k +d k-1 +•••+d • ω ⊙d + θ d k-2 +d k-3 +•••+1 ρ d k-1 +d k-2 +•••+1 • ωT + π k+1 T 0 α = 0, which leads to π k+1 T 0 = λθ d k-1 +d k-2 +•••+d+1 ρ d k +d k-1 +•••+d .
By a similar analysis to (10), we have

π k+1 = λθ d k-1 +d k-2 +•••+d+1 ρ d k +d k-1 +•••+d µ ω = θ d k-1 +d k-2 +•••+d+1 ρ d k +d k-1 +•••+d+1 • ω.
This completes the proof. Now, we compute the expected sojourn time T d that a tagged arriving customer spends in the supermarket model. For the PH service times, a tagged arriving customer is the kth customer in the corresponding queue with probability vector π ⊙d k-1 -π ⊙d k . When k ≥ 1, the head customer in the queue has been served, and so its service time is residual and is denoted as X R . Let X be of phase type with irreducible representation (α, T ). Then X R is of phase type with irreducible representation (ω, T ). Clearly, we have

E [X] = α (-T ) -1 e, E [X R ] = ω (-T ) -1 e.
Thus it is easy to see that the expected sojourn time of the tagged arriving customer is given by

E [T d ] = π d 0 -π ⊙d 1 e E [X] + ∞ k=1 π ⊙d k -π ⊙d k+1 e {E [X R ] + kE [X]} = ρ d θ (ω -α) (-T ) -1 e + α (-T ) -1 e 1 + ∞ k=1 θ d k -1 d-1 ρ d k+1 -d d-1 . (14)
When the arrival process and the service time distribution are Poisson and exponential, respectively, it is clear that α = ω = θ = 1 and α (-T )

-1 e = 1/µ, thus we have

E [T d ] = 1 µ ∞ k=0 ρ d k+1 -d d-1
, which is the same as Corollary 3.8 in Mitzenmacher [START_REF] Mitzenmacher | The power of two choices in randomized load balancing[END_REF].

In what follows we consider an interesting problem: When using the PH approximation, how many moments of the service time distribution are needed to obtain a better accuracy for computing the fixed point or the expected sojourn time. It is well-known from the theory of probability distributions that the first three moments are basic for analyzing such an accuracy. Also, we can construct a PH distribution of order 2 by using the first three moments. Telek and Heindl [START_REF] Telek | Matching moments for acyclic discrete and continuous phasetype distributions of second order[END_REF] provided a fitting procedure for matching a PH distribution of order 2 in terms of the first three moments exactly. It is necessary to list the fitting procedure as follows:

For a nonnegative random variable X, let m n = E [X n ], n ≥ 1. We take a PH distribution of order 2 with the canonical representation (α, T ), where α = (η, 1 -η) and

T = -ξ 1 ξ 1 0 -ξ 2 , 0 ≤ η ≤ 1 and 0 < ξ 1 ≤ ξ 2 .
Note that the three unknown parameters η, ξ 1 and ξ 2 can be obtained from the first three moments m 1 , m 2 and m 3 .

In Table 1, c 2 X = m 2 m 2 1 -1 which is the squared coefficient of variation. If the moments do not satisfy these conditions in Table 1, then we may analyze the following four cases:

(a.1) if m 2 < 1.5m 2 1 , then we take m 2 = 1.5m 2 1 ; (a.2) if 0.5 ≤ c 2 X ≤ 1 and m 3 < 3m 3 1 3c 2 X -1 + √ 2 1 -c 2 X 3 2
, then we take

m 3 = 3m 3 1 3c 2 X -1 + √ 2 1 -c 2 X 3 2
;

Table 1. Specific Bounds of the First Three Moments

Moment Condition Bounds

m1 0 < m1 < ∞ m2 1.5m 2 1 ≤ m2 m3 0.5 ≤ c 2 X ≤ 1 3m 3 1 " 3c 2 X -1 + √ 2 `1 -c 2 X ´3 2 " ≤ m3 ≤ 6m 3 1 c 2 X 1 < c 2 X 3 2 m 3 1 `1 + c 2 X ´2 < m3 < ∞ (a.3) if 0.5 ≤ c 2 X ≤ 1 and m 3 > 6m 3 1 c 2 X , then we take m 3 = 6m 3 1 c 2 X ; and (a.4) if 1 < c 2 X and m 3 ≤ 3 2 m 3 1 1 + c 2 X 2
, then we take

m 3 = 3 2 m 3 1 1 + c 2 X 2 . Let c = 3m 2 2 -2m 1 m 3 , d = 2m 2 1 -m 2 , b = 3m 1 m 2 -m 3 and a = b 2 -6cd.
If the first three moments satisfy their specific bounds shown in Table 1 or the exceptive four cases, then the three unknown parameters η, ξ 1 and ξ 2 can be computed in the following three cases.

(1) If c > 0, then

η = -b + 6m 1 d + √ a b + √ a , ξ 1 = b - √ a c , ξ 2 = b + √ a c .
(

) If c < 0, then η = b -6m 1 d + √ a -b + √ a , ξ 1 = b + √ a c , ξ 2 = b - √ a c . 2
(3) If c = 0, then

η = 0, ξ 1 > 0, ξ 2 = 1 m 1 .
From the above discussion, we can construct a PH distribution of order 2 to approximate an arbitrarily given general distribution under the same first three moments. In fact, such an approximation achieves a better accuracy in computations of the fixed point and the expected sojourn time.

For the PH distribution of order 2, we have

T + T 0 α = -ξ 1 ξ 1 0 -ξ 2 + 0 ξ 2 η 1 -η = -ξ 1 ξ 1 ξ 2 η -ξ 2 η , which leads to ω = ξ 2 η ξ 1 + ξ 2 η , ξ 1 ξ 1 + ξ 2 η and θ = ξ d 1 + ξ d 2 η d (ξ 1 + ξ 2 η) d .
Note that the PH distributions are dense in the set of all nonnegative random variables, we can numerically provide necessary understanding for the role played by the general service times in performance analysis of the supermarket model by means of the PH approximation of order 2. At the same time, all computations involved in the PH approximation of order 2 are very simple to implement.

In the remainder of this section, we provide two examples to illustrate that our approach is effective and efficient in analyzing the doubly exponential solution of supermarket models with non-exponential service requirements. Table 2 shows how the doubly exponential solution (π 1 to π 4 ) depends on the vectors α (1) and α (2), respectively. Let λ = 1, µ = 2.7500 and d = 2. Figure 2 provides a comparison for the expected sojourn times between the exponential service time and the PH service time when they have a same expected service time. This figure shows that the PH service time makes the lower expected sojourn time.

Exponential convergence to the fixed point

In this section, we study exponential convergence of the current location S (t) of the supermarket model to its fixed point π.

We provide some notation for comparison of two vectors. Let a = (a 1 , a 2 , a 3 , . . .)

and b = (b 1 , b 2 , b 3 , . . .). We write a ≺ b if a k < b k for some k ≥ 1 and a l ≤ b l for l = k, l ≥ 1; and a b if a k ≤ b k for all k ≥ 1.
The following proposition analyzes how the current location S (t) of the supermarket model can be affected by the initial point S (0) for t > 0.

Proposition 1 If S (0) S (0), then S (t) S (t). It follows from the system of differential vector equations ( 3), ( 4) and ( 5)

that for d ≥ 2 d dt S (t) = S ⊙d (t) Q 1 + S (t) Q 2 , (15)
where

Q 1 =      -λI λI -λI λI -λI λI . . . . . .      , Q 2 =      0 T 0 α T T 0 α T . . . . . .      . It follows from (15) that d dt S ⊙(1-d) (t) = eQ 1 + S ⊙(1-d) (t) Q 2 . Let W (t) = S ⊙(1-d) (t). Then d dt W (t) = eQ 1 + W (t) Q 2 .
Hence we obtain

W (t) = W (0) + e t 0 exp {Q 2 x} dx exp {Q 1 t} . Since W (t) = S ⊙(1-d) (t), we have W (0) = S ⊙(1-d) (0)
and

S (t) = S ⊙(1-d) (t) ⊙ 1 1-d = [W (t)] ⊙ 1 1-d .
Therefore, the solution to the system of differential vector equations ( 3), ( 4) and ( 5) is given by

S (t) = S ⊙(1-d) (0) + e t 0 exp {Q 2 x} dx exp {Q 1 t} ⊙ 1 1-d . (16)
Note that t 0 exp {Q 2 x}dx > 0 and exp {Q 1 t} > 0, it is easy to see from (16) that if S (0) S (0), then S (t) S (t) for t > 0. This completes the proof.

Based on Proposition 1, the following theorem shows that the fixed point π is an upper bound of the current location S (t) for all t ≥ 0.

Theorem 2 For the supermarket model, if there exists some k such that S k (0) = 0, then the sequence {S k (t)} has an upper bound sequence which decreases doubly exponentially for all t ≥ 0, that is, S (t) π for all t ≥ 0.

Proof: Let S k (0) = π k for k ≥ 1. Then for each k ≥ 1, S k (t) = S k (0) = π k for all t ≥ 0, since S (0) = S 1 (0) , S 2 (
0) , S 2 (0) , . . . is a fixed point in the supermarket model. If S k (0) = 0 for some k, then S k (0) ≺ S k (0) and S j (0) S j (0) for j = k, j ≥ 1, thus S (0) ≺ S (0). It is easy to see from Proposition 1 that S k (t) S k (t) = π k for all k ≥ 1 and t ≥ 0. Thus we obtain that for all k ≥ 1 and t ≥ 0

S k (t) ≤ θ d k-1 -1 d-1 ρ d k -1 d-1 • ω.
This completes the proof.

To show the exponential convergence, we define a Lyapunov function Φ (t) as

Φ (t) = ∞ k=1 w k [π k -S k (t)] e
in terms of the fact that S k (t) π k for k ≥ 1 and π 0 = S 0 (t) = 1, where {w k } is a positive scalar sequence with w k+1 ≥ w k ≥ w 1 = 1 for k ≥ 2.

The following theorem illustrates that the distance between the fixed point and the current location quickly comes close to zero with exponential convergence.

Theorem 3 For t ≥ 0, Φ (t) ≤ c 0 e -δt , where c 0 and δ are two positive constants. In this case, the potential function Φ (t) is exponentially convergent.

Proof: Note that Then it follows from (17) that In this case, it is easy to see that

Φ (t) = ∞ k=1 w k [π k -S k (t)
d dt Φ (t) = -{[(w 2 -w 1 )] c 1 (t) + w 1 [f 1 (t) -1] d 1 (t)} • [π 1 -S 1 (t)] e
w 2 ≥ 1 + δ + 1 -f 1 (t) c 1 (t)
and for k ≥ 2

w k+1 ≥ w k + δw k c k (t) + d k (t) c k (t) (w k -w k-1
) .

Thus we have

d dt Φ (t) ≤ -δ ∞ k=0
w k [π k -S k (t)] e = -δΦ (t) ,

Fig. 1 .

 1 Fig. 1. The supermarket model: each customer can probe the loading of d servers

  Example one (PH Distribution) Let λ = 1, d = 5, m = 3, α (1) = (1/3, 1/3, 1/3) and α (2) = (1/12, 7/12, 1/3),

  Example two (Expected Sojourn Time) We consider an m-order PH distribution with irreducible representation (α, T ), where m = 2, α = (1/2, 1/2) and

Fig. 2 .

 2 Fig. 2. The expected sojourn times correspond to the exponential and PH distributions

  k+1 -w k ) c k (t) + (w k-1 -w k ) d k (t)] • [π k -S k (t)] e.For a constant δ > 0, we takew 1 = 1, [(w 2 -w 1 )] c 1 (t) + w 1 [f 1 (t) -1] d 1 (t) ≥ δw 1 and (w k+1 -w k ) c k (t) + (w k-1 -w k ) d k (t) ≥ δw k .

Table 2 .

 2 The doubly exponential solution depends on the vector α

		α = ( 1 3 , 1 3 , 1 3 )	α = ( 1 12 , 7 12 , 1 3 )
	π1	(0.0741, 0.1358 , 0.2346)	(0.0602, 0.1728, 0.2531)
	π2 (5.619e-05, 1.030e-05, 1.779e-04 ) (7.182e-05, 2.063e-04, 3.020e-04)
	π3 (1.411e-20, 2.587e-20, 4.469e-20) (1.739e-19, 4.993e-19, 7.311e-19)
	π4 (1.410e-98, 2.586e-98, 4.466e-98) (1.444e-92, 4.148e-92, 6.074e-92)

  ] e, -w 1 [λS d 0 (t) α -λS ⊙d 1 (t) + S 1 (t) T + S 2 (t) T 0 α]e (t) -λS ⊙d k (t) + S k (t) T + S k+1 (t) T 0 α]e.By means of S 0 (t) = 1 and T e = -T 0 , we can obtaind dt Φ (t) = -w 1 [λ -λS ⊙d 1 (t) e -S 1 (t) T 0 + S 2 (t) T 0 ] (t) e -λS ⊙d k (t) e -S k (t) T 0 + S k+1 (t) T 0 ]. (17)We take some nonnegative constants c k (t) and d k (t) for k ≥ 1 such thatλ = f 1 (t) S 1 (t) T 0 , for k ≥ 1 λS ⊙d k (t) e = c k (t) [π k -S k (t)] e and S

	It follows from Equations (3) to (5) that
	d dt	∞ k=1 k=2 Φ (t) = -∞ k-1 -w k [λS ⊙d w k [λS ⊙d k-1
	we have	d dt	Φ (t) = -	∞ k=1	w k	d dt	S k (t) e.

k (t) T 0 = d k (t) [π k -S k (t)] e.

Acknowledgements

The work of Quan-Lin Li was supported by the National Science Foundation of China under grant No. 10871114, John C.S. Lui was supported by the RGC grant, and Yang Wang was supported by the National Science Foundation of China under grant No. 61001075.

which can leads to

This completes the proof.