arXiv:1704.06398v2 [] 24 Jul 2017

Tail sums of Wishart and GUE eigenvalues beyond the
bulk edge.
Iain M. Johnstone∗
Stanford University and Australian National University
July 26, 2017
Abstract
Consider the classical Gaussian unitary ensemble of size N and the real Wishart
ensemble WN (n, I). In the limits as N → ∞ and N/n → γ > 0, the expected number of eigenvalues that exit the upper bulk edge is less than one, 0.031 and 0.170
respectively, the latter number being independent of γ. These statements are consequences of quantitative bounds on tail sums of eigenvalues outside the bulk which
are established here for applications in high dimensional covariance matrix estimation.

1

Introduction

This paper develops some tail sum bounds on eigenvalues outside the bulk that are needed
for results on estimation of covariance matrices in the spiked model, Donoho et al. [2017].
This application is described briefly in Section 4. It depends on properties of the eigenvalues
of real white Wishart matrices, distributed as WN (n, I), which are the main focus of this
note.
Specifically, suppose that A ∼ WN (n, I), and that λ1 ≥ . . . ≥ λN are eigenvalues
of the sample covariance matrix n−1 A. In the limit N/n → γ > 0, it is well known
that the empirical distribution of {λi } converges to the Marcenko-Pastur law (see e.g.
Pastur and Shcherbina [2011, Corollary 7.2.5]), which is supported on an interval Iγ —
√
augmented with 0 if γ > 1 — having upper endpoint λ(γ) = (1 + γ)2 . We focus on the
eigenvalues λi that exit this “bulk” interval Iγ on the upper side. In statistical application,
such exiting eigenvalues might be mistaken for “signal” and so it is useful to have some
bounds on what can happen under the null hypothesis of no signal. Section 3 studies the
mean value behavior of quantities such as
N
X
TN =
[λi − λ(γ)]q+ ,
i=1

q≥0

Thanks to David Donoho and Matan Gavish for discussions and help. Peter Forrester provided an
important reference. Work on this paper was completed during a visit to the Research School of Finance,
Actuarial Studies and Statistics, A. N. U., whose hospitality and support is gratefully acknowledged. Work
supported in part by the U.S. National Science Foundation and the National Institutes of Health.
∗

1

which for q = 0 reduces to the number TN0 of exiting eigenvalues.
a.s.
It is well known that the largest eigenvalue λ1 → λ(γ) [Geman, 1980], and that closed
intervals outside the bulk support contain no eigenvalues for N large with probability
one [Bai and Silverstein, 1998]. However these and even large deviation results for λ1
[Majumdar and Vergassola, 2009] and TN0 [Majumdar and Vivo, 2012] seem not to directly
yield the information on E(TN ) that we need. Marino et al. [2014] looked at the variance
of TN0 using methods related to those of this note. Recently, Chiani [2017] has studied the
probability that all eigenvalues of Gaussian, Wishart and double Wishart random matrices
lie within the bulk, and derived universal limiting values of 0.6921 and 0.9397 in the real
and complex cases respectively.
In summary, the motivation for this note is high-dimensional covariance estimation, but
there are noteworthy byproducts: the asymptotic values of E(TN0 ) are perhaps suprisingly
small, and numerically for the Gaussian Unitary Ensemble (GUE), it is found that the
chance of even two exiting eigenvalues is very small, of order 10−6 .

2

The Gaussian Unitary Ensemble Case (GUE)

We begin with GUE to illustrate the methods in the simplest setting, and to note an error
in the literature. Recall that the Gaussian Unitary ensemble GUE(N) is the Gaussian
probability measure on the space of N × N Hermitian matrices with density proportional
to exp{− 21 NtrA2 }.
Theorem 1. Let λ1 , . . . , λN be eigenvalues of an N-by-N matrix from the GUE. Denote
by λ+ = 2 the upper edge of the Wigner semicircle, namely, the asymptotic density of the
eigenvalues. For q ≥ 0, let
TN =

N
X
i=1

(λi − λ+ )q+ .

(1)

Then, with a constant cq specified at (3) below,
E(TN ) = cq N −2q/3 (1 + o(1)).
In particular, for q = 0 and TN = #{i : λi > λ+ },
1
E(TN ) → c0 = √ ≈ 0.030629.
6 3π

(2)

Proof. We use the so-called one-point function and bounds due to Tracy and Widom [1994,
1996]. To adapt to their
let (yi )N
1 be the eigenvalues of GUE with joint density
Pnnotation,
2
2
proportional to exp(− 1 yi )∆ (y), where ∆(y) is the usual Vandermonde. In
√ this scaling
the eigenvalue bulk
p concentrates as the
psemi-circle on [−µN , µN ] with µN = 2N.
We have yi = N/2 λi and µN = N/2 λ+ , for λ+ = 2, so that
TN =

N
X
1

(λi − λ+ )q+ =

N
 2 q/2 X

N

2

1

(yi − µN )q+ .

From the determinantal structure of GUE, the marginal density of a single (unordered)
eigenvalue yi is given by the one-point function
N

−1

SN (y, y) = N

−1

N
−1
X

φ2k (y),

k=0

where φk (y) are the (Hermite) functions obtained by orthonormalizing y k e−y
 2 q/2 Z ∞
E(TN ) =
(y − µN )q SN (y, y)dy.
N
µN

2 /2

. Thus

Now introduce the TW scaling
τN = √

y = µN + τN x,

1
,
2N 1/6

and let Ai denote the Airy function. Tracy and Widom [1996, p 745-6] show that
SτN (x, x) = τN SN (µN + τN x, µN + τN x)
Z ∞
→ KA (x, x) =
Ai2 (x + z)dz,
0

with the convergence being dominated: SτN (x, x) ≤ M 2 e−2x . Consequently,
 2τ 2 q/2 Z ∞
N
xq SτN (x, x)dx
E(TN ) =
N
0
Z ∞
−2q/3
∼N
xq KA (x, x)dx.
0

In particular, E(TN ) = O(N −2q/3 ), and if q = 0, then E(TN ) converges to a positive
constant.
Integration by parts and Olver et al. [2010, 9.11.15] yield
Z ∞
Z ∞ Z ∞
q
cq =
x KA (x, x)dx =
xq
Ai2 (z)dzdx
0
0
x
Z ∞
(3)
1
2Γ(q + 1)
q+1
2
=
x Ai (x)dx = √
q+1 0
π12(2q+9)/6 Γ((2q + 9)/6)
√
For q = 0 the constant becomes c0 = 1/(6 3π).
Remarks. 1. Ullah [1983] states, in our notation, that the expected number of
eigenvalues above the bulk edge, E(TN ) ∼ 0.25N −1/2 . This claim cannot be correct: a
counterexample uses the limiting law F2 for y(1) = maxi yi of Tracy and Widom [1994]:
√
(4)
E(TN ) ≥ Pr(y(1) > 2N ) → 1 − F2 (0) = 0.030627.
We evaluated numerically in Mathematica the formulas (U3), (U6) and (U7) for p =
(2/N)E(TN ) given in Ullah [1983]. While numerical results from intermediate formula (U3)
3

Table 1: For
√ GUE(N), the probabilities pN (k) of exactly k eigenvalues exceeding the upper
bulk edge 2N, along with the expected number E(TN ), to be compared with limiting
value (2).
N
pN (1)
pN (2)
pN (3)
E(TN )
10
25
50
100
250
500

2.868 · 10−2
2.955 · 10−2
2.994 · 10−2
3.019 · 10−2
3.039 · 10−2
3.048 · 10−2

1.36 · 10−6
1.70 · 10−6
1.88 · 10−6
2.00 · 10−6
2.09 · 10−6
2.14 · 10−6

6.9 · 10−14
1.4 · 10−13
1.9 · 10−13
2.3 · 10−13
2.6 · 10−13
2.8 · 10−13

0.028681
0.029551
0.029944
0.030195
0.030392
0.030480

are consistent with our (2), neither those from (U6) nor those from the final result (U7)
are consistent with (U3), or indeed with each other!
2. The striking closeness of the right side of (4) to (2) led us to use the Matlab toolbox
of Bornemann [2010] to evaluate numerically
√
(n)
pN (k) = Pr( exactly k of {yi } > 2N) = E2 (k, J)
√
with J = ( 2N, ∞), in the notation of Bornemann [2010]. The results, in Table 1, confirm
that the probability of 2 or more eigenvalues exiting the bulk is very small, of order 10−6 ,
for all N. This is also suggested by the plots of the densities of y(1) , y(2) , . . . in the scaling
limit in Figure 4 of Bornemann [2010], which itself extends Figure 2 of Tracy and Widom
[1994].

3

The real Wishart case

Suppose λi are eigenvalues of n−1 XX ⊤ for X a N × n matrix with i.i.d. N(0, 1) entries.
√
Assume that γN = N/n → γ ∈ (0, 1]. Set λ(γ) = (1 + γ)2 .
We recall the scaling for the Tracy-Widom law from the largest eigenvalue λ1 :
λ1 = λ(γN ) + N −2/3 τ (γN )WN
where WN converges in distribution to W ∼ T W1 and τ (γ) =

√ √
γ( γ + 1)4/3 .

Theorem 2. (a) Suppose η(λ, c) ≥ 0 is jointly continuous in λ and c, and satisfies
η(λ, c) = 1
η(λ, c) ≤ Mλ

for λ ≤ λ(c)
for some M and all λ.

Suppose also that cN − γN = O(N −2/3 ). Then for q > 0,
E

X
N
i=1

[η(λi , cN ) − 1]
4

q



→ 0.

(5)

(b) Suppose cN − γN ∼ sσ(γ)N −2/3 , where σ(γ) = τ (γ)/λ′ (γ) = γ(1 +
X

Z
N
q
q
−2q/3
E
[λi − λ(cN )]+ ∼ τ (γ)N
i=1

√

γ)1/3 . Then

∞
s

(x − s)q+ K1 (x, x)dx.

(6)

where K1 is defined at (9) below.
(c) In particular, let Nn = #{i : λi ≥ λ(cN )} and suppose that cN − γN = o(N −2/3 ). Then
Z ∞
ENn → c0 =
K1 (x, x)dx ≈ 0.17.
0

Remarks. 1. Part (b) represents a sharpening of (5) that is relevant when η(λ) = η(λ, γ)
is Hölder continuous in λ near the bulk edge λ(γ),
η(λ) − η(λ(γ)) ∼ (λ − λ(γ))q+ .
The example q = 1/2 occurs commonly for optimal shrinkage rules η ∗ (λ) in Donoho et al.
[2017].
2. Section 4 explains why we allow cN to differ from γN .
Proof. Define
TN =

N
X

(
[η(λ, c) − 1]q
F (λ, c) =
[λ − λ(c)]q+

F (λi , cN ),

i=1

(a)
(b).

We adapt the discussion here to the notation used in Tracy and Widom [1998] and Johnstone
[2001]. Let (yi )N
1 = nλi be the eigenvalues of WN (n, I) with joint density function PN (y1 , . . . , yN )
with explicit form given, for example, in [Johnstone, 2001, eq. (4.1)]. We obtain
Z ∞
E(TN ) =
F (y/n, cN )R1 (y)dy,
0

R

where R1 (y1 ) = N (0,∞)N−1 PN (y1 , . . . , yN )dy2 · · · dyN is the one-point (correlation) function. It follows from Tracy and Widom [1998, p814–16] that
R1 (y) = T1 (y) = 21 tr KN (y, y)

(7)

where KN (x, y) is the 2 × 2 matrix kernel associated with PN , see e.g. [Tracy and Widom,
1998, eq. (3.1)]. It follows from Widom [1999] that
1
tr
2

KN (y, y) = S(y, y) + ψ(y)(ǫφ)(y) = S1 (y, y),

(8)

where the functions S(y, y ′), ψ(y) and φ(y) are defined in terms of orthonormalized Laguerre
polynomials in Widom [1999] and studied further in Johnstone [2001]. The function ǫ(x) =
1
sgnx and the operator ǫ denotes convolution with the kernel ǫ(x − y).
2

5

For convergence, introduce the Tracy-Widom scaling
y = µN + σN x,
1
2

and nh = n + 21 and define
p
√
1/3
σN = c(Nh /nh )Nh ,
µN = ( Nh + nh )2 ,
√
√
√
where c(γ) = (1 + γ)1/3 (1 + 1/ γ) = (1 + γ)1/3 λ′ (γ) We now rescale the scalar-valued
function (8):
S1τ (x, x) = σN S1 (µN + σN x, µN + σN x).
where we set Nh = N +

We can rewrite our target E(TN ) using (7), (8) and this rescaling in the form
Z ∞
E(TN ) =
F (ℓN (x), cN )S1τ (x, x)dx,
δN

where ℓN (x) = (µN +σN x)/n, δN = (nλ(cN )−µN )/σN and we used the fact that F (λ, c) = 0
for λ ≤ λ(c).
It follows from [Johnstone, 2001, eq. (3.9)] that


Z ∞
Z ∞
S1τ (x, x) = 2
φτ (x + u)ψτ (x + u)du + ψτ (x) cφ −
φτ (u)du .
0

x

It is shown in equations (3.7), 3.8) and Sec. 5 of that paper that
1
φτ (x), ψτ (x) → √ Ai(x)
2
and, uniformly in N and in intervals of x that are bounded below, that
φτ (x), ψτ (x) = O(e−x ).
√
Along with cφ → 1/ 2 (cf. App. A7 of same paper), this shows that


Z ∞
Z ∞
2
1
S1τ (x, x) → K1 (x, x) =
Ai (x + z)dz + 2 Ai(x) 1 −
Ai(z)dz > 0
0

(9)

x

with the convergence being dominated
S1τ (x, x) ≤ M 2 e−2x + M ′ e−x .

(10)

Before completing the argument for (a) – (c), we note it is easily checked that
n−1 µN = λ(γN ) + O(N −1 ),
so that
δN =

n
[λ(cN ) − λ(γN )] + O(N −1/3 ).
σN

If cN − γN = θN N −2/3 for θN = O(1) then
δN ∼

n −2/3
N
θN λ′ (γ) ∼ θN /σ(γ),
σN
6

(11)

since we have
N 2/3 σN /n ∼ σ(γ)λ′ (γ) = τ (γ).

(12)

In case (a), then, δN ≥ −A for some A. We then have ℓN (x) → λ(γ) for all x ≥ −A,
and so from joint continuity
η(ℓN (x), cN ) → η(λ(γ), γ) = 1,
and hence for all x ≥ −A,
F (ℓN (x), cN ) = [η(ℓN (x), cN ) − 1]q → 0

(13)

The convergence is dominated since the assumption η(λ, c) ≤ Mλ implies that |F (ℓN (x), cN )| ≤
C(1 + |x|q ). Hence the convergence (13) along with (10) and the dominated convergence
theorem implies (5).
For case (b),
Z ∞
2q/3
N
E(TN ) =
[N 2/3 (ℓN (x) − λ(cN ))]q+ S1τ (x, x)dx.
δN

Observe that
N 2/3 (λ(γN ) − λ(cN )) ∼ N 2/3 λ′ (γ)(γN − cN ) ∼ −sτ (γ),
and so from (11) and (12), we have
N 2/3 (ℓN (x)−λ(cN )) = O(N −1/3 )+N 2/3 (λ(γN )−λ(cN ))+N 2/3 n−1 σN x ∼ τ (γ)(x−s). (14)
In addition, from (14), we have
N 2/3 |ℓN (x) − λ(cN )| ≤ M(1 + |x|),
so that the convergence is dominated and (6) is proven.
For case (c), we have only to evaluate
Z ∞
Z ∞
Z
1
c0 =
K1 (x, x)dx =
KA (x, x)dx + 4
0

0

∞

G′ (x)dx = I1 + I2 ,
0

I1 was evaluated in the previous section and G(x) = [1 −
Rwhere
∞
Ai(z)dz = 1/3, from Olver et al. [2010, 9.10.11], we obtain
0
4I2 = G(∞) − G(0) = 1 − (2/3)2 = 5/9,

with the result

5
1
≈ 0.031 + 0.139 = 0.16952.
c0 = √ +
6 3π 36

7

R∞
x

Ai(z)dz]2 . Since

4

Application to covariance estimation

We indicate how Theorem 2 is applied to covariance estimation in the spiked model studied
in Donoho et al. [2017]. Consider a sequence of statistical problems indexed by dimension
p and sample size n. In the nth problem X̌ ∼ Np (0, Σ) where p = pn sastisfies pn /n →
γ ∈ (0, 1] and the population covariance matrix Σ = Σp has fixed ordered eigenvalues
ℓ1 ≥ . . . ≥ ℓr > 1 for all n, and then ℓr+1 = . . . = ℓpn = 1.
Suppose that the sample covariance matrix Š = Šn,pn has eigenvalues λ̌1 ≥ . . . ≥ λ̌p
and corresponding eigenvectors v1 , . . . , vp . Consider shrinkage estimators of the form
Σ̂η =

p
X

η(λ̌j , cp )vj vj⊤ ,

(15)

j=1

where η(λ, c) is a continuous bulk shrinker, that is, satisfies the conditions (a) of Theorem
2. Without loss of generality, as explained in the reference cited, we may also assume that
λ → η(λ, c) is non-decreasing. In the spiked model, the typical choice for cp in practice
would be to set cp = p/n, and we adopt this choice below.
It is useful to analyse an “oracle” or “rank-aware” variant of (15) which takes advantage
of the assumed structure of Σp , especially the fixed rank r of Σp − I:
Σ̂η,r =

r
X

η(λ̌j , cp )vj vj⊤

+

p
X

vj vj⊤ .

j=r+1

j=1

The error in estimation of Σ using Σ̂ is measured by a loss function Lp (Σ, Σ̂). One seeks
conditions under which the losses Lp (Σ, Σ̂η ) and Lp (Σ, Σ̂η,r ) are asymptotically equivalent.
They consider a large class of loss functions which satisfy a Lipschitz condition which
implies that, for some q,
|Lp (Σ, Σ̂η ) − Lp (Σ, Σ̂η,r )| ≤ C(ℓ1 , η(λ̌1))

p
X

[η(λ̌j , cp ) − 1]q .

j=r+1

Suppose now that Π : Rp → Rp−r is a projection on the span of the p − r unit eigenvectors
of Σ. Let X = ΠX̌ and let λ1 ≥ · · · ≥ λp−r denote the eigenvalues of n−1 XX ⊤ . By the
Cauchy interlacing Theorem (e.g. [Bhatia, 1997, p. 59]), we have
λ̌j ≤ λj−r

for r + 1 ≤ j ≤ p,

(16)

where the (λi )p−r
i=1 are the eigenvalues of a white Wishart matrix Wp−r (n, I). From the
monotonicity of η,
p−r
p
X
X
q
[η(λi , cp ) − 1]q .
(17)
[η(λ̌j , cp ) − 1] ≤
i=1

j=r+1

Now apply part (a) of Theorem 2 with the identifications
N ← p − r,

cN ← cp .
8

Clearly γN = N/n → γ and
cN − γ N =

N +r N
−
= O(N −2/3 ),
n
n

since r is fixed. We conclude that the right side of (17) and hence |Lp (Σ, Σ̂η ) − Lp (Σ, Σ̂η,r )|
converge to 0 in L1 and in probability.
Part (c) of Theorem 2 helps to give an example where the losses Lp (Σ, Σ̂η ) and Lp (Σ, Σ̂η,r )
−1
are not asymptotically equivalent. Indeed, let Lp (Σ, Σ̂η ) = kΣ̂−1
η − Σ k, with k · k denoting
matrix operator norm. Here the
shrinkage rule η = η ∗ (λ, c) is discontinuous at the
√ optimal
upper bulk edge λ(c) = (1 + c)2 :
η ∗ (λ, c) = 1
∗

η (λ, c) → 1 +

√

c

for λ ≤ λ(c)

for λ ↓ λ(c).

Proposition 3 of Donoho et al. [2017] shows that
D

−1
−1
−1
kΣ̂−1
η − Σ k − kΣ̂η,r − Σ k → W,

(18)

where W has a two point distribution (1 − π)δ0 + πδw with non-zero
√ probability
π = Pr(T W1 > 0) at location w = f (ℓ+ ) − f (ℓr ), where ℓ+ = 1 + c and the function

1/2
c(ℓ − 1)
f (ℓ) =
ℓ(ℓ − 1 + γ)

is strictly decreasing for ℓ ≥ ℓ+ .
Part (c) of Theorem 2, along with interlacing inequality (16), is used in the proof to
establish that Nn = #{i ≥ r + 1 : λin > λ+ (cn )}, the number of noise eigenvalues exiting
the bulk, is bounded in probability.

5

Final Remarks

It is apparent that the same methods will show that the value of c0 for the Gaussian
Orthogonal Ensemble will be the same as for the real Wishart (Laguerre Orthogonal Ensemble), and similarly that the value of c0 for the white complex Wishart (Laguerre Unitary
Ensemble) will agree with that for GUE.
Some natural questions are left for further work. First, the evaluation of c0 for values
of β other than 1 and 2, and secondly universality, i.e. that the limiting constants do not
require the assumption of Gaussian matrix entries.
Finally, this article appears in a special issue dedicated to the memory of Peter Hall.
Hall’s many contributions to high dimensional data have been reviewed by Samworth [2016].
However, it seems that Peter did not publish specifically on problems connected with the
application of random matrix theory to statistics — the exception that proves the rule of
his extraordinary breadth and depth of interests. Nevertheless the present author’s work
on this specific topic, as well as on many others, has been notably advanced by Peter’s
support — academic, collegial and financial – in promoting research visits to Australia and
contact with specialists there in random matrix theory, particularly at the University of
Melbourne, Peter’s academic home since 2006.
9

References
Z. D. Bai and Jack W. Silverstein. No eigenvalues outside the support of the limiting spectral distribution of large-dimensional sample covariance matrices. Annals of Probability,
26(1):316–345, 1998. ISSN 0091-1798.
Rajendra Bhatia. Matrix Analysis, volume 169 of Graduate Texts in Mathematics. SpringerVerlag, New York, 1997. ISBN 0-387-94846-5.
F. Bornemann. On the numerical evaluation of distributions in random matrix theory:
a review. Markov Processes and Related Fields, 16(4):803–866, 2010. ISSN 1024-2953.
arXiv:0904.1581.
M. Chiani. On the probability that all eigenvalues of Gaussian, Wishart, and double
Wishart random matrices lie within an interval. IEEE Transactions on Information
Theory, 63(7):4521–4531, 2017.
DLMF. NIST Digital Library of Mathematical Functions. http://dlmf.nist.gov/, Release
1.0.9 of 2014-08-29, 2014. Online companion to Olver et al. [2010].
David Donoho, Matan Gavish, and Iain M. Johnstone. Optimal shrinkage of eigenvalues
in the spiked covariance model. arxiv:1311.0851v3; in press, Annals of Statistics, 2017.
Stuart Geman. A limit theorem for the norm of random matrices. Annals of Probability,
8:252–261, 1980.
Iain M. Johnstone. On the distribution of the largest eigenvalue in principal components
analysis. Annals of Statistics, 29:295–327, 2001.
Satya N. Majumdar and Massimo Vergassola. Large deviations of the maximum eigenvalue
for Wishart and Gaussian random matrices. Physical Review Letters, 102:060601, Feb
2009.
Satya N. Majumdar and Pierpaolo Vivo. Number of relevant directions in principal component analysis and Wishart random matrices. Physical Review Letters, 108:200601, May
2012.
Ricardo Marino, Satya N. Majumdar, Grégory Schehr, and Pierpaolo Vivo. Phase transitions and edge scaling of number variance in Gaussian random matrices. Physical Review
Letters, 112:254101, Jun 2014.
F. W. J. Olver, D. W. Lozier, R. F. Boisvert, and C. W. Clark, editors. NIST Handbook
of Mathematical Functions. Cambridge University Press, New York, NY, 2010. Print
companion to DLMF.
Leonid Pastur and Mariya Shcherbina. Eigenvalue Distribution of Large Random Matrices,
volume 171 of Mathematical Surveys and Monographs. American Mathematical Society,
Providence, RI, 2011. ISBN 978-0-8218-5285-9.
Richard J. Samworth. Peter Hall’s work on high-dimensional data and classification. Annals
of Statistics, 44(5):1888–1895, 2016. ISSN 0090-5364.
10

Craig A. Tracy and Harold Widom. Level-spacing distributions and the Airy kernel. Communications in Mathematical Physics, 159:151–174, 1994.
Craig A. Tracy and Harold Widom. On orthogonal and symplectic matrix ensembles.
Communications in Mathematical Physics, 177:727–754, 1996.
Craig A. Tracy and Harold Widom. Correlation functions, cluster functions, and spacing
distributions for random matrices. Journal of Statistical Physics, 92:809–835, 1998.
N Ullah. Number of energy levels outside Wigner’s semicircle. Journal of Physics A:
Mathematical and General, 16(18):L767, 1983.
H. Widom. On the relation between orthogonal, symplectic and unitary ensembles. Journal
of Statistical Physics, 94:347–363, 1999.

11

