Approximating the Limiting Quicksort Distribution

James Allen Fill^{1}

Department of Mathematical Sciences

The Johns Hopkins University

            jimfill@jhu.edu and http://www.mts.jhu.edu/~fill/

and Svante Janson

Department of Mathematics

Uppsala University

       svante.janson@math.uu.se and http://www.math.uu.se/~svante/

ABSTRACT

The limiting distribution of the normalized number of comparisons used by Quicksort to sort an array of n numbers is known to be the unique fixed point with zero mean of a certain distributional transformation S. We study the convergence to the limiting distribution of the sequence of distributions obtained by iterating the transformation S, beginning with a (nearly) arbitrary starting distribution. We demonstrate geometrically fast convergence for various metrics and discuss some implications for numerical calculations of the limiting Quicksort distribution. Finally, we give companion lower bounds which show that the convergence is not faster than geometric.

AMS 2000 subject classifications. Primary 68W40; secondary 68P10, 60E05, 60E10, 60F05.

Key words and phrases. Quicksort, characteristic function, density, moment generating function, sorting algorithm, coupling, Fourier analysis, Kolmogorv–Smirnov distance, total variation distance, integral equation, numerical analysis, d_{p}-metric.

Date. January 15, 2001.

1Research supported by NSF grant DMS–9803780, and by the Acheson J. Duncan Fund for the Advancement of Research in Statistics.

1 Introduction and summary

The Quicksort algorithm of Hoare [9] is “one of the fastest, the best-known, the most generalized, … and the most widely used algorithms for sorting an array of numbers” [4]. Quicksort is the standard sorting procedure in Unix systems, and in a special issue of Computing in Science & Engineering, guest editors Jack Dongarra and Francis Sullivan ([3]; see also [10]) chose Quicksort as one of the ten algorithms “with the greatest influence on the development and practice of science and engineering in the 20th century.” Our goal in this introductory section is to review briefly some of what is known about the analysis of Quicksort and to summarize how this paper advances that analysis.

The Quicksort algorithm for sorting an array of n numbers is extremely simple to describe. If n = 0 or n = 1, there is nothing to do. If n ≥ 2, pick a number uniformly at random from the given array. Compare the other numbers to it to partition the remaining numbers into two subarrays. Then recursively invoke Quicksort on each of the two subarrays.

Let X_{n} denote the (random) number of comparisons required (so that X_{0} = 0). Then X_{n} satisfies the distributional recurrence relation

X_{n} =^{L} X_{U_{n}−1} + X_{n−U_{n}}^{∗} + n − 1, n ≥ 1,

where =^{L} denotes equality in law (i.e., in distribution), and where, on the right, U_{n} is distributed uniformly on the set {1, …, n}, X_{j}^{∗} =^{L} X_{j}, and

U_{n}; X_{0}, …, X_{n−1}; X_{0}^{∗}, …, X_{n−1}^{∗}

are all independent.

As is well known and quite easily established, for n ≥ 0 we have

μ_{n} := E X_{n} = 2(n + 1)H_{n} − 4n ∼ 2n ln n,

where H_{n} := ∑_{k=1}^{n} k^{−1} is the nth harmonic number and ∼ denotes asymptotic equivalence. It is also routine to compute explicitly the standard deviation of X_{n} (see Exercise 6.2.2-8 in [12]), which turns out to be ∼ n\sqrt{7 − \frac{2}{3}π^{2}}.

Consider the normalized variate

Y_{n} := (X_{n} − μ_{n})/n,  n ≥ 1.  (1.1)

Régnier [14] showed using martingale arguments that Y_{n} → Y in distribution, with Y satisfying the distributional identity

Y =^{L} UY + (1 − U)Y^{∗} + g(U) =: h_{Y,Y^{∗}}(U),  (1.2)

where

g(u) := 2u ln u + 2(1 − u) ln(1 − u) + 1,  (1.3)

and where, on the right of =^{L} in (1.2), U, Y, and Y^{∗} are independent, with Y^{∗} =^{L} Y and U ∼ unif(0, 1). Rösler [15] showed that (1.2) characterizes the limiting law L(Y), in the precise sense that F := L(Y) is the unique fixed point of the operator

G = L(V) ↦ SG := L(UV + (1 − U)V^{∗} + g(U))  (1.4)

(in what should now be obvious notation) subject to

E V = 0, Var V < ∞.

[The fixed points of G with finite mean are the translates L(Y + c) with c constant, but there are other fixed points without mean; see [7] for a complete characterization.]

Rösler [15] showed that the moment generating function of the limiting distribution L(Y) is everywhere finite. We have studied the limiting distribution further in [6], showing that L(Y) has a density f which is infinitely differentiable, and that each derivative f^{(k)}(y) is bounded and decays as y → ±∞ more rapidly than any power of |y|^{−1}. (This improves an earlier result by Tan and Hadjicostas [18].)

The purpose of the present paper is to study the convergence to the limiting distribution L(Y) of the sequence of distributions obtained by iterating Rösler’s operator S in (1.4), beginning with a (nearly) arbitrary starting distribution. To fix notation, we let Z_{0} be an arbitrary random variable, and F_{0} := L(Z_{0}) its distribution. We define, for n ≥ 1,

Z_{n} := h_{Z_{n−1},Z_{n−1}^{∗}}(U),

with Z_{n−1}^{∗} =^{L} Z_{n−1} and Z_{n−1}, Z_{n−1}^{∗}, and U independent; in other words,

F_{n} := L(Z_{n}) = S^{n}F_{0}, n ≥ 0.

Let ∥X∥_{2} := (E X^{2})^{1/2} denote the L^{2}-norm, and let d_{2} denote the metric on the space of probability distributions with finite variance defined by

d_{2}(G_{1}, G_{2}) := min ∥X_{1} − X_{2}∥_{2},  (1.5)

taking the minimum over all pairs of random variables X_{1} and X_{2} (defined on the same probability space) with L(X_{1}) = G_{1} and L(X_{2}) = G_{2}. Note that, using the coupling with X_{1} and X_{2} independent, for any G_{1} and G_{2} each with zero mean and finite variance,

d_{2}(G_{1}, G_{2}) ≤ (E X_{1}^{2} + E X_{2}^{2})^{1/2} ≤ ∥X_{1}∥_{2} + ∥X_{2}∥_{2},  (1.6)

when L(X_{1}) = G_{1} and L(X_{2}) = G_{2}. Rösler [15] showed that if Z_{0} has mean 0 and finite variance, then F_{n} → F in the d_{2}-distance with a geometric rate:

d_{2}(F_{n}, F) ≤ (2/3)^{n/2}d_{2}(F_{0}, F) ≤ (2/3)^{n/2}(Var Z_{0} + σ^{2})^{1/2},  (1.7)

where

σ^{2} := Var Y = 7 − \frac{2}{3}π^{2} ≐ 0.42.  (1.8)

Our main interest is to show similar estimates for other measures of the distance between F_{n} and F.

We will show in Section 3, using estimates of the characteristic functions given in [6] and Section 2, that the distribution F_{n} has a bounded, continuous density function f_{n}, at least as soon as n ≥ 3, and that, if Z_{0} has mean 0 and finite variance, then f_{n} converges uniformly to f, with a geometric rate of convergence, as n → ∞. We further show geometrically fast convergence in the total variation and Kolmogorov–Smirnov distances, too.

In Section 4 we give bounds for the moment generating functions of Y and of Z_{n}. In Section 5 we show that if Z_{0} has mean 0 and a finite moment generating func-

tion ψ_{0}, then the moment generating function ψ_{n} of F_{n} is finite and converges uniformly on compact intervals to the moment generating function of Y, again with a geometric rate of convergence. We study in particular the cases Z_{0} = 0 and Z_{0} normally distributed with zero mean and sufficiently large variance; it turns out that in these cases ψ_{n}(λ) converges monotonically.

In Section 6, we discuss some implications for numerical calculations of the limiting Quicksort distribution F, showing how explicit and arbitrarily small error bounds can be obtained.

Finally, in Sections 7–8 we give some companion lower bounds, showing that the convergence is not faster than geometrical for several different metrics. We also show geometrically fast convergence in the d_{p} metric for any finite p.

Remark 1.1. The mode and rate of convergence of the distribution of the actual normalized Quicksort variables Y_{n} of (1.1) to the limit F is a quite different matter, which will be studied in another paper [8].

2 Bounds on the characteristic functions

In [6] we gave bounds on the characteristic function of Y. The same method yields, more generally, bounds on the characteristic function of Z_{n} for arbitrary Z_{0}. We write ϕ_{X}(t) := E e^{itX} for any random variable X.

Theorem 2.1. For every real p ≥ 0 there is a constant 0 < c_{p} < ∞ such that for any Z_{0} and any n > p + 1, the characteristic function ϕ_{Z_{n}}(t) satisfies

|ϕ_{Z_{n}}(t)| ≤ c_{p}|t|^{−p} for all t ∈ R.  (2.1)

The best possible constants c_{p} satisfy c_{0} = 1, c_{1/2} ≤ 2, c_{3/4} ≤ \sqrt{8π}, c_{1} ≤ 4π, c_{3/2} < 187, c_{5/2} < 103215, c_{7/2} < 197102280, and the relation

c_{p+1} ≤ 2^{p+1}c_{p}^{1+(1/p)}p/(p − 1),  p > 1;  (2.2)

moreover, at least if we restrict (2.1) to n ≥ p + 2,

c_{p} ≤ 2^{p^{2}+6p},  p > 0.  (2.3)

[The bounds on the constants c_{p} obtained here are the same as for the special case Z_{0} =^{L} Y (whence Z_{n} =^{L} Y for every n) in [6]. However, there is no reason to believe that our method yields the best possible bounds, and the best constants for the special case in [6] may be smaller than the best constants in Theorem 2.1 here.]

Proof. The proof is almost identical to the proof of the special case in [6], so we will omit some details. For any random variable Z, we abuse notation slightly and denote by SZ the random variable h_{Z,Z^{∗}}(U) = UZ + (1 − U)Z^{∗} + g(U) where U, Z, and Z^{∗} are independent, with Z^{∗} =^{L} Z and U ∼ unif(0, 1); thus SZ is a random variable with the distribution SL(Z). By conditioning on U, we obtain the fundamental relation

ϕ_{SZ}(t) = ∫_{0}^{1}ϕ_{Z}(ut) ϕ_{Z}((1 − u)t) e^{itg(u)} du,  t ∈ R,  (2.4)

and thus the estimate

|ϕ_{SZ}(t)| ≤ ∫_{0}^{1}|ϕ_{Z}(ut)| |ϕ_{Z}((1 − u)t)| du.  (2.5)

To complete the proof, we give a series of lemmas.

Lemma 2.2. For any real numbers y and z, the random variable h_{y,z}(U) defined by (1.2) satisfies

|E e^{ith_{y,z}(U)}| ≤ 2|t|^{−1/2}.

Proof. This follows by a method of van der Corput [2, 13, 6], using little more than the fact that h_{y,z} is convex with h′′_{y,z} ≥ 8 on (0, 1).

Lemma 2.3. For any random variable Z and real t, we have |ϕ_{SZ}(t)| ≤ 2|t|^{−1/2}.

Proof. Lemma 2.2 yields

|ϕ_{SZ}(t)| = |E e^{ith_{Z,Z^{∗}}(U)}| ≤ E|E(e^{ith_{Z,Z^{∗}}(U)} | Z, Z^{∗})| ≤ 2|t|^{−1/2}.

Returning to our sequence (Z_{n}), the preceding lemma applies to all elements except Z_{0}, i.e.,

|ϕ_{Z_{n}}(t)| ≤ 2|t|^{−1/2},  n ≥ 1,  (2.6)

which yields the case p = 1/2 of Theorem 2.1. We improve the exponent by induction, using (2.5).

Lemma 2.4. Let 0 < p < 1. If |ϕ_{Z}(t)| ≤ c_{p}|t|^{−p}, t ∈ R, then

|ϕ_{SZ}(t)| ≤ \frac{[Γ(1 − p)]^{2}}{Γ(2 − 2p)} c_{p}^{2}|t|^{−2p}.

Proof. By (2.5) and the hypothesis,

|ϕ_{SZ}(t)| ≤ ∫_{0}^{1}c_{p}^{2}|ut|^{−p}|(1 − u)t|^{−p} du = c_{p}^{2}|t|^{−2p} ∫_{0}^{1}u^{−p}(1 − u)^{−p} du,

and the result follows by evaluating the beta integral.

In particular, using (2.6), Lemma 2.4 yields

|ϕ_{Z_{n}}(t)| ≤ 4π|t|^{−1},  n ≥ 2.  (2.7)

This proves (2.1) for p = 1, with c_{1} ≤ 4π. Since |ϕ_{Z_{n}}(t)| ≤ 1, for any p ≤ 1 we trivially have |ϕ_{Z_{n}}(t)| ≤ |ϕ_{Z_{n}}(t)|^{p}, which by (2.7) establishes (2.1) for all p ≤ 1 with c_{p} ≤ (4π)^{p}; applying Lemma 2.4 again, we obtain (2.1) for all p < 2. Somewhat better numerical bounds are obtained for 1/2 < p < 1 by taking a geometric average between the cases p = 1/2 and p = 1; this yields c_{p} ≤ 2^{2p}π^{2p−1}, 1/2 ≤ p ≤ 1. In particular, we have c_{3/4} ≤ \sqrt{8π}, and thus, by Lemma 2.4, c_{3/2} ≤ 8π^{1/2}[Γ(1/4)]^{2} < 186.4 < 187.

Lemma 2.5. Let p > 1. If |ϕ_{Z}(t)| ≤ c_{p}|t|^{−p}, t ∈ R, then

|ϕ_{SZ}(t)| ≤ 2^{p+1}c_{p}^{1+(1/p)} \frac{p}{p − 1}|t|^{−(p+1)}.

Proof. This is similar to the proof of Lemma 2.4, substituting the hypothesis (and the trivial |ϕ_{Z}| ≤ 1) into (2.5), but the estimate of the integral is slightly more complicated; for details see [6].

Lemma 2.5 completes, by induction, the proof of (2.1) and the estimate (2.2).

The bound for c_{3/2} obtained above and (2.2) now yield (using Maple) first c_{5/2} < 103215 and then c_{7/2} < 197102280. These bounds and (2.2) further yield

c_{p} ≤ 2^{p^{2}+5p},  p = k + \frac{3}{2},  (2.8)

for integers k ≥ 0; again see [6] for details. To obtain (2.3) if p > 1/2, let p_{1} := ⌈p−\frac{1}{2}⌉+\frac{1}{2}. Then, by (2.1) and (2.8), provided n ≥ p + 2 > p_{1} + 1,

|ϕ_{Z_{n}}(t)|^{1/p} ≤ |ϕ_{Z_{n}}(t)|^{1/p_{1}} ≤ 2^{p_{1}+5}|t|^{−1} ≤ 2^{p+6}|t|^{−1}.

The case p ≤ 1/2 follows similarly from (2.6), which completes the proof of Theorem 2.1.

Remark 2.6. A variety of other bounds are possible. For example, if we begin with the inequality (2.7) and use (2.5), we can easily derive the following result:

|ϕ_{Z_{n}}(t)| ≤ \frac{32π^{2}}{t^{2}} (ln(\frac{t}{4π} ) + 2) ≤ \frac{32π^{2} ln t}{t^{2}} for all t ≥ 1.72 and n ≥ 3.  (2.9)

3 Convergence of densities

It is easily checked that the random variable h_{y,z}(U) is absolutely continuous for every fixed y and z, and thus, by mixing, SZ is absolutely continuous for every Z. In other words, for any Z_{0}, the random variables Z_{n} have densities for all n ≥ 1; cf. [18]. These densities may be unbounded and discontinuous, at least for n = 1, as is seen in the case Z_{0} ≡ 0. However, we now can show that for n ≥ 3, at least, no such irregularities occur.

Theorem 3.1. If n ≥ 3, then Z_{n} has a bounded continuous density function f_{n}, for any Z_{0}. More generally, if k ≥ 0, then f_{n} is k times continuously differentiable for all n ≥ k +3, and there exists a constant C_{k} independent of Z_{0} and n (with n ≥ k +3) such that |f_{n}^{(k)}(x)| ≤ C_{k}, x ∈ R. Explicitly, |f_{n}(x)| ≤ 16 when n ≥ 5, and |f′_{n}(x)| ≤ 2466 when n ≥ 6.

Proof. Theorem 2.1 shows, in particular, that as soon as n ≥ 3,

|ϕ_{Z_{n}}(t)| ≤ min(1, 187|t|^{−3/2}),

and thus ϕ_{Z_{n}} is integrable. This implies, as is well-known (see e.g., [5, Theorem XV.3.3]) that Z_{n} has a bounded continuous density f_{n} given by the Fourier inversion formula

f_{n}(x) = \frac{1}{2π} ∫_{−∞}^{∞}e^{−itx} ϕ_{Z_{n}}(t) dt,  x ∈ R.  (3.1)

Moreover, using Theorem 2.1 with p = k + \frac{3}{2}, we see that t^{k}ϕ_{Z_{n}}(t) is also integrable when n ≥ k + 3, which by a standard argument shows that f_{n} is k times differentiable, with

f_{n}^{(k)}(x) = \frac{1}{2π} ∫_{−∞}^{∞}(−it)^{k} e^{−itx} ϕ_{Z_{n}}(t) dt,  x ∈ R;  (3.2)

and thus

sup_{x} |f_{n}^{(k)}(x)| ≤ \frac{1}{2π} ∫_{−∞}^{∞}|t|^{k} |ϕ_{Z_{n}}(t)| dt,  (3.3)

where the latter integral can be estimated using Theorem 2.1 with p = k + \frac{3}{2}. The argument above yields the bound

|f_{n}(x)| ≤ \frac{1}{2π} ∫_{−∞}^{∞}min(1, 187|t|^{−3/2}) dt = \frac{3}{π} 187^{2/3} < 31.3,  n ≥ 3.  (3.4)

To obtain better numerical bounds we combine Theorem 2.1 for p = 0, 1/2, 3/2, 1, 5/2, 7/2 and (2.9) (for t in different intervals; see [6] for details); this yields, provided n ≥ 5, f_{n}(x) ≤ \frac{1}{2π} ∫ |ϕ_{n}| < 15.3; similarly, invoking also (2.1) with p = 9/2, f′_{n}(x) ≤ \frac{1}{2π} ∫ |t||ϕ_{n}(t)| dt < 2465.9 for n ≥ 6.

Theorem 3.2. Suppose that E Z_{0} = 0 and Var Z_{0} < ∞. Then the density functions f_{n} of Theorem 3.1 converge uniformly to the (smooth) density function f of Y at a geometric rate:

sup_{x} |f_{n}(x) − f(x)| = O(r^{n}) for every fixed r > (2/3)^{1/2}.

Explicitly, for any p > 1 and n > p + 1,

sup_{x} |f_{n}(x) − f(x)| ≤ \frac{A}{2π} (\frac{2c_{p}}{A} )^{2/(p+1)} \frac{p + 1}{p − 1}(\frac{2}{3})^{( \frac{1}{2} − \frac{1}{p+1} )n},  (3.5)

where A := (Var Z_{0} + σ^{2})^{1/2} and c_{p} is as in Theorem 2.1. In particular,

sup_{x} |f_{n}(x) − f(x)| ≤ 2297A(\frac{2}{3})^{5n/18} < 2297A(0.8935)^{n},  n ≥ 5.  (3.6)

Moreover,

sup_{x} |f_{n}(x) − f(x)| ≤ \frac{128A}{π}(\frac{2}{3})^{(n/2)−3.7\sqrt{n}},  n ≥ 3.  (3.7)

Proof. By the Fourier inversion formula (3.1),

|f_{n}(x) − f(x)| ≤ \frac{1}{2π} ∫_{−∞}^{∞}|ϕ_{Z_{n}}(t) − ϕ_{Y} (t)| dt.  (3.8)

In order to estimate the right hand side, note that for any random variables X and Y,

|ϕ_{X}(t) − ϕ_{Y} (t)| ≤ E|e^{itX} − e^{itY} | ≤ E|tX − tY | ≤ |t| ∥X − Y ∥_{2};

since the characteristic functions here depend on the marginal distributions only, this and the definition (1.5) yield

|ϕ_{X}(t) − ϕ_{Y} (t)| ≤ |t|d_{2}(L(X), L(Y)).

In particular, with δ_{n} := d_{2}(F_{n}, F),

|ϕ_{Z_{n}}(t) − ϕ_{Y} (t)| ≤ |t|δ_{n}.  (3.9)

Further, for any p > 1 and n > p + 1, Theorem 2.1 yields the estimate

|ϕ_{Z_{n}}(t) − ϕ_{Y} (t)| ≤ |ϕ_{Z_{n}}(t)| + |ϕ_{Y} (t)| ≤ 2c_{p}|t|^{−p}.

Consequently, for any T > 0,

∫_{−∞}^{∞}|ϕ_{Z_{n}}(t) − ϕ_{Y} (t)| dt ≤ ∫_{−T}^{T} δ_{n}|t| dt + ∫_{|t|>T} 2c_{p}|t|^{−p} dt

= δ_{n}T^{2} + 4\frac{c_{p}}{p − 1}T^{1−p}.

For given n and p, the optimal choice here is T := (2c_{p}/δ_{n})^{1/(p+1)}, giving the bound

∫_{−∞}^{∞}|ϕ_{Z_{n}}(t) − ϕ_{Y} (t)| dt ≤ \frac{p + 1}{p − 1}(2c_{p})^{2/(p+1)}δ_{n}^{1−(2/(p+1))}.  (3.10)

With (3.8) and the estimate (1.7), this yields (3.5). Choosing p = 7/2 and evaluating the constants numerically, using A ≥ σ > 0.648, we obtain (3.6).

To obtain the final estimate, we use (2.3) and observe that, for p ≥ 2,

(\frac{2c_{p}}{A} )^{2/(p+1)} ≤ (\frac{2^{p^{2}+6p+1}}{σ} )^{2/(p+1)} = 2^{2(p+5)}(\frac{2^{−4}}{σ} )^{2/(p+1)} ≤ 2^{2(p+5)}(1 − \frac{2}{p + 1}),

which by (3.5) yields that for n ≥ p + 2 ≥ 4,

sup_{x} |f_{n}(x) − f(x)| ≤ \frac{A}{2π} 2^{2p+10}(\frac{2}{3})^{( \frac{1}{2} − \frac{1}{p+1} )n}.

Choosing the optimal p := [n ln(3/2)/(2 ln 2)]^{1/2} − 1, we find (3.7) [with the constant (8(ln 2)/ ln(3/2))^{1/2} < 3.69812 multiplying \sqrt{n}], at least when n ≥ 31. For 3 ≤ n ≤ 30, (3.7) follows trivially from (3.4), since the right hand side of (3.7) then is larger than

To test out Theorem 3.2 numerically, choose Z_{0} ≡ 0, so that A = σ ≐ 0.648. For n = 100, (3.6) yields the bound 0.0192; for n ≥ 177, (3.7) is better, and yields for example 3.21 × 10^{−6} for n = 177, 2.07 × 10^{−6} for n = 180, and 1.07 × 10^{−7} for n = 200.

Remark 3.3. Similarly, using (3.2), we obtain geometric uniform convergence of the first derivatives, and of any higher derivatives, of the density functions.

Remark 3.4. Suppose that Z_{0} has finite moments of all orders. Then, by Lemma 7.2 below, E|Z_{n}|^{p} is finite and stays bounded in n, for each real 0 ≤ p < ∞. It follows that the characteristic functions ϕ_{Z_{n}} are infinitely differentiable with derivatives bounded uniformly in n. If we apply both Theorem 2.1 and (3.9) to |ϕ_{Z_{n}}(t) − ϕ_{Y} (t)| and take the geometric mean of the resulting bounds, we find, for n > 2p + 2,

|ϕ_{Z_{n}}(t) − ϕ_{Y} (t)| ≤ [2c_{2p+1}|t|^{−2p}d_{2}(F_{n}, F)]^{1/2} .

It follows easily by induction on k, using [6, Lemma 2.10], that in fact, for every real p ≥ 0 and integer k ≥ 0, there is a constant c_{p,k} [depending on L(Z_{0})] such that for all n > 2^{k+1}p + 2 we have, with ρ_{k} := (2/3)^{2^{−k−2}} < 1,

sup_{t∈R} |t|^{p} |ϕ_{Z_{n}}^{(k)}(t) − ϕ_{Y}^{(k)}(t)| ≤ c_{p,k}ρ_{k}^{n}.

Omitting details, since the Fourier transform is continuous on the Schwartz space [17]

S := {f : sup_{t} |t|^{p}|f^{(k)}(t)| < ∞ for all p, k ≥ 0},

it follows that for each k and p, |x|^{p}f_{n}^{(k)}(x) converges uniformly to |x|^{p}f^{(k)}(x) with geometric rate.

Theorem 3.2 treats uniform approximation of f by f_{n}, using the norm ∥f_{n} − f∥_{∞} := sup_{x} |f_{n}(x) − f(x)|. We now turn to studying the error in the L^{1}-norm ∥f_{n} − f∥_{1} := ∫_{−∞}^{∞}|f_{n} − f|.

Note first that, because ∫_{−∞}^{∞}f_{n} = ∫_{−∞}^{∞}f = 1,

\frac{1}{2} ∫_{−∞}^{∞}|f_{n}(x) − f(x)| dx = ∫_{−∞}^{∞}(f(x) − f_{n}(x))^{+} dx,

and that this coincides with the total variation distance

d_{TV}(F_{n}, F) := sup_{A⊆R} |P(Z_{n} ∈ A) − P(Y ∈ A)|;  (3.11)

moreover, it dominates the Kolmogorov–Smirnov distance

d_{KS}(F_{n}, F) := sup_{x∈R} |P(Z_{n} ≤ x) − P(Y ≤ x)| ≤ d_{TV}(F_{n}, F).  (3.12)

Theorem 3.5. Suppose that E Z_{0} = 0 and Var Z_{0} < ∞. Then the total variation and Kolmogorov–Smirnov distances between F_{n} and F converge geometrically to 0: d_{KS}(F_{n}, F) ≤ d_{TV}(F_{n}, F) = O(r^{n}) for every fixed r > (2/3)^{1/2}. Explicitly, for any n ≥ 1,

d_{KS}(F_{n}, F) ≤ d_{TV}(F_{n}, F) ≤ 135An(\frac{2}{3})^{(n/2)−3.7\sqrt{n}}.  (3.13)

Proof. For any a ∈ (0, 1),

d_{TV}(F_{n}, F) = ∫_{−∞}^{∞}(f(x) − f_{n}(x))^{+} dx ≤ ∥f_{n} − f∥_{∞}^{1−a} ∫_{−∞}^{∞}f(x)^{a} dx,  (3.14)

where ∥f_{n} − f∥_{∞} is estimated in Theorem 3.2. The final integral can be estimated by Hölder’s inequality: for any b > 0

\begin{aligned}∫_{−∞}^{∞}f(x)^{a} dx & = ∫_{−∞}^{∞}f(x)^{a}e^{ab|x|} · e^{−ab|x|} dx \\ & ≤ (∫_{−∞}^{∞}f(x)e^{b|x|} dx)^{a}(∫_{−∞}^{∞}e^{−ab|x|/(1−a)} dx)^{1−a} \\ & ≤ [ψ(b) + ψ(−b)]^{a} (\frac{2(1 − a)}{ab})^{1−a} \\ & = \frac{2}{ab}(ab\frac{ψ(b) + ψ(−b)}{2})^{a}(1 − a)^{1−a}, \end{aligned}  (3.15)

where ψ(λ) := E e^{λY} is the moment generating function of Y. Rösler [15] proved that ψ(λ) is finite for all λ; thus ∫ f^{a} < ∞ for every a ∈ (0, 1), and the first claim follows by (3.14) and Theorem 3.2.

For (3.13) we choose b = 1/3, for which it will be shown in Theorem 4.1 below that ψ(±b) ≤ exp(1/9) < 1.2, and thus (3.15) implies ∫_{−∞}^{∞}f^{a} < 2/(ab) = 6/a. Denoting the right hand side of (3.7) by B, we thus obtain from (3.14) and (3.7), observing that B ≥ (3/2)^{−n/2},

d_{TV}(F_{n}, F) ≤ \frac{6}{a}B^{1−a} ≤ \frac{6}{a}(3/2)^{an/2}B.

We optimize by taking a := 2/(n ln(3/2)) and obtain the following bound (for n ≥ 5, so that a < 1; smaller n are trivial since d_{TV} ≤ 1):

d_{TV}(F_{n}, F) ≤ 3enB ln(3/2) = \frac{384e ln(3/2)}{π}An(\frac{2}{3})^{(n/2)−3.7\sqrt{n}}.

Remark 3.6. If we are content with a weaker explicit bound, we can avoid invoking estimates of ψ by using moments of Y instead. For example,

∫_{−∞}^{∞}f(x)^{1/2} dx ≤ (∫_{−∞}^{∞}f(x)(σ^{2} + x^{2}) dx)^{1/2}(∫_{−∞}^{∞} \frac{dx}{σ^{2} + x^{2}} )^{1/2} = (2πσ)^{1/2} < 2.1

and thus

d_{KS}(F_{n}, F) ≤ d_{TV}(F_{n}, F) ≤ 2.1∥f_{n} − f∥_{∞}^{1/2}.

4 Bounds on moment generating functions

Letting ψ_{Z}(λ) := E e^{λZ} denote the moment generating function of a random variable Z, we find in analogy with (2.4) the relation

ψ_{SZ}(λ) = ∫_{0}^{1}ψ_{Z}(uλ) ψ_{Z}((1 − u)λ) e^{λg(u)} du,  λ ∈ R.  (4.1)

In particular, it follows that if ψ_{Z}(λ) is finite for all λ, then so is ψ_{SZ}(λ).

Rösler [15] proved that the moment generating function ψ_{Y} is everywhere finite and that for every L ≥ 0 there is a constant K_{L} such that

ψ_{Y} (λ) ≤ e^{K_{L}λ^{2}},  |λ| ≤ L.  (4.2)

Moreover, it is implicit in the proof that

if ψ_{Z}(λ) ≤ e^{K_{L}λ^{2}} for |λ| ≤ L, then ψ_{SZ}(λ) ≤ e^{K_{L}λ^{2}} for |λ| ≤ L.  (4.3)

Note that (4.3) implies by induction that if we choose Z_{0} ≡ 0, then ψ_{Z_{n}}(λ) ≤ e^{K_{L}λ^{2}}, |λ| ≤ L, for every n, and thus (4.2) follows by Fatou’s lemma. More generally, if (4.3) holds and ψ_{Z_{0}}(λ) ≤ e^{K_{L}λ^{2}}, |λ| ≤ L, then by induction ψ_{Z_{n}}(λ) ≤ e^{K_{L}λ^{2}}, |λ| ≤ L, for every n.

Rösler did not give explicit values of the constants K_{L}, but such values can be obtained from his proof as follows. [Actually, Rösler [15] treated the somewhat more complicated case of the variables Y_{n} of (1.1); see [8] for explicit constants in that case. In our case there are some simplifications leading to better constants. Moreover, we introduce some deviations from Rösler’s proof designed to improve our bounds.]

Theorem 4.1. Let L_{0} ≐ 5.018 be the largest root of e^{L} = 6L^{2}. Then (4.2) and (4.3) hold with

K_{L} = \begin{cases}1, & L ≤ 0.42, \\ 12, & 0.42 < L ≤ L_{0}, \\ 2L^{−2}e^{L}, & L_{0} < L, \end{cases}

or any larger number. In particular, we can always take K_{L} = max(2L^{−2}e^{L}, 12).

For λ ≤ 0, we can obtain much better estimates. [For (4.3), we restrict to λ ≤ 0 in both the assumption and the conclusion.]

Theorem 4.2. We have (4.2) and (4.3) for λ ≤ 0 with

K_{L} = \begin{cases}0.5, & L ≤ 0.62, \\ 1.25, & 0.62 < L, \end{cases}

or any larger number. In particular, we can always take K_{L} = 1.25 for λ ≤ 0.

Proof of Theorems 4.1 and 4.2. If ψ_{Z}(λ) ≤ e^{Kλ^{2}} for |λ| ≤ L, then by (4.1), for |λ| ≤ L,

ψ_{SZ}(λ) ≤ ∫_{0}^{1} e^{Kλ^{2}[u^{2}+(1−u)^{2}]+λg(u)} du = e^{Kλ^{2}}Ee^{λg(U)−2Kλ^{2}U(1−U)}.

Hence, (4.3) holds with K_{L} = K if (and only if)

f_{K}(λ) := Ee^{λg(U)−2Kλ^{2}U(1−U)} ≤ 1,  when |λ| ≤ L.  (4.4)

Similarly, (4.3) holds with K_{L} = K for λ ≥ 0 (respectively, for λ ≤ 0) if (4.4) holds for 0 ≤ λ ≤ L (resp., for −L ≤ λ ≤ 0). Clearly, f_{K}(λ) decreases as K increases, and thus if some K satisfies (4.4), then so does any larger K.

Following Rösler, we argue differently for small and large L in order to find a K satisfying (4.4). For small L we use a Taylor expansion. By straightforward differentiations,

f_{K}(0) = 1,
f′_{K}(0) = Eg(U) = 0,
f′′_{K}(0) = E(g(U)^{2} − 4KU(1 − U)) = \frac{1}{3}σ^{2} − \frac{2}{3}K,
f′′′_{K}(λ) = E[((g(U) − 4KλU(1 − U))^{3} − 12KU(1 − U)(g(U) − 4KλU(1 − U)))
× exp(λg(U) − 2Kλ^{2}U(1 − U))].

We write the last formula as f′′′_{K}(λ) = E[X(U, λ)] and note that 0 ≤ U(1 − U) ≤ 1/4 and −η ≤ g(U) ≤ 1, where

η := −g(\frac{1}{2}) = 2 ln 2 − 1 ≐ 0.386.

Consider first λ ≥ 0. By Taylor’s formula, for 0 ≤ λ ≤ L,

\begin{aligned}f_{K}(λ) & ≤ 1 + \frac{1}{2}λ^{2}f′′_{K}(0) + \frac{1}{6}λ^{3} sup_{0≤ξ≤L} f′′′_{K}(ξ) \\ & ≤ 1 + \frac{1}{6}λ^{2}(σ^{2} − 2K + L sup_{0≤ξ≤L} f′′′_{K}(ξ)) \end{aligned}

so (4.4) is satisfied for λ ≥ 0 provided

L sup_{0≤ξ≤L} f′′′_{K}(ξ) ≤ 2K − σ^{2}.  (4.5)

If g(U) ≥ 0, we find

X(U, λ) ≤ (1 + 3K^{2}L)e^{L}, 0 ≤ λ ≤ L;

while if g(U) ≤ 0, we find

X(U, λ) ≤ 3K(η + KL), 0 ≤ λ ≤ L.

For K ≥ 1, in either case, because 3η > 1,

X(U, λ) ≤ (3Kη + 3K^{2}L)e^{L}, 0 ≤ λ ≤ L,

and thus

L sup_{0≤ξ≤L} f′′′_{K}(ξ) ≤ L(3Kη + 3K^{2}L)e^{L}.

It is readily checked that this is less than 2K − σ^{2} so that (4.5) holds, for K = 1 and L = 0.42.

For larger L, we begin by another crude estimate. Let W =^{L} U/2 be uniformly distributed on (0, 1/2). Then, by |g(U)| ≤ 1 and symmetry,

f_{K}(λ) ≤ e^{|λ|}E exp(−2Kλ^{2}U(1 − U)) = e^{|λ|}E exp(−2Kλ^{2}W (1 − W))

≤ e^{|λ|}E exp(−Kλ^{2}W ) = e^{|λ|} ∫_{0}^{1} exp(−Kλ^{2}u/2) du

= e^{|λ|} \frac{1 − exp(−Kλ^{2}/2)}{Kλ^{2}/2} =: g_{K}(λ). (4.6)

Note that that g_{K}, too, decreases if K is increased. Taking the logarithmic derivative, we find for λ > 0,

\begin{aligned}(ln g_{K}(λ))′ & = 1 − \frac{2}{λ} + Kλe^{−Kλ^{2}/2}(1 − exp(−Kλ^{2}/2))^{−1} \\ & = 1 − \frac{2}{λ} + \frac{Kλ}{e^{Kλ^{2}/2} − 1}. \end{aligned}  (4.7)

For λ ≥ 2, this is evidently positive, and thus g_{K} then is increasing. Hence, if K ≥ K̃ := 2L^{−2}e^{L}, then

g_{K}(λ) ≤ g_{K}(L) ≤ gK̃(L) = 1 − exp(−e^{L}) < 1, 2 ≤ λ ≤ L.

For smaller λ, we take K = 12, and check numerically that g_{12}(0.42) < 1. Moreover,

e^{Kλ^{2}/2} − 1 = e^{6λ^{2}} − 1 ≥ 6λ^{2} + 18λ^{4}

and further, if 1/3 ≤ λ ≤ 1,

(1 − \frac{λ}{2})^{−1} ≤ 1 + λ ≤ 1 + 3λ^{2}

and thus

\frac{Kλ}{e^{Kλ^{2}/2} − 1} ≤ \frac{12λ}{6λ^{2}(1 + 3λ^{2})} ≤ \frac{2}{λ}(1 − \frac{λ}{2}) = \frac{2}{λ} − 1.

Hence, (4.7) shows that g_{12} is decreasing on [1/3, 1], and thus

g_{12}(λ) ≤ g_{12}(0.42) < 1, 0.42 ≤ λ ≤ 1.

Finally,

g_{12}(λ) ≤ \frac{1}{6} \frac{e^{λ}}{λ^{2}} ≤ \frac{e}{6} < 1, 1 ≤ λ ≤ 2.

Combining these estimates we find that if K ≥ max(12, 2L^{−2}e^{L}), then f_{K}(λ) ≤ g_{K}(λ) < 1 whenever 0.42 ≤ λ ≤ L, while f_{K}(λ) ≤ f_{1}(λ) ≤ 1 for 0 ≤ λ ≤ 0.42, and thus (4.4) holds when λ ≥ 0.

We have also shown that K = 12 will do for L ≤ 2 and λ ≥ 0; since 2L^{−2}e^{L} is increasing for L ≥ 2, and thus less than 12 for 2 ≤ L < L_{0} but larger than 12 for L > L_{0}, Theorem 4.1 for λ ≥ 0 follows.

For λ ≤ 0, we again use Taylor’s formula for small |λ|; arguing as above we see that (4.4) holds for λ ≤ 0 provided

L sup_{−L≤ξ≤0}(−f′′′_{K}(ξ)) ≤ 2K − σ^{2}.  (4.8)

It is easily checked numerically that max_{0≤u≤1} u(1 − u)g(u) < 0.033. It follows that

X(u, λ) ≥ (−η^{3} − 0.396K − 3K^{2}L)e^{ηL}, −L ≤ λ ≤ 0.

Hence, (4.8) holds and (4.4) is satisfied for λ ≤ 0 provided

(η^{3} + 0.396K + 3K^{2}L)Le^{ηL} ≤ 2K − σ^{2}.

It is readily checked that this holds for K = 0.5 and L = 0.62.

For larger L we argue as follows. The function h(u) := g(u) + 4ηu(1 − u) satisfies

h′′(u) = \frac{2}{u(1 − u)} − 8η ≥ 8 − 8η > 0, 0 < u < 1.

Hence h is convex, and since h′(1/2) = 0,

h(u) ≥ h(\frac{1}{2}) = 0, 0 ≤ u ≤ 1.

Consequently, if λ ≤ 0 and K|λ| ≥ 2η, then

λg(U) − 2Kλ^{2}U(1 − U) ≤ λh(U) ≤ 0

and thus f_{K}(λ) ≤ 1. Choosing K = 2η/0.62 < 1.247, this shows f_{K}(λ) ≤ 1 for λ ≤ −0.62, while f_{K}(λ) ≤ f_{0.5}(λ) ≤ 1 for −0.62 ≤ λ ≤ 0 by the preceding case.

This completes the proof of both theorems.

If we just want a bound on ψ_{Y}, (4.2) and Theorems 4.1 and 4.2 can be stated more simply as follows (ignoring the better bounds obtained for small λ).

Corollary 4.3. With L_{0} as in Theorem 4.1,

ψ_{Y} (λ) ≤ \begin{cases}e^{1.25λ^{2}}, & λ ≤ 0, \\ e^{12λ^{2}}, & 0 ≤ λ ≤ L_{0}, \\ e^{2e^{λ}}, & λ ≥ L_{0}. \end{cases}

In particular, ψ_{Y} (λ) ≤ exp(max(12λ^{2}, 2e^{λ})).

The bound e^{2e^{λ}} is very large even for moderately large λ, but the next result shows that ψ_{Y} (λ) really is of essentially this size. In particular, it follows that ln ln ψ_{Y} (λ) ∼ λ as λ → +∞.

Theorem 4.4. If γ < 2/e, then for sufficiently large λ,

ψ_{Y} (λ) ≥ exp(γλ^{−1}e^{λ}).

Proof. Since a moment generating function is convex and ψ′_{Y} (0) = EY = 0, ψ_{Y} is increasing on [0, ∞). Moreover, g is decreasing on [0, 1/2]. Hence, if 0 ≤ δ ≤ 1/2, the integrand in (4.1) with Z = Y is for 0 ≤ u ≤ δ at least ψ_{Y} (0)ψ_{Y} ((1 − δ)λ)e^{λg(δ)} and the same holds for 1 − δ ≤ u ≤ 1 by symmetry. Consequently,

ψ_{Y} (λ) ≥ 2 ∫_{0}^{δ} ψ_{Y} (uλ) ψ_{Y} ((1 − u)λ) e^{λg(u)} du ≥ 2δψ_{Y} ((1 − δ)λ)e^{λg(δ)}, 0 ≤ δ ≤ 1/2.

(4.9)

Let a > 1/2 be a constant to be determined later and choose δ := ae^{−λ}. Then g(δ) = 1 − O(λe^{−λ}) and thus by (4.9), for λ ≥ ln(2a),

ψ_{Y} (λ) ≥ 2ae^{−O(λ^{2}e^{−λ})}ψ_{Y} (λ − aλe^{−λ}).

If 0 < ε < 2a, there thus exists A such that for λ ≥ A,

ψ_{Y} (λ) ≥ (2a − ε)ψ_{Y} (λ − aλe^{−λ}).

Given λ ≥ A, let λ_{0} := λ and define inductively λ_{n+1} := λ_{n} − aλ_{n}e^{−λ_{n}}, n ≥ 0. Let N be the smallest integer with λ_{N} < A. Then ψ_{Y} (λ_{n}) ≥ (2a − ε)ψ_{Y} (λ_{n+1}), n = 0, …, N − 1, and thus

ψ_{Y} (λ) = ψ_{Y} (λ_{0}) ≥ (2a − ε)^{N}ψ_{Y} (λ_{N}) ≥ (2a − ε)^{N}.

It remains to estimate N from below. Since e^{x} is increasing,

∫_{λ_{n+1}}^{λ_{n}} e^{x} dx ≤ e^{λ_{n}}(λ_{n} − λ_{n+1}) = aλ_{n} ≤ aλ

and thus

Naλ ≥ ∫_{λ_{N}}^{λ_{0}} e^{x} dx ≥ ∫_{A}^{λ} e^{x} dx = e^{λ} − e^{A}.

Consequently,

ln ψ_{Y} (λ) ≥ N ln(2a − ε) ≥ \frac{ln(2a − ε)}{a}λ^{−1}(e^{λ} − e^{A}), λ ≥ A.

We choose a = e/2, which maximizes ln(2a)/a. Then ln(2a)/a = 2/e, we may choose ε so small that ln(2a − ε)/a > γ, and the result follows.

As is well known, bounds on the moment generating function yield bounds on the tails of the distribution.

Theorem 4.5. If y ≥ 2e^{L_{0}} = 12L_{0}^{2} ≐ 302.1, then

P(Y ≥ y) ≤ exp(−y(ln y − 1 − ln 2)).

Proof. For λ ≥ L_{0}, by Corollary 4.3,

P(Y ≥ y) ≤ e^{−λy}Ee^{λY} ≤ exp(2e^{λ} − yλ),

and the result follows by taking λ = ln(y/2).

Remark 4.6. The same estimate holds for every Z_{n} provided, say, ψ_{Z_{0}}(λ) ≤ exp(12λ^{2}); for example, when Z_{0} ≡ 0.

Theorem 4.4 suggests that the true size of P(Y ≥ y) is (for large y) not much smaller than the upper bound in Theorem 4.5. Indeed, Knessl and Szpankowski [11] have found (assuming an as yet unverified regularity hypothesis) a much more precise formula for the asymptotics of P(Y ≥ y) which is of the order exp(−y[ln y + ln ln y + O(1)]).

For the left tail, Corollary 4.3 similarly implies P(Y ≤ y) ≤ exp(−y^{2}/5) for y ≤ 0, but this result is much weaker than the doubly exponential decay found by Knessl and Szpankowski [11].

5 Geometric rate of convergence for moment generating functions

Theorem 5.1. Suppose that Z_{0} has mean zero and an everywhere finite moment generating function ψ_{Z_{0}}. Then ψ_{Z_{n}}(λ) → ψ_{Y} (λ) at a geometric rate for every fixed λ ∈ R. Explicitly, if L ≥ 0 and K_{L} are such that (4.2) and (4.3) hold, and if moreover

ψ_{Z_{0}}(λ) ≤ e^{K_{L}λ^{2}},  |λ| ≤ L,  (5.1)

then, for every n ≥ 0 and |λ| ≤ L/2,

\begin{aligned}|ψ_{Z_{n}}(λ) − ψ_{Y} (λ)| & ≤ (Var Z_{0} + σ^{2})^{1/2}|λ|(ψ_{Z_{n}}(2λ) + ψ_{Y} (2λ))^{1/2}(2/3)^{n/2} \\ & ≤ 2^{1/2}(Var Z_{0} + σ^{2})^{1/2}|λ|e^{2K_{L}λ^{2}}(2/3)^{n/2}. \end{aligned}

Of course, if the hypotheses in the first sentence of the theorem’s statement are met, then, given L ≥ 0, (5.1) holds for some K_{L} < ∞, which by Theorem 4.1 may be chosen so large that (4.2) and (4.3) hold, too.

Proof. By (4.3) and induction, the estimate (5.1) holds for every ψ_{Z_{n}}. Fix n ≥ 0 and consider the optimal d_{2}-coupling of (the laws of) Z_{n} and Y. Then for λ ∈ [−L/2, L/2] we have, using the mean value theorem and the Cauchy–Schwarz inequality,

|E e^{λZ_{n}} − E e^{λY} | ≤ E |e^{λZ_{n}} − e^{λY} |
≤ E (|λ||Z_{n} − Y |e^{max(λZ_{n},λY)})
≤ |λ|(E|Z_{n} − Y |^{2})^{1/2}(E e^{2 max(λZ_{n},λY)})^{1/2}
≤ |λ|(E|Z_{n} − Y |^{2})^{1/2}(E e^{2λZ_{n}} + E e^{2λY} )^{1/2}.

By the optimality of the coupling and (1.7),

(E|Z_{n} − Y |^{2})^{1/2} = d_{2}(F_{n}, F) ≤ (Var Z_{0} + σ^{2})^{1/2}(2/3)^{n/2},

and by (4.2) and (5.1) for ψ_{Z_{n}}

E e^{2λZ_{n}} + E e^{2λY} ≤ 2e^{K_{L}(2λ)^{2}},

whence the result follows.

Note further that the operator ψ_{Z} ↦ ψ_{SZ} given by (4.1) is monotone, in the specific sense that if ψ_{Z}(λ) ≤ ψ_{W} (λ) for |λ| ≤ L, then also ψ_{SZ}(λ) ≤ ψ_{SW} (λ) for |λ| ≤ L. In particular, by induction, if ψ_{Z_{0}}(λ) ≤ ψ_{Z_{1}}(λ) for |λ| ≤ L, then ψ_{Z_{n}}(λ) increases monotonically to its limit ψ_{Y} (λ) for |λ| ≤ L. Likewise, if ψ_{Z_{0}}(λ) ≥ ψ_{Z_{1}}(λ) for |λ| ≤ L, then ψ_{Z_{n}}(λ) decreases to ψ_{Y} (λ) for |λ| ≤ L.

We give two simple special cases.

Corollary 5.2. Suppose that Z_{0} ≡ 0. Then ψ_{Z_{n}}(λ) increases monotonically to ψ_{Y} (λ) for every fixed λ. If L ≥ 0 and K_{L} are such that (4.2) and (4.3) hold, then, for every n ≥ 0 and |λ| ≤ L/2,

\begin{aligned}0 ≤ ψ_{Y} (λ) − ψ_{Z_{n}}(λ) & ≤ 2^{1/2}σ|λ| (ψ_{Y} (2λ))^{1/2} (2/3)^{n/2} \\ & ≤ 2^{1/2}σ|λ|e^{2K_{L}λ^{2}}(2/3)^{n/2}. \end{aligned}

Proof. Since E Z_{1} = 0, by Jensen’s inequality ψ_{Z_{1}}(λ) ≥ 1 = ψ_{Z_{0}}(λ), and the monotonicity follows. In particular, ψ_{Z_{n}}(2λ) ≤ ψ_{Y} (2λ), and since (5.1) trivially is satisfied, the result follows from Theorem 5.1.

Corollary 5.3. Suppose L ≥ 0 and K_{L} are such that (4.2) and (4.3) hold, and let Z_{0} ∼ N(0, 2K_{L}). Then ψ_{Z_{n}}(λ) decreases monotonically to ψ_{Y} (λ) for every fixed λ with |λ| ≤ L, and, for every n ≥ 0 and |λ| ≤ L/2,

\begin{aligned}0 ≤ ψ_{Z_{n}}(λ) − ψ_{Y} (λ) & ≤ (4K_{L} + 2σ^{2})^{1/2}|λ| (ψ_{Z_{n}}(2λ))^{1/2} (2/3)^{n/2} \\ & ≤ (4K_{L} + 2σ^{2})^{1/2}|λ|e^{2K_{L}λ^{2}}(2/3)^{n/2}. \end{aligned}

Proof. Since ψ_{Z_{0}}(λ) = e^{K_{L}λ^{2}}, (4.3) yields ψ_{Z_{1}}(λ) ≤ ψ_{Z_{0}}(λ), and the monotonicity follows. The estimate thus follows from Theorem 5.1.

6 On numerical calculations

The preceding results make it possible, in principle at least, to calculate the density, distribution, characteristic, and moment generating functions of Y numerically, with provable arbitrarily high accuracy.

To begin, the results of earlier sections show that it suffices to start with a suitable L(Z_{0}), for example unit mass at 0 or a normal distribution, and then calculate the corresponding quantity for Z_{n}, for a large n that can be determined. The distribution of Z_{n} can be calculated recusively; for the characteristic and moment generating functions we have the recurrence relations (2.4) and (4.1), while for the density functions we have the following recursion:

Theorem 6.1. If n ≥ 0 is arbitrary and Z_{0} has a bounded continuous density function f_{0}, or if Z_{0} is arbitrary and n ≥ 3, then Z_{n} and Z_{n+1} have bounded continuous density functions f_{n} and f_{n+1} satisfying the identity

f_{n+1}(x) = ∫_{u=0}^{1}∫_{z∈R}f_{n}(z) f_{n} (\frac{x − g(u) − (1 − u)z}{u} ) \frac{1}{u} dz du,  x ∈ R,  (6.1)

with g(·) given by (1.3).

Proof. Our proof (similar to that of Theorem 4.1 in [6]) is by induction on n ≥ 0 in the first case and on n ≥ 3, using Theorem 3.1 to get started, in the second case. We may therefore assume as our induction hypothesis that f_{n} is bounded and continuous. It is easily checked that, for each 0 < u < 1, the inner integral

h_{u}(x) := ∫_{z∈R}f_{n}(z) f_{n} (\frac{x − g(u) − (1 − u)z}{u} ) \frac{1}{u} dz

is a density function for the random variable

uZ_{n} + (1 − u)Z_{n}^{∗} + g(u),  (6.2)

and, using dominated convergence, that h_{u} is bounded and continuous. Indeed, h_{u}(x) ≤ (sup f_{n})/u, and since h_{u} = h_{1−u} by symmetry in (6.2), h_{u}(x) ≤ 2 sup f_{n}, uniformly in u and x. It follows, by dominated convergence again, that x ↦ f_{n+1}(x) = ∫_{0}^{1}h_{u}(x) du is a bounded continuous density for Z_{n+1}.

The integrals in (2.4), (4.1), or (6.1) have to be computed numerically—as does the integral of f_{n} to get F_{n}—but that can be done with arbitrary precision since the results above provide bounds for the integrands and their derivatives. [The function g(u) has an unbounded derivative as u → 0 or u → 1, but that can be handled by truncating the interval.] Consequently, to calculate ϕ_{Z_{n}}(t) with given precision for a given t, it suffices to know ϕ_{Z_{n−1}}(t_{k}) with another given precision for a finite number of points t_{k}, which can be done recursively. (However, a brute force recursion along these lines seems to require too many numerical integrations to be practical if we want reasonably high provable accuracy.)

Remark 6.2. To calculate the density f_{n} numerically, it might be better to compute ϕ_{Z_{n}} recursively by (2.4) and then use (3.1), instead of using the recursion (6.1) directly. This is both because (6.1) is a double integral and because we have the simple bounds |ϕ_{n}| ≤ 1 and |ϕ_{n}^{(k)}| ≤ E|Z_{n}|^{k}, k ≥ 1.

7 The metrics d_{p} and a lower bound on d_{2}(F_{n}, F)

At (1.7) we recalled Rösler’s fundamental result

d_{2}(F_{n}, F) = O(ρ^{n})

with ρ := (2/3)^{1/2}. The question naturally arises as to whether there is a lower bound that matches at least to the extent that

d_{2}(F_{n}, F) = Ω(r^{n})

for some r > 0. Of course, the answer is negative without any further restrictions, since if F_{0} = F then F_{n} = F for every n. However, our main result of this section asserts that this is the only exception, at least among distributions F_{0} with finite moments of all orders:

Theorem 7.1. If F_{0} ≠ F has finite moments of all orders, then there exists r > 0 (depending on F_{0}) so that

d_{2}(F_{n}, F) = Ω(r^{n}).

Our arguments for Theorem 7.1 will require use of metrics d_{p} generalizing (1.5). So we will warm up in Section 7.1 by recalling the definition of, and two useful facts about, d_{p} and in Section 7.2 by extending the upper bound result (1.7) to d_{p} for p ≥ 1. Then in Section 7.3 we will prove a sharpened version of Theorem 7.1 (namely, Theorem 7.7).

7.1 The metrics d_{p}

For real 1 ≤ p < ∞, let ∥X∥_{p} := (E|X|^{p})^{1/p} denote the L^{p}-norm, and let d_{p} denote the metric on the space of probability distributions with finite L^{p}-norm defined by

d_{p}(F, G) := min ∥X − Y ∥_{p},

taking the minimum, as at (1.5), over all couplings of L(X) = F and L(Y) = G. It is worth noting that there is a coupling [namely, X = F^{−1}(U) and Y = G^{−1}(U) for U uniform and a suitable definition of the inverse probability transform F^{−1}] that achieves the minimum simultaneously for each 1 ≤ p < ∞ (assuming F and G have finite moments of all orders): see [1].

We begin with two elementary facts that will be useful later. The proof of the first fact (Lemma 7.2) shows that S is a contraction for the d_{p}-metric.

Lemma 7.2. Consider real 1 ≤ p < ∞. The d_{p}-distance from the limiting Quicksort distribution F does not increase when the operator S of (1.4) is applied. Therefore, d_{p}(F_{n}, F) is nonincreasing, and hence bounded, in n if E |Z_{0}|^{p} < ∞.

Proof. With a slight abuse of notation, we find, for Z with any law,

d_{p}(SZ, Y) = d_{p}(SZ, SY) ≤ ∥U(Z − Y) + (1 − U)(Z^{∗} − Y^{∗})∥_{p},  (7.1)

coupling (Y, Z) optimally and (Y^{∗}, Z^{∗}) optimally and choosing U, (Y, Z), and (Y^{∗}, Z^{∗}) to be independent. In calculating the L^{p}-norm value on the right in (7.1), condition on U and then use subadditivity of L^{p}-norm together with independence to bound that value by ∥Z − Y ∥_{p} = d_{p}(Z, Y). This establishes the first assertion: d_{p}(SZ, Y) ≤ d_{p}(Z, Y).

Therefore, d_{p}(F_{n}, F) = d_{p}(S^{n}Z_{0}, Y) is nonincreasing, and hence bounded by d_{p}(Z_{0}, Y), which is bounded by ∥Z_{0}∥_{p} + ∥Y ∥_{p} < ∞ if E |Z_{0}|^{p} < ∞.

Remark 7.3. Conversely, if ∥SZ∥_{p} < ∞, then E|uZ + (1 − u)Z^{∗}|^{p} < ∞ for some u ∈ (0, 1), and thus E|Z|^{p} < ∞ too. Hence, if E|Z_{0}|^{p} = ∞, then E|Z_{n}|^{p} = ∞ and d_{p}(Z_{n}, Y) = ∞ for all n.

Lemma 7.4. For real 2 < p < q < ∞ we have, for any F and G,

d_{p}(F, G) ≤ d2^{\frac{2(q−p)}{p(q−2)}} (F, G) × dq^{\frac{q(p−2)}{p(q−2)}} (F, G).

Proof. Using the common optimal coupling for d_{2} and d_{q}, this is immediate from the inequality

∥X∥_{p} ≤ ∥X∥2^{\frac{2(q−p)}{p(q−2)}} ∥X∥q^{\frac{q(p−2)}{p(q−2)}} ,

which in turn follows from the fact [16, Exercise 4(b) of Chapter 3] that ln ∥X∥_{p}^{p} is convex in p ∈ (0, ∞).

7.2 Geometric rate of convergence in each metric d_{p}

Under suitable conditions, we can establish a geometric rate of convergence for d_{p}(F_{n}, F) for any real 1 ≤ p < ∞. We begin with an elementary lemma.

Lemma 7.5. If p ≥ 2, then for all x, y ≥ 0,

(x + y)^{p} ≤ x^{p} + y^{p} + c_{p}(x^{p−1}y + xy^{p−1}),

where c_{p} := p(p − 1)2^{p−2}.

Proof.

\begin{aligned}(x + y)^{p} − x^{p} − y^{p} & = ∫_{0}^{y} p((x + t)^{p−1} − t^{p−1}) dt = ∫_{t=0}^{y} ∫_{u=0}^{x} p(p − 1)(t + u)^{p−2} du dt \\ & ≤ p(p − 1)xy(x + y)^{p−2} ≤ p(p − 1)xy2^{p−2}(x^{p−2} + y^{p−2}). \end{aligned}

Theorem 7.6. Let p_{0} ≐ 6.557 be the largest positive solution to

(\frac{2}{p_{0} + 1})^{1/p_{0}} = (\frac{2}{3})^{1/2}  (7.2)

and let, for any ε > 0,

β_{p} := \begin{cases}(\frac{2}{3})^{1/2}, & 1 ≤ p < p_{0}, \\ (\frac{2}{3})^{1/2} + ε, & p = p_{0}, \\ (\frac{2}{p+1})1/p & p > p_{0}. \end{cases}  (7.3)

Thus, for p ≥ 2 except p = p_{0}, β_{p} = max((2/3)^{1/2}, (2/(p + 1))^{1/p}). Then, for any Z_{0} with zero mean and finite variance, and every p ≥ 1 such that E|Z_{0}|^{p} < ∞, there exists a constant α_{p} < ∞ [depending on L(Z_{0})] such that

d_{p}(Z_{n}, Y) ≤ α_{p}β_{p}^{n}.  (7.4)

Proof. First we note that (7.2) can be written (3/2)^{p_{0}/2} = (p_{0} + 1)/2. One root of this equation is 2, and since (3/2)^{p/2} is convex, with derivative less than 1/2 at p = 2, it follows that the equation has two positive roots, 2 and p_{0} > 2, and that (2/3)^{p/2} > 2/(p + 1) for 2 < p < p_{0}, while (2/3)^{p/2} < 2/(p + 1) for p > p_{0}.

Next we note that (7.4) holds for p ≤ 2, with

α_{p} := (Var Z_{0} + σ^{2})^{1/2}, p ≤ 2,

by (1.7) and the inequality d_{p} ≤ d_{2}, p ≤ 2. We then proceed by induction on ⌊p⌋. For the induction step, suppose that p > 2 and that Z_{0} has zero mean and satisfies ∥Z_{0}∥_{p} < ∞. By the induction hypothesis, there exist constants 0 < α_{q} < ∞, 1 ≤ q ≤ p − 1, such that

d_{q}(Z_{n}, Y) ≤ α_{q}β_{q}^{n} for all q ≤ p − 1 and n ≥ 0.  (7.5)

Using our usual coupling of Z_{n} and Y in terms of the optimal coupling of Z_{n−1} and Y, we find easily by Lemma 7.5, for n ≥ 1 [with (Z_{n−1}, Y), (Z_{n−1}^{∗}, Y^{∗}), and U independent],

d_{p}^{p}(Z_{n}, Y) ≤ E|U(Z_{n−1} − Y) + (1 − U)(Z_{n−1}^{∗} − Y^{∗})|^{p}

≤ E(U|Z_{n−1} − Y | + (1 − U)|Z_{n−1}^{∗} − Y^{∗}|)^{p}

≤ E(U^{p}|Z_{n−1} − Y |^{p}) + E((1 − U)^{p}|Z_{n−1}^{∗} − Y^{∗}|^{p})

+ c_{p}E(U^{p−1}(1 − U)|Z_{n−1} − Y |^{p−1}|Z_{n−1}^{∗} − Y^{∗}|)
+ c_{p}E(U(1 − U)^{p−1}|Z_{n−1} − Y ||Z_{n−1}^{∗} − Y^{∗}|^{p−1})

= \frac{2}{p + 1}d_{p}^{p}(Z_{n−1}, Y) + \frac{2c_{p}}{p(p + 1)}d_{1}(Z_{n−1}, Y) d_{p−1}^{p−1}(Z_{n−1}, Y).

So by induction on n it follows that

d_{p}^{p}(Z_{n}, Y) ≤ (\frac{2}{p + 1})^{n}d_{p}^{p}(Z_{0}, Y) + \frac{c_{p}}{p} ∑_{i=0}^{n−1}(\frac{2}{p + 1})^{n−i}d_{1}(Z_{i}, Y)d_{p−1}^{p−1}(Z_{i}, Y).

By the induction hypothesis (7.5) this yields, for some a_{1}, a_{2} < ∞ (depending on p),

d_{p}^{p}(Z_{n}, Y) ≤ a_{1}(\frac{2}{p + 1})^{n} + a_{2} ∑_{i=0}^{n−1}(\frac{2}{p + 1})^{n−i}(β_{1}β_{p−1}^{p−1})^{i}.  (7.6)

Let γ := β_{1}β_{p−1}^{p−1}. We break our treatment into three cases:

(i) If γ > 2/(p + 1), we write the sum in (7.6) as

γ^{n} ∑_{i=0}^{n−1} (\frac{2}{p + 1}γ^{−1})^{n−i} < (1 − \frac{2}{p + 1}γ^{−1})^{−1}γ^{n},

and thus (7.6) shows that (7.4) holds with β_{p}^{p} = γ.

(ii) If γ < 2/(p + 1), we write the sum in (7.6) as

(\frac{2}{p + 1})^{n} ∑_{i=0}^{n−1} (γ \frac{p + 1}{2})^{i} < (1 − γ \frac{p + 1}{2})^{−1}(\frac{2}{p + 1})^{n},

and thus (7.4) holds with β_{p}^{p} = 2/(p + 1).

(iii) If γ = 2/(p + 1), the sum in (7.6) equals n(2/(p + 1))^{n}. Consequently, (7.4) holds with any β_{p} > (2/(p + 1))^{1/p}.

It remains to verify that this yields the β_{p} given in (7.3).

First, if 2 < p < p_{0}, then the induction hypothesis yields γ = β_{1}β_{p−1}^{p−1} = (2/3)^{p/2} > 2/(p + 1), so case (i) gives β_{p}^{p} = (2/3)^{p/2}. Similarly, for p = p_{0}, γ = (2/3)^{p/2} = 2/(p + 1) and (iii) shows that any β_{p} > (2/(p + 1))^{1/p} = (2/3)^{1/2} will do.

For p_{0} < p < p_{0} + 1, we have γ = β_{1}β_{p−1}^{p−1} = (2/3)^{p/2} < 2/(p + 1), so case (ii) yields β_{p}^{p} = 2/(p + 1). The same applies for p = p_{0} + 1, since again (2/3)^{p/2} < 2/(p + 1) and we thus may choose ε so small that γ = (\frac{2}{3})^{1/2}((\frac{2}{3})^{1/2} + ε)^{p−1} < 2/(p + 1).

Finally, for p > p_{0} + 1, we have

γ = β_{1}β_{p−1}^{p−1} = (\frac{2}{3})^{1/2} \frac{2}{p} < \frac{2}{p + 1},

since p/(p+1) is increasing and equals (2/3)^{1/2} when p = (\sqrt{3/2}−1)^{−1} = 2(\sqrt{6}−2)^{−1} = \sqrt{6} + 2 < 5 < p_{0}. Hence case (ii) applies.

7.3 Lower bounds

The main goal of this subsection is to establish Theorem 7.1, or rather the sharper Theorem 7.7 below. Since (as noted in Section 4) the limiting Quicksort distribution F has everywhere finite moment generating function, it is uniquely determined by its moments. Hence if F_{0} ≠ F has finite moments of all orders, then E Z_{0}^{j} ≠ E Y^{j} for some integer j ≥ 1.

Theorem 7.7. Suppose F_{0} ≠ F has finite moments of all orders, and let p be the smallest positive integer such that E Z_{0}^{p} ≠ E Y^{p}. Then, for any 0 < r < (\frac{2}{p+1})^{p/2},

d_{2}(F_{n}, F) = Ω(r^{n}).

(The implicit multiplicative constant depends on both F_{0} and r.)

Remark 7.8. The cases p = 1 and p = 2 are a bit special, and in these cases we claim that Theorem 7.7 holds even with r = (\frac{2}{p+1})^{p/2}, i.e., with r = 1 and r = 2/3, respectively, and without the assumption that F_{0} has finite moments.

We may and shall assume that E Z_{0}^{2} < ∞, since otherwise d_{2}(F_{n}, F) = ∞.

First, p = 1 when E Z_{0} ≠ E Y = 0; in this case Z_{n} converges in distribution to Y + E Z_{0} and not to Y, and thus inf d_{2}(Z_{n}, Y) > 0, i.e., the theorem holds with r = 1. Indeed, we have the sharper result that

d_{2}(Z_{n}, Y) = d_{2}(Y + E Z_{0}, Y) + O(d_{2}(Z_{n}, Y + E Z_{0})) = |E Z_{0}| + O((2/3)^{n/2}).

Next, p = 2 when E Z_{0} = 0 but Var Y ≠ Var Z_{0}; in this case Theorem 7.9 shows that the result holds with r = 2/3. Even in this case we have a gap between the lower bound Ω((2/3)^{n}) and Rösler’s upper bound O((2/3)^{n/2}); it is an open problem to find the rate of approximation more precisely.

We prove Theorem 7.7 using the following Theorem 7.9, which is a similar lower bound for the d_{p}-metric.

Theorem 7.9. Let p ≥ 1 be an integer, and suppose that E Z_{0}^{j} = E Y^{j} for integers 1 ≤ j ≤ p − 1, and that E Z_{0}^{p} exists and is finite but fails to equal E Y^{p}. Then

d_{p}(F_{n}, F) = Ω ((\frac{2}{p+1})^{n}) .

Theorem 7.9 is, in turn, a simple consequence of the following two elementary lemmas.

The first lemma demonstrates a sense in which the value of p in Theorems 7.7 and 7.9 persists from F_{0} to each F_{n}; the second gives a general lower bound on d_{p} in terms of discrepancy in pth moments.

Lemma 7.10. Let p ≥ 1 be an integer, and suppose for n = 0 that E Z_{n}^{j} = E Y^{j} for integers 1 ≤ j ≤ p − 1, and that E Z_{n}^{p} exists and is finite but fails to equal E Y^{p}. Then for every n ≥ 0 the same is true and, moreover,

E Z_{n}^{p} − E Y^{p} = (\frac{2}{p+1})^{n} (E Z_{0}^{p} − EY^{p}) .

Proof. If E|Z|^{m} < ∞, then, with Z, Z^{∗} =^{L} Z, and U independent, by (1.4) and a trinomial expansion we have

\begin{aligned}E(SZ)^{m} & = ∑_{j+k≤m} \frac{m!}{j! k! (m − j − k)!}E(U^{j}Z^{j}(1 − U)^{k}(Z^{∗})^{k}g(U)^{m−j−k}) \\ & = ∑_{j+k≤m} \frac{m!}{j! k! (m − j − k)!}E(U^{j}(1 − U)^{k}g(U)^{m−j−k}) EZ^{j} EZ^{k}. \end{aligned}

We apply this with m = 1, …, p for both Z = Z_{n−1} and Z = Y, and note that by induction on n, all terms in the sum with j ≤ p − 1 and k ≤ p − 1 coincide for the two choices of Z. Hence, EZ_{n}^{m} = EY^{m} for 1 ≤ m ≤ p − 1 and

EZ_{n}^{p} − EY^{p} = (EU^{p} + E(1 − U)^{p})(EZ_{n−1}^{p} − EY^{p}) = \frac{2}{p + 1}(EZ_{n−1}^{p} − EY^{p}),

and the result follows.

Lemma 7.11. Let p ≥ 1 be an integer. Then, for any F and G,

d_{p}(F, G) ≥ \frac{|E X^{p} − E Y^{p}|}{∑_{j=0}^{p−1} ∥X∥_{p}^{p−1−j} ∥Y ∥_{p}^{j}}

with X ∼ F and Y ∼ G (and 0^{0} := 1).

Proof. Let (X, Y) be an optimal coupling of F and G. If p = 1, then

d_{1}(F, G) = ∥X − Y ∥_{1} = E |X − Y | ≥ |E X − E Y |,

as desired. If p ≥ 2, we employ the factorization

X^{p} − Y^{p} = (X − Y) ∑_{j=0}^{p−1} X^{p−1−j}Y^{j},

whence

|E X^{p} − E Y^{p}| ≤ E|X^{p} − Y^{p}|
≤ ∥X − Y ∥_{p} ∥∑_{j=0}^{p−1} X^{p−1−j}Y^{j}∥_{p/(p−1)}

≤ d_{p}(F, G) ∑_{j=0}^{p−1} ∥X^{p−1−j}Y^{j}∥_{p/(p−1)} ,  (7.7)

where at the second inequality we have employed Hölder’s inequality and at the third we have invoked the optimality of the coupling. Another application of Hölder’s inequality, this time with conjugate exponents \frac{p−1}{p−1−j} and \frac{p−1}{j}, yields

∥X^{p−1−j}Y^{j}∥_{p/(p−1)} ≤ ∥X∥_{p}^{p−1−j} ∥Y ∥_{p}^{j}  (7.8)

for 1 ≤ j ≤ p − 2, and (7.8) is trivially an equality when j = 0 or j = p − 1. Combining (7.7) and (7.8) and rearranging, we obtain the desired result.

Proof of Theorem 7.9. By Lemma 7.2, we have the bound

∥Z_{n}∥_{p} ≤ d_{p}(Z_{n}, Y) + ∥Y ∥_{p} ≤ d_{p}(Z_{0}, Y) + ∥Y ∥_{p} ≤ ∥Z_{0}∥_{p} + 2∥Y ∥_{p}.

Thus Lemmas 7.11 and 7.10 yield the explicit bound

d_{p}(Z_{n}, Y) ≥ \frac{|E Z_{0}^{p} − E Y^{p}|}{∑_{j=0}^{p−1} (∥Z_{0}∥_{p} + 2∥Y ∥_{p})^{p−1−j} ∥Y ∥_{p}^{j}} (\frac{2}{p + 1})^{n}.

Proof of Theorem 7.7. The cases p = 1 and p = 2 follow immediately from Theorem 7.9; see also Remark 7.8.

When p ≥ 3, fix q > p. By Lemmas 7.4 and 7.2 (the latter applied to d_{q}), for some C_{q} we have

d_{p}(F_{n}, F) ≤ C_{q} d2^{\frac{2(q−p)}{p(q−2)}} (F_{n}, F)

and thus Theorem 7.9 implies Theorem 7.7 with r = (\frac{2}{p+1})^{\frac{p(q−2)}{2(q−p)}} . By taking q sufficiently

large, we obtain the result for any r < (\frac{2}{p+1})^{p/2}.

We have assumed in Theorems 7.1 and 7.7 that F_{0} has finite moments of all orders. What happens if this fails? If E |Z_{0}|^{p} = ∞ for some p > 2, then d_{p}(Z_{n}, Y) = ∞ for all n, but what can be said about d_{2}(Z_{n}, Y)? It seems reasonable to conjecture that we have at least as large d_{2}-distance in this case as in the nicer case with all moments finite, and that Theorems 7.1 and 7.7 hold for all F_{0} ≠ F. Unfortunately, we have not been able to prove this, but we offer the following partial result.

Theorem 7.12. If d_{2}(F_{n}, F) = O(r^{n}) for every r > 0, then F_{0} = F. Consequently, if F_{0} ≠ F, there exists r > 0 such that d_{2}(F_{n}, F) > r^{n} for infinitely many values of n.

Proof. We may assume that EZ_{0} = 0 and EZ_{0}^{2} = EY^{2} (in particular, EZ_{0}^{2} < ∞), because otherwise d_{2}(F_{n}, F) = Ω(r^{n}) with r = 2/3: see Remark 7.8. By induction then EZ_{n} = 0 and EZ_{n}^{2} = EY^{2} < 1 for every n: see Lemma 7.10.

As usual, let Z, Z^{∗} =^{L} Z, and U be independent. If |Z| ≥ 2x, |Z^{∗}| ≤ 2, and \frac{2}{3} ≤ U ≤ 1, where x ≥ 5, then

|SZ| = |UZ + (1 − U)Z^{∗} + g(U)| ≥ \frac{2}{3}|Z| − \frac{1}{3}|Z^{∗}| − 1 ≥ \frac{4}{3}x − \frac{5}{3} ≥ x.

Thus,

P(|SZ| ≥ x) ≥ P(|Z| ≥ 2x) · P(|Z| ≤ 2) · \frac{1}{3}, x ≥ 5.

If further EZ^{2} ≤ 1, and thus by Chebyshev’s inequality P(|Z| ≤ 2) = 1 − P(|Z| > 2) ≥ 1 − \frac{1}{4} = \frac{3}{4}, this yields

P(|SZ| ≥ x) ≥ \frac{1}{4}P(|Z| ≥ 2x), x ≥ 5.

Hence, by induction on n and our assumption on the first two moments of Z_{0}, for any x ≥ 5,

P(|Z_{n}| ≥ x) ≥ 4^{−n}P(|Z_{0}| ≥ 2^{n}x), n ≥ 0,

and in particular

P(|Z_{n}| ≥ 2^{n}) ≥ 4^{−n}P(|Z_{0}| ≥ 4^{n}),  n ≥ 3.  (7.9)

Now suppose that d_{2}(Z_{n}, Y) = O(r^{n}). Using an optimal coupling between Z_{n} and Y, and the fact that Y has moments of all orders, we find

\begin{aligned}P(|Z_{n}| ≥ 2^{n}) & ≤ P(|Z_{n} − Y | ≥ 2^{n−1}) + P(|Y | ≥ 2^{n−1}) \\ & ≤ 2^{2−2n}d_{2}^{2}(Z_{n}, Y) + P(|Y | ≥ 2^{n−1}) = O(2^{−2n}r^{2n}). \end{aligned}

Combining this with (7.9), we obtain, for n ≥ 3,

P(|Z_{0}| ≥ 4^{n}) ≤ 4^{n}P(|Z_{n}| ≥ 2^{n}) = O(r^{2n}),

which implies that E|Z_{0}|^{p} < ∞ for every p > 0 such that 4^{p}r^{2} < 1.

Consequently, if d_{2}(Z_{n}, Y) = O(r^{n}) for every r > 0, then E|Z_{0}|^{p} < ∞ for every p > 0, and Theorem 7.7 applies to yield F_{0} = F.

Remark 7.13. Our proof of Theorem 7.12, combined with the proof of Theorem 7.7, shows that if F_{0} ≠ F and p (assumed ≥ 3 here) is the smallest positive integer such that either E Z_{0}^{p} does not exist or E Z_{0}^{p} ≠ E Y^{p}, then d_{2}(Z_{n}, Y) > r^{n} for infinitely many values of n for any 0 < r < r_{p}, with r_{p} := 2^{−q}, where q is the unique solution in (p, ∞)

to 2^{\frac{2q(q−p)}{p(q−2)}} = \frac{p+1}{2}.

8 Other lower bounds

In Section 7 we showed that convergence of the iterates F_{n} to F in the d_{2}-metric is not faster than geometric. In this final section we show likewise that the convergence is not faster than geometric in the other metrics we have considered in this paper. We again assume that F_{0} ≠ F has finite moments of all orders. (Without this hypothesis, we can prove partial results by the method used in the proof of Theorem 7.12, but we do not know whether the full results hold.)

8.1 Kolmogorov–Smirnov and total variation distances

Recall the definitions of Kolmogorov–Smirnov distance and total variation distance given at (3.11) and (3.12), respectively. We begin with a simple lemma.

Lemma 8.1. Let p > 0. For any X ∼ F and Y ∼ G each with finite pth absolute moment, if K := d_{KS}(F, G), then, for any 0 ≤ M < ∞,

|E (X^{p}; X > 0) − E (Y^{p}; Y > 0)| ≤ KM^{p} + E (X^{p}; X > M) + E (Y^{p}; Y > M)

and, if p is an integer,

|E X^{p} − E Y^{p}| ≤ 2KM^{p} + E (|X|^{p}; |X| > M) + E (|Y |^{p}; |Y | > M) .

Proof. Define X_{M} := min(X^{+}, M), where X^{+} = max(X, 0), and similarly Y_{M}. Then

0 ≤ E(X^{p}; X > 0) − EX_{M}^{p} ≤ E(X^{p}; X > M)

and similarly for Y, while

\begin{aligned}|E X_{M}^{p} − E Y_{M}^{p} | & = |∫_{0}^{M} px^{p−1}P(X > x) dx − ∫_{0}^{M} px^{p−1}P(Y > x) dx| \\ & ≤ ∫_{0}^{M} px^{p−1}|P(X > x) − P(Y > x)| dx ≤ KM^{p}. \end{aligned}

Together, these yield the first inequality.

The second follows by applying the first to (X, Y) and to (−X, −Y) and summing or subtracting, depending on the parity of p.

Theorem 8.2. Suppose F_{0} ≠ F has finite moments of all orders, and let p be defined as in Theorem 7.7. Then, for any 0 < r < 2/(p + 1),

d_{TV}(F_{n}, F) ≥ d_{KS}(F_{n}, F) = Ω(r^{n}).

(The implicit multiplicative constant depends on both F_{0} and the choice of r.)

Proof. Let K_{n} := d_{KS}(F_{n}, F) > 0. If we apply Lemma 8.1 and then use Lemma 7.10, we find, for any q ≥ p,

\begin{aligned}|E Z_{0}^{p} − E Y^{p}| (\frac{2}{p+1})^{n} & ≤ 2K_{n}M^{p} + E (|Z_{n}|^{p}; |Z_{n}| > M) + E (|Y |^{p}; |Y | > M) \\ & ≤ 2K_{n}M^{p} + M^{−(q−p)}E|Z_{n}|^{q} + M^{−(q−p)}E|Y |^{q}, \end{aligned} (8.1)

for any 0 ≤ M < ∞. It follows from Lemma 7.2 that E|Z_{n}|^{q} ≤ C_{q}, for some C_{q} not depending on n. Choosing M = K_{n}^{−1/q} thus gives, with c := |E Z_{0}^{p} − E Y^{p}| > 0,

c(\frac{2}{p+1})^{n} ≤ 2K_{n}^{1−(p/q)} + 2C_{q}K_{n}^{1−(p/q)},

and thus K_{n} = Ω(r^{n}) with r = (\frac{2}{p+1})^{q/(q−p)}. The result follows, since r → 2/(p + 1) as q → ∞.

8.2 Density functions and characteristic functions

We immediately obtain results for the density functions f_{n}, which by Theorem 3.1 exist at least for n ≥ 3 and by Theorem 3.2 converge uniformly, at a geometric rate, to the density f of Y.

Corollary 8.3. Suppose F_{0} ≠ F has finite moments of all orders, and let p be defined as in Theorem 7.7. Then, for any 0 < r < 2/(p + 1),

∫_{−∞}^{∞} |f_{n}(x) − f(x)| dx = Ω(r^{n})  (8.2)

and

sup_{x} |f_{n}(x) − f(x)| = Ω(r^{n}).  (8.3)

Proof. The estimate (8.2) follows from Theorem 8.2, because (whenever f_{n} exists) ∫_{−∞}^{∞} |f_{n}(x) − f(x)| dx = 2d_{TV}(F_{n}, F).

The estimate (8.3) follows from Theorem 8.2 using inequality (3.14) and the discussion following it.

Similarly, we have a geometric lower bound for the L^{1}(R) and L^{∞}(R) distances of the characteristic functions.

Corollary 8.4. Suppose F_{0} ≠ F has finite moments of all orders, and let p be defined as in Theorem 7.7. Then, for any 0 < r < 2/(p + 1),

∫_{−∞}^{∞} |ϕ_{Z_{n}}(t) − ϕ_{Y} (t)| dt = Ω(r^{n}).  (8.4)

and

sup_{t} |ϕ_{Z_{n}}(t) − ϕ_{Y} (t)| = Ω(r^{n}).  (8.5)

Proof. The estimate (8.4) is immediate from Corollary 8.3 and inequality (3.8). Next, (8.4) for some r = r_{0} and Theorem 2.1 imply (8.5) for any r < r_{0} by the argument used to show (3.10) in the proof of Theorem 3.2.

It is not too hard to extend Corollaries 8.3 and 8.4 to any L^{q}(R) distance, 1 ≤ q ≤ ∞.

8.3 Moment generating functions

Finally, we consider lower bounds for the convergence of moment generating functions. We assume for simplicity that Z_{0} has an everywhere finite moment generating function, and know by Theorem 5.1 that then ψ_{Z_{n}} converges to ψ_{Y} pointwise, and uniformly on compact sets, with geometric rate. For lower bounds we first note that if F_{0} ≠ F and p is as in Theorem 7.7, then the derivatives of ψ_{Z_{n}} and ψ_{Y} at the origin, which equal the corresponding moments of Z_{n} and Y, by Lemma 7.10 coincide up to order p − 1, while the pth derivatives differ by c(\frac{2}{p+1})^{n} with c = EZ_{0}^{p} − EY^{p} ≠ 0. For λ close to the origin, this and a Taylor expansion shows that |ψ_{Z_{n}}(λ) − ψ_{Y} (λ)| = Ω(|λ|^{p}(\frac{2}{p+1})^{n}), but the range of λ where we can prove that this is valid depends on n.

Indeed, there is no general lower bound for |ψ_{Z_{n}}(λ) − ψ_{Y} (λ)| for a fixed λ, since there may be points λ where ψ_{Z_{n}}(λ) and ψ_{Y} (λ) coincide “accidentally”. For example, suppose that Z_{0} is bounded with E Z_{0} = 0 and Var Z_{0} > Var Y. By induction, the same holds for each Z_{n}: see Lemma 7.10 and note that g(U) is bounded. Consequently, for each n, Taylor’s formula shows that ψ_{Z_{n}}(λ) > ψ_{Y} (λ) for small positive λ, while ψ_{Z_{n}}(λ) = exp(O(λ)) and thus Theorem 4.4 shows that ψ_{Z_{n}}(λ) < ψ_{Y} (λ) for large λ. Hence there exists for every n at least one positive λ = λ_{n} such that ψ_{Z_{n}}(λ) = ψ_{Y} (λ). Nevertheless, such points have to be isolated, and if we consider the maximum deviation over an interval, we have a geometric lower bound.

Theorem 8.5. Suppose F_{0} ≠ F has everywhere finite moment generating function, and let (a, b) be a nonempty interval. Then there exists r > 0 such that sup_{a≤λ≤b} |ψ_{Z_{n}}(λ) − ψ_{Y} (λ)| = Ω(r^{n}).

Proof. We use the fact that the moment generating functions ψ_{Z_{n}} and ψ_{Y} are entire analytic functions in the complex plane C.

Let R := |a| + |b| + 1. There exists a (unique) function ω which is continuous on D_{R} := {z ∈ C : |z| ≤ R} and analytic in Ω := {z : |z| < R} \ [a, b] such that ω(z) = 0 for |z| = R and ω(z) = 1 for z ∈ [a, b]; this function is called harmonic measure and is probabilistically given by the probability that a Brownian motion starting at z hits [a, b] before it hits {z : |z| = R}.

Let f_{n}(z) := ψ_{Z_{n}}(z) − ψ_{Y} (z) and u_{n}(z) := ln |f_{n}(z)| ≥ −∞. For z ∈ D_{R},

|f_{n}(z)| ≤ |ψ_{Z_{n}}(z)| + |ψ_{Y} (z)| ≤ ψ_{Z_{n}}(R) + ψ_{Y} (R) + ψ_{Z_{n}}(−R) + ψ_{Y} (−R),

which by Theorem 5.1 is bounded by some constant A < ∞ (depending on Z_{0} but not on n). Let further δ_{n} := max_{a≤λ≤b} |f_{n}(λ)|; we may of course restrict attention to those values of n satisfying δ_{n} < 1. Now u_{n}(z) ≤ ln A for |z| = R and u_{n}(z) ≤ ln δ_{n} for z ∈ [a, b]; thus (since A ≥ 1)

u_{n}(z) ≤ ln A + (ln δ_{n})ω(z)  (8.6)

for every z ∈ ∂Ω. Since u_{n} is subharmonic and the right hand side is harmonic in Ω and continuous on its closure, (8.6) holds for every z ∈ \overline{Ω} = D_{R}, cf. [16, Theorems 17.3 and 17.4]. In particular, setting ε := inf_{|z|≤1} ω(z) > 0, we have

u_{n}(z) ≤ ln A + ε ln δ_{n}, |z| ≤ 1,

or

|f_{n}(z)| ≤ Aδ_{n}^{ε},  |z| ≤ 1.  (8.7)

Let p be as in Theorem 7.7. By (8.7) and Cauchy’s estimates [16, Theorem 10.26],

|f_{n}^{(p)}(0)| ≤ p! Aδ_{n}^{ε}.

Since by Lemma 7.10

|f_{n}^{(p)}(0)| = |EZ_{n}^{p} − EY^{p}| = Ω((\frac{2}{p+1})^{n}),

it follows that δ_{n} = Ω(r^{n}) with r = (\frac{2}{p+1})^{1/ε}.

References

[1] Cambanis, S., Simons, G., and Stout, W. Inequalities for Ek(X, Y) when the marginals are fixed. Z. Wahrscheinlichkeitstheorie verw. Gebiete 36 (1976), 285– 294.

[2] van der Corput, J. G. Zahlentheoretische Abschätzungen. Math. Ann. 84 (1921), 53–79.

[3] Dongarra, J. and Sullivan, F. Guest editors’ introduction: the top 10 algorithms. Computing in Science & Engineering 2 (2000).

[4] Eddy, W. F. and Schervish, M. J. How many comparisons does Quicksort use? J. Algorithms 19 (1995), 402–431.

[5] Feller, W. An Introduction to Probability Theory and its Applications. Vol. II. Second edition. Wiley, New York, 1971.

[6] Fill, J. A. and Janson, S. Smoothness and decay properties of the limiting Quicksort density function. In D. Gardy and A. Mokkadem, editors, Mathematics and Computer Science: Algorithms, Trees, Combinatorics and Probabilities, Trends in Mathematics, pages 53–64. Birkhäuser Verlag, 2000. Refereed article, available from

http://www.mts.jhu.edu/~fill/ or http://www.math.uu.se/~svante/.

[7] Fill, J. A. and Janson, S. A characterization of the set of fixed points of the Quicksort transformation. Electron. Comm. Probab. 5 (2000), 77–84 (electronic).

[8] Fill, J. A. and Janson, S. Quicksort asymptotics. In preparation.

[9] Hoare, C. A. R. Quicksort. Comput. J. 5 (1962), 10–15.

[10] JaJa, J. A perspective on quicksort. Computing in Science & Engineering 2 (2000).

[11] Knessl, C. and Szpankowski, W. Quicksort algorithm again revisited. Discrete Math. Theor. Comput. Sci. 3 (1999), 43–64.

[12] Knuth, D. E. The Art of Computer Programming. Volume 3. Sorting and searching. Second edition. Addison–Wesley, Reading, Mass., 1998.

[13] Montgomery, H. L. Ten Lectures on the Interface Between Analytic Number Theory and Harmonic Analysis. CBMS Reg. Conf. Ser. Math. 84, AMS, Providence, R.I., 1994.

[14] Régnier, M. A limiting distribution for quicksort. RAIRO Inform. Théor. Appl. 23 (1989), 335–343.

[15] Rösler, U. A limit theorem for ‘Quicksort’. RAIRO Inform. Théor. Appl. 25 (1991), 85–100.

[16] Rudin, W. Real and Complex Analysis. Second edition. McGraw–Hill, New York, [17] Schwartz, L. Théorie des Distributions. Second edition. Hermann, Paris, 1966.

[18] Tan, K. H. and Hadjicostas, P. Some properties of a limiting distribution in Quicksort. Statist. Probab. Lett. 25 (1995), 87–94.