Large Random Matrices and Applications to Machine Learning

Contents

• Introduction and overview
• 0. Review: eigenvalues, resolvents and Gaussian vectors
• 1. Spectral measures and the Stieltjes transform
• 2. Gaussian matrix calculus and the Poincaré–Nash inequality
• 3. The Marchenko–Pastur theorem
• 4. Structured models and deterministic equivalents
• 5. Multi-antenna (MIMO) channels: mutual information and ergodic capacity
• 6. Spiked models: outliers and eigenvectors
• 7. The largest eigenvalue and the Tracy–Widom law
• 8. Source detection in high dimension and spectrum sensing
• 9. G-estimation, power estimation and positioning
• 10. Array processing: angles of arrival, MUSIC and G-MUSIC
• 11. Community detection on large graphs
• 12. Kernel random matrices and kernel methods
• 13. Complements for self-study
• 14. Research directions and open problems
• Exercises
• Projects and assessment
• References

Introduction and overview

Audience: a student with a good M1 in mathematics, statistics or signal processing who wants to know why classical estimators fail when the dimension is comparable to the number of samples, and how to correct them.

Empirical first: draw 1000 samples of a 500-dimensional standard Gaussian vector. The true covariance is the identity, so every eigenvalue equals 1. Yet the eigenvalues of the sample covariance spread from 0.08 to 2.87, and its condition number is about 35 (code in ch. 0). This spread is not bad luck: the Marchenko–Pastur law predicts both ends within 0.05, and it biases every method that plugs the sample covariance in.

This course studies the spectra of large random matrices when both dimensions grow at the same rate. It builds the tools (Stieltjes transform, Gaussian calculus, deterministic equivalents), proves the two founding limit theorems, and then uses them to design and analyse estimators and detectors for wireless communications, array processing, graphs and kernel methods.

It shares tools with Kernel Methods for Machine Learning, Graphs in Machine Learning, Detection Theory and Industrial Applications and Theoretical Foundations of Deep Learning.

What you will be able to do

• Compute and invert the Stieltjes transform of a spectral measure, and use it to prove the convergence of empirical spectral distributions.
• Apply Gaussian integration by parts and the Poincaré–Nash inequality to compute expectations and bound variances of traces of resolvents.
• State and prove the Marchenko–Pastur theorem for Gaussian data, and compute its support, density and moments.
• Write the deterministic equivalents of sample covariance, variance-profile and signal-plus-noise models, and solve their fixed-point equations numerically.
• Predict the ergodic capacity of a large multi-antenna channel and explain why it grows linearly with the number of antennas.
• Locate the outlier eigenvalues of a spiked model, compute the alignment of its eigenvectors, and explain the phase transition.
• Use the Tracy–Widom law to set the threshold of a test based on the largest eigenvalue.
• Derive the generalized likelihood ratio test for source detection and predict its power in the high-dimensional regime.
• Build consistent G-estimators of population eigenvalues and powers, and use them to locate a transmitter.
• Correct the MUSIC algorithm for large arrays, detect communities in a dense graph, and analyse kernel spectral clustering with the same spiked-model tools.

Official objectives and syllabus (2026–2027)

The course description on master-mva.com announces an introduction to the theory of large random matrices, and in particular to its applications to signal processing and learning. The only prerequisite it lists is a course in probability.

It states that the methodological and the applied aspects are presented separately on purpose, and that the course integrates the two. The first part covers the theory: the basic techniques, the Marchenko–Pastur and Tracy–Widom theorems, and the structured and spiked models used in statistics.

The second part lists six applications, from multi-antenna channels to kernel methods. The organisation announced is 10 sessions of 3 hours and one exam session. Assessment counts the lab reports for 50% and a final exam, a written test or a project, for 50%.

The Université Paris-Saclay programme table lists the course as “Random matrix theory”, with 24 h of lectures and 6 h of labs (5 ECTS). Its course detail gives only the assessment type, “Evaluation Continue non Intégrale”, and no programme of its own.

Every official topic has a section in these notes. The extensions beyond the syllabus are the Wigner semicircle law used as a first example (ch. 2), the support of the limit spectrum (ch. 4), and two self-study chapters: complements (general entries, moments, free probability, fluctuations of linear statistics, shrinkage) in ch. 13 and research directions in ch. 14.

Where this course fits in the programme

Official prerequisite: a course in probability. In practice the course also uses linear algebra (spectral theorem, min-max), complex analysis at the level of holomorphic functions and residues, and Gaussian vectors; ch. 0 reviews them.

• Refreshers — Probability and Statistics, Linear Algebra, Analysis and Optimization and Scientific Python cover ch. 0 and the code.
• Semester 1 courses that prepare it — Introduction to Statistical Learning (PCA, ridge regression, high-dimensional risk), Computational Statistics (Monte Carlo, used to check every limit), and Theoretical Foundations of Deep Learning (double descent, whose precise curves come from the tools of ch. 4 and ch. 14).
• Related courses of semester 2 — Kernel Methods for Machine Learning (kernel matrices, ch. 12), Graphs in Machine Learning (spectral clustering, stochastic block models, ch. 11), Detection Theory and Industrial Applications (likelihood ratio tests, ch. 8) and Interactions: Complexity, Stochasticity, Extrema and Rare Events (extremes; the Tracy–Widom law of ch. 7 is the extreme-value law of a strongly correlated system).
• After the course — the Research Internship: high-dimensional statistics, wireless communications and the theory of learning all use these tools.

Course outline

• 0. Review — eigenvalue inequalities, resolvents, Schur complements, quadratic forms of random vectors.
• 1. Spectral measures and the Stieltjes transform — the transform, its inversion, the continuity theorem, rank perturbations.
• 2. Gaussian matrix calculus and the Poincaré–Nash inequality — integration by parts, variance bounds, the semicircle law as a first application.
• 3. The Marchenko–Pastur theorem — the limit spectrum of a sample covariance matrix, its edges and moments.
• 4. Structured models and deterministic equivalents — population covariance, variance profile, signal plus noise.
• 5. MIMO channels — mutual information and ergodic capacity of large multi-antenna channels.
• 6. Spiked models — outlier eigenvalues, eigenvector alignment and the phase transition.
• 7. The largest eigenvalue and the Tracy–Widom law — fluctuations at the edge and how to compute the law.
• 8. Source detection and spectrum sensing — the generalized likelihood ratio test and its limits.
• 9. G-estimation, power estimation and positioning — consistent estimators of population quantities, applied to locating a transmitter.
• 10. Array processing — MUSIC, its failure with large arrays, and the G-MUSIC correction.
• 11. Community detection on large graphs — the spectral detectability threshold of the stochastic block model.
• 12. Kernel random matrices — why kernel matrices linearize in high dimension, and what it implies for kernel spectral clustering.
• 13. Complements for self-study — general entries, the moment method, free probability, fluctuations, shrinkage.
• 14. Research directions and open problems — a self-study tour of current work.

Session plan

30 h in class = 10 sessions of 3 h. The volume of 30 h (24 h of lectures, 6 h of labs) comes from the Université Paris-Saclay programme table, and master-mva.com announces 10 sessions of 3 hours plus one exam session, which agrees. The split below is a proposal: neither source gives the content of each session.

The 6 h of labs are placed as four 1.5-hour lab slots, in sessions 3, 5, 8 and 10 (Practice column); their reports are the 50% of lab work in the official assessment. The exam session is not counted in the 30 h.

Level 1 = foundations, 2 = core results and algorithms, 3 = advanced and research-level material. The level never goes down from one session to the next.

Self-study, outside the 30 h (about 7 h in total): ch. 13 Complements for self-study (≈ 4 h) · ch. 14 Research directions and open problems (≈ 3 h).

How to study this module

• Read ch. 0 before session 1 if the min-max theorem or Schur complements feel distant: the proofs use them without warning.
• After each session, do the exercises of its row, then run the code of the chapter and change one parameter (the ratio of dimensions, the spike strength, the sample size) to see the limit appear.
• Check every limit theorem by simulation before trusting it at finite size: the code blocks show how close the asymptotic formulas already are for matrices of size 50 to 1000.
• The code uses NumPy and SciPy only, with synthetic data; every block runs in a few seconds on a laptop.
• The lab slots of sessions 3, 5, 8 and 10 use the code of chapters 1–3, 5, 8–9 and 10–12; the projects of the last section extend them.

0. Review: eigenvalues, resolvents and Gaussian vectors

Every chapter manipulates the eigenvalues of large Hermitian matrices through their resolvents, and averages quadratic forms of random vectors. This chapter collects the M1-level facts used later, with short proofs.

The Linear Algebra, Analysis and Optimization refresher and the Probability and Statistics refresher cover the same ground at a slower pace.

The phenomenon that motivates the course takes five lines of NumPy: white data, a sample covariance, and eigenvalues that are far from 1.

import numpy as np

rng = np.random.default_rng(0)
N, n = 500, 1000                      # dimension, number of samples
X = rng.standard_normal((N, n))       # white data: true covariance = I
S = X @ X.T / n                       # sample covariance
lam = np.linalg.eigvalsh(S)
c = N / n
print(f"smallest / largest eigenvalue: {lam.min():.3f} / {lam.max():.3f}")
print(f"Marchenko-Pastur edges       : {(1 - c**0.5)**2:.3f} / {(1 + c**0.5)**2:.3f}")
print(f"condition number of S        : {lam.max() / lam.min():.1f}")

The output gives eigenvalues between 0.082 and 2.868, against the predicted edges 0.086 and 2.914 of ch. 3, and a condition number of 34.9 instead of 1.

Eigenvalue inequalities

Throughout, A and B are Hermitian N \times N matrices, with eigenvalues in decreasing order \lambda_1(A) \geq \dots \geq \lambda_N(A). The spectral norm is \|A\| = \max_k |\lambda_k(A)|.

Resolvents, Sherman–Morrison and Schur complements

The Schur complement formula inverts a matrix by blocks: if M = \begin{pmatrix} a & b^* \\ b & D \end{pmatrix} with a scalar and D invertible, the top-left entry of the inverse is (M^{-1})_{11} = \frac{1}{a - b^* D^{-1} b} . It gives the diagonal entries of a resolvent in terms of the resolvent of a smaller matrix, which is the starting point of the general-entries proofs of ch. 13.

Quadratic forms of random vectors

The workhorse of the course is the concentration of x^* A x around \mathrm{tr}\, A when x has independent standardized entries and A is independent of x.

1. Spectral measures and the Stieltjes transform

The object of the course is the distribution of the eigenvalues of a large matrix. The Stieltjes transform turns it into an analytic function that equals the normalized trace of the resolvent, and that is what the proofs compute.

The Stieltjes transform also controls finite-rank perturbations, which is the first hint of the spiked models of ch. 6.

So the ESD cannot see a perturbation of fixed rank. Centring the data, for instance, subtracts a rank-one matrix from the sample covariance and does not change its limit (Exercise 2). Individual eigenvalues, however, can move far away: that is the subject of ch. 6.

The code recovers the density of a sample covariance spectrum from its Stieltjes transform, at three values of \varepsilon:

import numpy as np

rng = np.random.default_rng(1)
N, n = 400, 1600
c = N / n
X = rng.standard_normal((N, n))
lam = np.linalg.eigvalsh(X @ X.T / n)

def stieltjes(z):                     # m_N(z) = (1/N) tr (S - z I)^{-1}
    return np.mean(1.0 / (lam - z))

a, b = (1 - c**0.5) ** 2, (1 + c**0.5) ** 2
def mp_density(x):
    return np.sqrt(max((b - x) * (x - a), 0.0)) / (2 * np.pi * c * x)

print("x     eps=0.1  eps=0.03  eps=0.01  MP density")
for x in [0.5, 1.0, 1.5, 2.0]:
    approx = [stieltjes(x + 1j * eps).imag / np.pi for eps in (0.1, 0.03, 0.01)]
    print(f"{x:<5} " + "  ".join(f"{v:7.3f}" for v in approx) + f"  {mp_density(x):9.3f}")

With \varepsilon = 0.1 the estimate is smoothed (0.669 against 0.842 at x = 0.5, near the left edge). With \varepsilon = 0.01, close to the eigenvalue spacing, it starts to follow individual eigenvalues (0.231 against 0.211 at x = 2). The useful range is between the spacing 1/N and the scale of the density.

2. Gaussian matrix calculus and the Poincaré–Nash inequality

For matrices with Gaussian entries, two tools replace most of the combinatorics of random matrix theory: integration by parts computes expectations of traces of resolvents, and the Poincaré–Nash inequality shows that their variance is negligible. This approach goes back to Pastur, and Hachem, Loubaton, Najim and coauthors used it systematically for wireless channels.

Integration by parts and the Poincaré–Nash inequality

Derivatives of resolvents

In the rest of the chapter and in ch. 3, X is an N \times n matrix with independent \mathcal{N}(0, 1) entries and columns x_1, \dots, x_n. We set \Sigma = \frac{1}{n} X X^\top, \qquad Q(z) = (\Sigma - z I_N)^{-1}, \qquad m_N(z) = \frac{1}{N} \mathrm{tr}\, Q(z) .

The variance is of order N^{-2}, not N^{-1}: \mathrm{tr}\, Q is a sum of N terms but fluctuates by O(1) only. This rigidity of spectra is quantified by the central limit theorem of ch. 13.

First application: the semicircle law

A Wigner matrix is the simplest model. Wigner (1955, 1958) found its limit law by the moment method (ch. 13); the Gaussian method proves it in a few lines.

The code checks the two orders of magnitude of the proof: the bias of \mathbb{E}\, m_N is of order 1/N and the variance of order 1/N^2.

import numpy as np

rng = np.random.default_rng(2)

def goe(N):
    G = rng.standard_normal((N, N))
    return (G + G.T) / np.sqrt(2 * N)          # off-diagonal variance 1/N

z = 1.0 + 0.5j
m_sc = (-z + np.sqrt(z - 2) * np.sqrt(z + 2)) / 2   # branch with Im m > 0
print(f"semicircle m(z) = {m_sc:.4f}")
for N in [50, 100, 200, 400]:
    ms = np.array([np.mean(1 / (np.linalg.eigvalsh(goe(N)) - z)) for _ in range(200)])
    print(f"N={N:3d}: mean m_N = {ms.mean():.4f}, "
          f"|bias| = {abs(ms.mean() - m_sc):.1e}, N^2 Var = {N**2 * ms.var():.3f}")

The bias falls from 3.9 \times 10^{-3} to 3.2 \times 10^{-4} as N goes from 50 to 400, and N^2\, \mathrm{Var}\, m_N stays near 1.7.

3. The Marchenko–Pastur theorem

The sample covariance matrix is the basic object of multivariate statistics. When the dimension N and the sample size n grow together, its spectrum does not concentrate at the population eigenvalues: it spreads over an interval whose shape the Marchenko–Pastur theorem gives.

Statement

The data are n independent vectors x_1, \dots, x_n of \mathbb{R}^N with independent entries of mean 0 and variance 1, so that the population covariance is I_N.

As in ch. 2, X = [x_1, \dots, x_n] and \Sigma = \frac{1}{n} X X^\top. We write c_n = N/n and assume N, n \to \infty, \qquad c_n \to c \in (0, \infty) .

The law itself

The bulk is not centred at 1 in a symmetric way: the density is larger on the left, and for c > 1 a fraction 1 - 1/c of the eigenvalues is exactly 0. The data-dependent spread has a simple first consequence: \Sigma is badly conditioned even when the population covariance is the identity.

Extreme eigenvalues

Weak convergence says nothing about a vanishing fraction of eigenvalues. A separate theorem controls the extreme ones, and ch. 6 and ch. 7 rely on it.

Ch. 13 shows how the finite fourth moment enters, and why without it the largest eigenvalue can escape to infinity while the ESD still converges.

The code checks the law, its first two moments and its universality: Gaussian and \pm 1 entries give the same limit.

import numpy as np
from scipy.integrate import quad

rng = np.random.default_rng(3)

def mp_cdf_factory(c):
    a, b = (1 - c**0.5) ** 2, (1 + c**0.5) ** 2
    dens = lambda x: np.sqrt(max((b - x) * (x - a), 0)) / (2 * np.pi * c * x)
    atom = max(1 - 1 / c, 0)
    def cdf(x):
        if x < a:
            return atom if x >= 0 else 0.0
        return atom + quad(dens, a, min(x, b))[0]
    return cdf

for N, n in [(300, 1200), (600, 1200), (1200, 600)]:
    c = N / n
    X = rng.standard_normal((N, n))            # Gaussian entries
    B = rng.choice([-1.0, 1.0], size=(N, n))   # Rademacher entries
    cdf = mp_cdf_factory(c)
    for name, Z in [("gauss", X), ("rademacher", B)]:
        lam = np.sort(np.linalg.eigvalsh(Z @ Z.T / n))
        grid = np.linspace(1e-9, (1 + c**0.5) ** 2 + 0.2, 200)
        emp = np.searchsorted(lam, grid, side="right") / N
        ks = max(abs(e - cdf(g)) for e, g in zip(emp, grid))
        print(f"c={c:.2f} {name:10s}: mean {lam.mean():.3f} (1), "
              f"2nd moment {np.mean(lam**2):.3f} ({1 + c:.2f}), sup|F_N - F| = {ks:.3f}")

In all six cases the Kolmogorov distance to \mu_c is below 0.01, and the second moments match 1 + c within 0.005, including c = 2 with its atom at 0.

4. Structured models and deterministic equivalents

Real data are not white. This chapter treats the three structured models of the official syllabus: a population covariance matrix, a variance profile, and a deterministic signal plus noise. In each case the limit spectrum is described by a fixed-point equation, and a deterministic equivalent replaces the random resolvent by a deterministic matrix.

Population covariance matrix

The observations are y_j = C^{1/2} x_j, where the x_j are as in ch. 3 and C is a deterministic positive semidefinite N \times N matrix with bounded norm. The sample covariance is \hat C = \frac{1}{n} \sum_{j=1}^n y_j y_j^* = \frac{1}{n}\, C^{1/2} X X^* C^{1/2} .

With C = I, \bar Q = m I with \delta = c m, and the equation becomes the Marchenko–Pastur equation of Theorem 3.1. The word “equivalent” matters: when the spectrum of C has no limit, \bar Q has no limit either, but Q - \bar Q still vanishes in the sense of the theorem.

When the ESD of C converges to a law H, the limit is described by the companion transform. Let \underline{m} be the Stieltjes transform of the limit ESD of the n \times n matrix \frac{1}{n} Y^* Y, which has the same non-zero eigenvalues as \hat C.

This is the equation of Marchenko and Pastur (1967) in the form of Silverstein and Bai (1995). Its explicit inverse gives the support.

When H has several atoms and c is small, the support splits into clusters, one per population eigenvalue. Bai and Silverstein (1999) show that the sample eigenvalues then split exactly as the population ones do. This separation condition is what ch. 9 needs to estimate each population eigenvalue.

import numpy as np

rng = np.random.default_rng(4)
N, n = 1000, 10000                    # c = N/n = 0.1
c = N / n
t, w = np.array([1.0, 4.0]), np.array([0.5, 0.5])   # population spectrum H
pop = np.repeat(t, N // 2)

def det_equiv(z, iters=500):
    """Fixed point delta = (1/n) tr C (C/(1+delta) - z)^{-1}; returns (1/N) tr Qbar(z)."""
    delta = 1.0 + 0j
    for _ in range(iters):
        delta = 0.5 * delta + 0.5 * np.sum(pop / (pop / (1 + delta) - z)) / n
    return np.mean(1 / (pop / (1 + delta) - z))

def support_edges(c, t, w):
    """Edges of the limit support: x(u) = -1/u + c sum w t/(1+tu) is increasing off the support."""
    u = -np.logspace(-3, 3, 200000)
    x = -1 / u + c * np.sum(w * t / (1 + np.outer(u, t)), axis=1)
    dx = 1 / u**2 - c * np.sum(w * t**2 / (1 + np.outer(u, t)) ** 2, axis=1)
    flips = np.nonzero(np.diff(dx > 0))[0]
    return np.sort(x[flips])

Y = np.sqrt(pop)[:, None] * rng.standard_normal((N, n))   # columns ~ N(0, C)
lam = np.linalg.eigvalsh(Y @ Y.T / n)
for z in [0.5 + 0.1j, 2.0 + 0.1j, 4.0 + 0.1j]:
    print(f"z={z}: m_N(z) = {np.mean(1 / (lam - z)):.4f}, deterministic = {det_equiv(z):.4f}")
print("support edges (theory):", np.round(support_edges(c, t, w), 3))
print(f"sample eigenvalue range: [{lam.min():.3f}, {lam.max():.3f}], "
      f"fraction below 2: {np.mean(lam < 2):.3f}")

The empirical and deterministic Stieltjes transforms agree to 5 \times 10^{-4} at z = 2 + 0.1i, in the gap, and to within 10^{-2} inside the clusters, where \Im z = 0.1 resolves fewer eigenvalues. The sample eigenvalues range over [0.552, 5.997], inside the predicted support up to the edge fluctuations of ch. 7.

Variance profile

In a wireless channel or a bipartite network, each entry may have its own variance. Let Y_{ij} = \frac{\sigma_{ij}}{\sqrt{n}} X_{ij} with bounded deterministic \sigma_{ij} and standardized independent X_{ij}.

With \sigma_{ij} = 1, t_i = m_c and \tilde t_j = \underline m: the system collapses to the Marchenko–Pastur equation. When \sigma_{ij}^2 = d_i depends on the row only, it is the population model with C = \mathrm{diag}(d_i) (Exercise 11). The proof follows the scheme of Theorem 3.1, with one equation per row and one per column; Hachem, Loubaton and Najim (2007) also give the deterministic equivalents of the mutual information used in ch. 5.

Signal plus noise

In the information-plus-noise model, a deterministic signal matrix A is observed in white noise: \Sigma = \frac{1}{n} (A + \sigma X)(A + \sigma X)^* . Its columns are not centred, so it is not a sample covariance matrix, but it is the right model for deterministic sources in ch. 8–10.

With A = 0 and \sigma = 1 the equation is again the Marchenko–Pastur equation. The three models share one structure: a self-consistent equation for a few scalar functions, solved numerically by damped fixed-point iterations, as in det_equiv above.

5. Multi-antenna (MIMO) channels: mutual information and ergodic capacity

A link with n transmit and N receive antennas carries more information than a single-antenna link, and the question is how much more. The mutual information is a log-determinant of a random matrix, so its large-system limit is a functional of the spectra of ch. 3 and ch. 4.

When the transmitter does not know H, Telatar (1999) shows that for Rayleigh fading the isotropic Gaussian input is optimal. The ergodic capacity is then the average \mathbb{E}\, I_N(\rho) over the fading, the rate achievable by coding over many independent channel realizations.

The large-system limit

By Proposition 5.2, \frac{1}{N} I_N(\rho) is the integral of \log(1 + \rho \lambda) against the ESD of \frac{1}{n} HH^*, and the integrand is unbounded only logarithmically. Theorem 3.1 therefore gives the limit of the capacity per receive antenna, and the deterministic equivalent gives it in closed form.

For R = I and c = N/n the limit per receive antenna is the Shannon transform of the Marchenko–Pastur law: \frac{1}{N}\, \bar I_N(\rho) = \log\Big(1 + \frac{\rho}{1 + \delta}\Big) + \frac{1}{c} \log(1 + \delta) - \frac{\delta}{c\, (1 + \delta)} = \int \log(1 + \rho \lambda)\, \mu_c(d\lambda) .

Foschini and Gans (1998) and Telatar (1999) observed this linear growth; the closed form of the Shannon transform of \mu_c appears in Verdú and Shamai (1999), for random spreading in CDMA.

At high SNR, \mathcal{V}_c(\rho) \approx \min(1, 1/c) \log \rho up to a constant: the slope \min(N, n) of I_N in \log \rho is the multiplexing gain. Receive correlation lowers the constant but not the slope as long as R is invertible.

The code compares Monte Carlo averages with Theorem 5.3, for independent and correlated antennas (R_{ij} = r^{|i - j|}), down to N = n = 2:

import numpy as np

rng = np.random.default_rng(6)

def det_equiv_mi(R, n, rho, iters=500):
    """Deterministic equivalent of E log det(I + rho/n R^{1/2} X X^* R^{1/2}), in nats."""
    N, s2 = R.shape[0], 1 / rho
    ev = np.linalg.eigvalsh(R)
    d = 1.0
    for _ in range(iters):
        d = 0.5 * d + 0.5 * np.sum(ev / (ev / (1 + d) + s2)) / n
    return (np.sum(np.log(1 + ev / (s2 * (1 + d))))
            + n * np.log(1 + d) - n * d / (1 + d))

def monte_carlo_mi(R, n, rho, reps=4000):
    N = R.shape[0]
    L = np.linalg.cholesky(R)
    out = []
    for _ in range(reps):
        X = (rng.standard_normal((N, n)) + 1j * rng.standard_normal((N, n))) / np.sqrt(2)
        H = L @ X
        out.append(np.linalg.slogdet(np.eye(N) + rho / n * H @ H.conj().T)[1])
    return np.mean(out), np.std(out)

rho = 10.0                                   # SNR of 10 dB
for N in [2, 4, 8]:
    for r in [0.0, 0.8]:                     # receive correlation R_ij = r^|i-j|
        R = r ** np.abs(np.subtract.outer(np.arange(N), np.arange(N)))
        mc, sd = monte_carlo_mi(R, N, rho)
        de = det_equiv_mi(R, N, rho)
        print(f"N=n={N}, r={r}: Monte Carlo {mc / np.log(2):6.3f} bits "
              f"(sd {sd / np.log(2):.2f}), deterministic equivalent {de / np.log(2):6.3f} bits")

Even for N = n = 2 the error is below 0.1 bit, and for N = n = 8 it is about 0.015 bit. The standard deviation of I_N stays between 1.0 and 1.3 bits whatever N: the mutual information fluctuates by O(1), as Proposition 2.4 suggests, and ch. 13 gives the Gaussian limit that predicts outage probabilities.

6. Spiked models: outliers and eigenvectors

A low-rank signal does not change the limit spectrum (Proposition 1.6), but it can push a few eigenvalues out of the bulk. Whether it does depends on its strength, with a sharp threshold. This phase transition decides whether principal component analysis, detection or clustering can see a signal at all.

The proofs reduce the location of an outlier to a scalar equation involving a quadratic form of the resolvent. The key input is that such forms behave like normalized traces.

The spiked Wigner model

The spiked covariance model

Two consequences matter for every application of the course.

• Below the threshold, the top sample eigenvector carries asymptotically no information about u: no method based on it can estimate the signal direction.
• Above it, the sample eigenvalue overestimates 1 + \ell and the eigenvector is only partly aligned. Both biases are known functions of \ell and c, so they can be inverted.

The inversion is explicit: \ell + c/\ell = \hat\lambda - 1 - c, and the larger root is \hat\ell = \frac{(\hat\lambda - 1 - c) + \sqrt{(\hat\lambda - 1 - c)^2 - 4c}}{2} . Ch. 10 uses it to correct the MUSIC algorithm.

import numpy as np

rng = np.random.default_rng(7)
N, n = 400, 800
c = N / n                                   # threshold sqrt(c) = 0.707
u = np.zeros(N); u[0] = 1.0                 # spike direction

print(" ell   lambda_max  theory   |<u,u_hat>|^2  theory")
for ell in [0.3, 0.6, 1.0, 2.0, 4.0]:
    lam_max, align = [], []
    for _ in range(20):
        Y = rng.standard_normal((N, n))
        Y[0] *= np.sqrt(1 + ell)            # population covariance I + ell u u^T
        lam, vec = np.linalg.eigh(Y @ Y.T / n)
        lam_max.append(lam[-1]); align.append((vec[:, -1] @ u) ** 2)
    if ell > np.sqrt(c):
        th_l, th_a = (1 + ell) * (1 + c / ell), (1 - c / ell**2) / (1 + c / ell)
    else:
        th_l, th_a = (1 + np.sqrt(c)) ** 2, 0.0
    print(f"{ell:4.1f}   {np.mean(lam_max):9.3f}  {th_l:6.3f}   {np.mean(align):12.3f}  {th_a:6.3f}")

Above the threshold the averages match the limits within 0.07 for the eigenvalues and 0.02 for the alignments. Near it (\ell = 0.6 and 1.0) the finite-size effects are visible: the alignment is 0.05 instead of 0 and 0.35 instead of 0.33, because the transition is smoothed over a window of width of order N^{-1/3}.

7. The largest eigenvalue and the Tracy–Widom law

Tests and detectors compare the largest eigenvalue with a threshold, so they need its fluctuations, not only its limit. They are of order n^{-2/3}, much smaller than the n^{-1/2} of classical statistics, and their law is the Tracy–Widom distribution.

Tracy and Widom (1994, 1996) derived these laws for the Gaussian unitary and orthogonal ensembles. Their means are about -1.77 (\beta = 2) and -1.21 (\beta = 1), their standard deviations about 0.90 and 1.27.

Both laws are asymmetric, with a heavier right tail: their skewness is positive, about 0.22 for \beta = 2 and 0.29 for \beta = 1.

To leading order in the exponent, F_\beta(s) decays like e^{-\beta |s|^3/24} as s \to -\infty. On the right, 1 - F_\beta(s) decays only like e^{-\frac{2\beta}{3} s^{3/2}}.

Since \sigma_{Nn} / n is of order n^{-2/3}, the largest eigenvalue of \Sigma = \frac{1}{n} X X^\top is (1 + \sqrt{c})^2 + O(n^{-2/3}).

The limit is universal. Soshnikov (2002) proved it for symmetric sub-Gaussian entries when n - N = O(N^{1/3}), Péché (2009) removed this restriction, and Pillai and Yin (2014) treated general entries with enough moments. El Karoui (2007) extended it to complex data with a general population covariance, with centring and scaling computed from the equation of ch. 4.

In the spiked model of ch. 6, the transition also changes the fluctuations. For complex Gaussian data, Baik, Ben Arous and Péché (2005) showed that the top eigenvalue keeps Tracy–Widom fluctuations below the threshold, and has Gaussian fluctuations of order n^{-1/2} above it.

Computing the law

The Fredholm determinant is the practical route. Bornemann (2010, Mathematics of Computation) showed that a Gauss–Legendre discretization with a few dozen nodes gives the distribution function to near machine precision, and his review of the same year tabulates the moments of F_1 and F_2. For \beta = 1, the formula of Ferrari and Spohn (2005) writes F_1(s) = \det(I - B_s)_{L^2(0, \infty)}, \qquad B_s(x, y) = \frac{1}{2}\, \mathrm{Ai}\Big(\frac{x + y}{2} + s\Big) .

import numpy as np
from scipy.special import airy

# Tracy-Widom (beta = 1) CDF by Nystrom discretization of a Fredholm determinant
# (Ferrari-Spohn formula F1(s) = det(I - B_s) on L^2(0, inf), B_s(x, y) = Ai((x+y)/2 + s)/2)
xg, wg = np.polynomial.legendre.leggauss(60)
L = 12.0                                     # truncate (0, inf) to (0, L)
x, w = L * (xg + 1) / 2, L * wg / 2

def F1(s):
    B = 0.5 * airy((x[:, None] + x[None, :]) / 2 + s)[0]
    sw = np.sqrt(w)
    return np.linalg.det(np.eye(len(x)) - sw[:, None] * B * sw[None, :])

s = np.linspace(-8, 6, 1401)
F = np.array([F1(t) for t in s])
dens = np.gradient(F, s)
mean = np.trapezoid(s * dens, s)
sd = np.sqrt(np.trapezoid((s - mean) ** 2 * dens, s))
q95 = s[np.searchsorted(F, 0.95)]
print(f"TW1: mean {mean:.3f}, sd {sd:.3f}, 95% quantile {q95:.3f}")

# Largest eigenvalue of real white Wishart matrices, Johnstone (2001) centring and scaling
rng = np.random.default_rng(8)
N, n = 100, 400
mu = (np.sqrt(n - 1) + np.sqrt(N)) ** 2
sigma = (np.sqrt(n - 1) + np.sqrt(N)) * (1 / np.sqrt(n - 1) + 1 / np.sqrt(N)) ** (1 / 3)
stats = []
for _ in range(2000):
    X = rng.standard_normal((N, n))
    stats.append((np.linalg.eigvalsh(X @ X.T)[-1] - mu) / sigma)
stats = np.array(stats)
print(f"simulation N={N}, n={n}: mean {stats.mean():.3f}, sd {stats.std():.3f}, "
      f"P(stat > {q95:.2f}) = {np.mean(stats > q95):.3f}")

The computed law has mean -1.206, standard deviation 1.267 and 95% quantile 0.980.

For N = 100 and n = 400, the standardized simulated eigenvalue has mean -1.29 and standard deviation 1.21, and exceeds 0.98 in 4.1% of the runs: a test at nominal level 5% is slightly conservative at this size.

8. Source detection in high dimension and spectrum sensing

A cognitive radio may use a frequency band only when the licensed (primary) user is silent. It must decide, from n samples received on N antennas, whether a source is present. With N and n of the same order and an unknown noise level, the tests of ch. 6 and ch. 7 are the right tools.

The test

Under the null hypothesis \mathcal{H}_0 the samples are noise; under \mathcal{H}_1 a single source with unknown channel h \in \mathbb{C}^N is present: \mathcal{H}_0 : y_t = \sigma w_t, \qquad \mathcal{H}_1 : y_t = h\, s_t + \sigma w_t, \qquad t = 1, \dots, n , with independent standard complex Gaussian s_t and w_t, and unknown \sigma^2. The population covariance is R = h h^* + \sigma^2 I under \mathcal{H}_1.

The signal-to-noise ratio is \rho = \|h\|^2 / \sigma^2, and the sample covariance is \hat R = \frac{1}{n} \sum_t y_t y_t^*.

T_N is pivotal: multiplying the data by a constant does not change it, so its law under \mathcal{H}_0 does not depend on \sigma^2, and the threshold can be set without knowing the noise level.

Three practical conclusions follow.

• Threshold. For a false-alarm probability \alpha, reject when T_N > \mu_{Nn}/n + q_{1-\alpha}\, \sigma_{Nn}/n, with the Tracy–Widom quantile q_{1-\alpha} (real data use F_1, as in the code below).
• Detectability. A source with \rho < \sqrt{c} is asymptotically undetectable by this test: its power tends to the level \alpha. With N = 50 antennas and n = 100 samples, the limit is \rho = 0.707, or -1.5 dB.
• Noise uncertainty. The energy detector, which compares \frac{1}{N} \mathrm{tr}\, \hat R with a threshold computed from \sigma^2, is optimal when \sigma^2 is known. Under an error on \sigma^2, its false-alarm rate tends to 0 or 1: this fragility is behind the “SNR wall” of Tandra and Sahai (2008).

Other eigenvalue-based detectors use the ratio \hat\lambda_{\max} / \hat\lambda_{\min}, whose limit under \mathcal{H}_0 is (1 + \sqrt{c})^2 / (1 - \sqrt{c})^2 (Cardoso, Debbah, Bianchi and Najim, 2008; Zeng and Liang, 2009). Estimating the number of sources uses the same ingredients, one eigenvalue at a time (Nadakuditi and Edelman, 2008; Kritchman and Nadler, 2009).

The code compares the GLRT with an energy detector calibrated for \sigma^2 = 1, when the noise variance is exactly 1 and when it is 1.1, a 0.4 dB error:

import numpy as np

rng = np.random.default_rng(9)
N, n = 50, 100                               # sensors, samples: c = 0.5
c = N / n
q95 = 0.98                                   # 95% quantile of TW1 (ch. 7)
mu = (np.sqrt(n - 1) + np.sqrt(N)) ** 2 / n
sig = (np.sqrt(n - 1) + np.sqrt(N)) * (1 / np.sqrt(n - 1) + 1 / np.sqrt(N)) ** (1 / 3) / n
thr_glrt = mu + q95 * sig                    # threshold on T = lambda_max / (tr R / N)
thr_energy = 1 + 1.645 * np.sqrt(2 / (N * n))  # energy test calibrated for sigma^2 = 1

def statistics(snr, noise_var=1.0):
    h = np.ones(N) * np.sqrt(snr * noise_var / N)          # ||h||^2 / sigma^2 = snr
    Y = np.outer(h, rng.standard_normal(n)) + np.sqrt(noise_var) * rng.standard_normal((N, n))
    R = Y @ Y.T / n
    return np.linalg.eigvalsh(R)[-1] / (np.trace(R) / N), np.trace(R) / N

for snr, nv in [(0.0, 1.0), (0.0, 1.1), (0.4, 1.0), (0.8, 1.0), (1.2, 1.0), (2.0, 1.0)]:
    T = np.array([statistics(snr, nv) for _ in range(1000)])
    print(f"snr={snr:.1f}, noise variance={nv}: P(GLRT detects) = {np.mean(T[:, 0] > thr_glrt):.3f}, "
          f"P(energy detects) = {np.mean(T[:, 1] > thr_energy):.3f}, mean T = {T[:, 0].mean():.3f}")
print(f"GLRT threshold {thr_glrt:.3f}; (1+sqrt c)^2 = {(1 + c**0.5)**2:.3f}; "
      f"detectability threshold sqrt(c) = {c**0.5:.3f}")

The GLRT keeps its false-alarm rate near 0.04 whether the noise variance is 1 or 1.1, while the energy detector jumps from 0.046 to 1.000. Below the threshold (\rho = 0.4) the GLRT detects in 5.6% of the runs, about its level. Above it, its power rises to 0.55 at \rho = 1.2 and 0.96 at \rho = 2, against 0.33 and 0.61 for the energy detector.

9. G-estimation, power estimation and positioning

The spectrum of \hat C is a distorted image of the spectrum of C, and ch. 4 gives the distortion explicitly. G-estimation, introduced by Girko, inverts it: it builds estimators of population quantities that are consistent when N and n grow together. This chapter derives the main estimators and uses them to locate a transmitter from received powers.

Two simple G-estimators

For C = I and c = 0.5, the naive estimate of \frac{1}{N}\mathrm{tr}\, C^{-1} = 1 is close to 2: plugging \hat C^{-1} into a linear discriminant or a Mahalanobis distance doubles it.

Estimating population eigenvalues

Suppose that C has K distinct eigenvalues \gamma_1 < \dots < \gamma_K with known multiplicities N_1, \dots, N_K, N_k / N fixed, and c < 1.

Let \hat\lambda_1 \leq \dots \leq \hat\lambda_N be the sample eigenvalues, and \mathcal{N}_k the indices of the k-th group of N_k of them.

The naive estimator, the average of the sample eigenvalues of the k-th cluster, is biased: the clusters spread asymmetrically (Example 4.4), and the largest population eigenvalue is overestimated while the smallest is underestimated.

Power estimation and positioning by trilateration

A sensor with N antennas observes a transmitter through an M-dimensional subspace. As a teaching model, take the received covariance R = \sigma^2 I_N + P\, U U^*, \qquad U^* U = I_M , where P is the received power. The population eigenvalues are \sigma^2 + P, with multiplicity M, and the noise level \sigma^2 for the rest. Theorem 9.3 estimates both, hence the power.

With a path-loss model P = P_0\, d^{-\alpha} and known P_0 and \alpha, the received power gives the distance d = (P_0 / P)^{1/\alpha}.

Three sensors at known positions a_1, a_2, a_3 then locate the transmitter x. Subtracting the first circle equation \|x - a_1\|^2 = d_1^2 from the others gives a linear system: 2\, (a_s - a_1)^\top x = d_1^2 - d_s^2 + \|a_s\|^2 - \|a_1\|^2, \qquad s = 2, 3 .

Couillet, Silverstein, Bai and Debbah (2011) treat the realistic version, where several sources with their own antennas and random channels share the band. Their estimator of each source power comes from the same contour-integral argument.

import numpy as np

rng = np.random.default_rng(10)

def mestre(lam, n, sizes):
    """Mestre (2008) estimator of population eigenvalues with known multiplicities (N < n)."""
    lam = np.sort(lam)
    s = np.sqrt(lam)
    mu = np.sort(np.linalg.eigvalsh(np.diag(lam) - np.outer(s, s) / n))
    out, start = [], 0
    for Nk in sizes:
        out.append(n / Nk * np.sum(lam[start:start + Nk] - mu[start:start + Nk]))
        start += Nk
    return np.array(out)

N, n, M = 64, 256, 16                  # antennas, samples, signal dimension
sizes = [N - M, M]                     # noise eigenvalue 1 (sigma^2), signal eigenvalue 1 + P
anchors = np.array([[0.0, 0.0], [100.0, 0.0], [0.0, 100.0]])
target = np.array([40.0, 30.0])
P0, alpha = 9000.0, 2.0                # received power P = P0 / d^alpha (noise power 1)

def received_power_estimates(P):
    pop = np.repeat([1.0, 1.0 + P], sizes)
    Y = np.sqrt(pop)[:, None] * rng.standard_normal((N, n))
    lam = np.sort(np.linalg.eigvalsh(Y @ Y.T / n))
    naive = lam[-M:].mean() - lam[:-M].mean()
    g = mestre(lam, n, sizes)
    return naive, g[1] - g[0]

def trilaterate(d):
    A = 2 * (anchors[1:] - anchors[0])
    b = d[0]**2 - d[1:]**2 + np.sum(anchors[1:]**2, 1) - np.sum(anchors[0]**2)
    return np.linalg.lstsq(A, b, rcond=None)[0]

dist = np.linalg.norm(anchors - target, axis=1)
P = P0 / dist**alpha
print("true distances:", np.round(dist, 1), " received SNR:", np.round(P, 2))
err = {"naive": [], "G-estimator": []}
for _ in range(200):
    est = np.array([received_power_estimates(p) for p in P])      # rows: anchors
    for k, name in enumerate(err):
        d_hat = (P0 / np.maximum(est[:, k], 1e-3)) ** (1 / alpha)
        err[name].append(np.linalg.norm(trilaterate(d_hat) - target))
for name, e in err.items():
    print(f"{name:12s}: median position error {np.median(e):5.2f} m, mean {np.mean(e):5.2f} m")

The naive error is a bias, not noise: it does not shrink with more runs, and it grows with c = N/n. The G-estimator removes it, and its remaining error decreases like 1/N.

10. Array processing: angles of arrival, MUSIC and G-MUSIC

An array of antennas estimates the directions of the sources it receives. The MUSIC algorithm reads them from the noise subspace of the sample covariance. With large arrays and few snapshots, the sample subspaces are biased (ch. 6), MUSIC loses its resolution, and a random-matrix correction restores it.

MUSIC replaces \Pi by the projector \hat\Pi on the N - K smallest sample eigenvectors, and returns the K deepest local minima of \hat\eta(\theta) = a(\theta)^* \hat\Pi\, a(\theta).

It is consistent when n \to \infty with N fixed.

The large-array regime

Let the number of sources K be fixed while N/n \to c. The matrix R/\sigma^2 is then a spiked covariance (ch. 6) whose spikes \ell_1, \dots, \ell_K are the eigenvalues of A P A^* / \sigma^2, with eigenvectors u_k.

When two sources are close, their two spikes are unequal (the steering vectors are nearly collinear) and the smaller one has a small \zeta. The two minima of \hat\eta merge, and MUSIC sees one source.

This estimator is the spiked form studied by Hachem, Loubaton, Mestre, Najim and Vallet (2013). The original G-MUSIC of Mestre and Lagunas (2008) weights all sample eigenvectors and does not assume few sources; Vallet, Loubaton and Mestre (2012) treat deterministic source signals with the signal-plus-noise model of ch. 4.

import numpy as np
from scipy.signal import find_peaks

rng = np.random.default_rng(11)

def steer(theta_deg, N):                     # half-wavelength ULA, unit-norm steering vectors
    m = np.arange(N)[:, None]
    return np.exp(1j * np.pi * m * np.sin(np.deg2rad(np.atleast_1d(theta_deg)))) / np.sqrt(N)

def doa_trial(N, n, thetas, P, grid):
    K, c = len(thetas), N / n
    S = np.sqrt(np.array(P))[:, None] * (rng.standard_normal((K, n)) + 1j * rng.standard_normal((K, n))) / np.sqrt(2)
    W = (rng.standard_normal((N, n)) + 1j * rng.standard_normal((N, n))) / np.sqrt(2)
    Y = steer(thetas, N) @ S + W
    lam, U = np.linalg.eigh(Y @ Y.conj().T / n)
    lam, U = lam[::-1], U[:, ::-1]
    x = lam[:K] / lam[K:].mean() - 1 - c      # invert lambda = (1 + l)(1 + c / l)
    ell = (x + np.sqrt(np.maximum(x**2 - 4 * c, 0))) / 2
    zeta = np.maximum((1 - c / ell**2) / (1 + c / ell), 0.05)   # eigenvector shrinkage
    proj = np.abs(U[:, :K].conj().T @ steer(grid, N)) ** 2
    out = {}
    for name, cost in [("MUSIC", 1 - proj.sum(0)), ("G-MUSIC", 1 - (proj / zeta[:, None]).sum(0))]:
        pk, _ = find_peaks(-cost)
        out[name] = np.sort(grid[pk[np.argsort(cost[pk])[:K]]]) if len(pk) >= K else None
    return out

N, n, grid = 40, 80, np.linspace(-30, 50, 4001)
for thetas, P in [([10.0, 16.0], [3.0, 3.0]), ([10.0, 12.0], [3.0, 3.0])]:
    hits = {"MUSIC": [], "G-MUSIC": []}
    for _ in range(200):
        for name, est in doa_trial(N, n, thetas, P, grid).items():
            ok = est is not None and np.all(np.abs(est - thetas) < abs(thetas[1] - thetas[0]) / 2)
            hits[name].append(est - thetas if ok else None)
    for name, h in hits.items():
        good = np.array([e for e in h if e is not None])
        rmse = np.sqrt(np.mean(good**2)) if len(good) else np.nan
        print(f"sources at {thetas}: {name:7s} resolves {len(good) / len(h):.2f} of trials, "
              f"RMSE when resolved {rmse:.3f} deg")

The floor 0.05 on \hat\zeta_k guards against spikes estimated near the threshold, where the correction would divide by almost zero. When a spike is truly below \sqrt{c}, no subspace method can recover its direction (ch. 6), and G-MUSIC does not change that.

11. Community detection on large graphs

Spectral clustering splits a graph with the top eigenvectors of its adjacency matrix. On a random graph with planted communities, the adjacency matrix is a low-rank signal plus a Wigner-like noise, so ch. 6 says exactly when the eigenvectors see the communities.

Off the diagonal, the mean of the adjacency matrix is \mathbb{E} A = \bar p\, \mathbf{1}\mathbf{1}^\top + \frac{p_{\mathrm{in}} - p_{\mathrm{out}}}{2}\, \sigma\sigma^\top . After subtracting the first term and scaling, the centred adjacency matrix decomposes as B = \frac{A - \bar p\, \mathbf{1}\mathbf{1}^\top}{\sqrt{n \bar p (1 - \bar p)}} = W + \theta\, v v^\top + E, \qquad v = \frac{\sigma}{\sqrt{n}} . Here W has independent centred entries of variance \frac{1}{n}(1 + O(n^{-1/2})), and E collects the diagonal correction, of norm O(n^{-1/2}).

Nadakuditi and Newman (2012) pointed out this transition for community detection. Its message is that a spectral method fails below the threshold, even though the communities exist.

The sign of \hat v_i classifies vertex i. The entries of \sqrt{n}\, \hat v are approximately \sqrt{1 - 1/\theta^2}\, \sigma_i plus independent Gaussian noise of variance 1/\theta^2, a heuristic that the simulation confirms. The misclassification rate is then \varepsilon(\theta) \approx \Phi\big(-\sqrt{\theta^2 - 1}\big) , with \Phi the standard normal distribution function.

import numpy as np
from scipy.stats import norm

rng = np.random.default_rng(12)
n, pbar = 2000, 0.1
sigma = np.repeat([1.0, -1.0], n // 2)                 # community labels

print("theta  top eig  theory  overlap^2  theory  error rate  theory")
for theta in [0.8, 1.2, 1.5, 2.0, 3.0]:
    delta = 2 * theta * np.sqrt(pbar * (1 - pbar) / n)  # p_in - p_out
    P = pbar + delta / 2 * np.outer(sigma, sigma)
    U = rng.random((n, n))
    A = np.triu(U < P, 1).astype(float)
    A = A + A.T                                          # undirected, no self-loops
    p_hat = A.sum() / (n * (n - 1))
    B = (A - p_hat) / np.sqrt(n * p_hat * (1 - p_hat))   # centred, scaled adjacency
    lam, vec = np.linalg.eigh(B)
    v = vec[:, -1]
    overlap = (v @ sigma) ** 2 / n
    err = min(np.mean(np.sign(v) != sigma), np.mean(np.sign(v) != -sigma))
    if theta > 1:
        th = (theta + 1 / theta, 1 - 1 / theta**2, norm.cdf(-np.sqrt(theta**2 - 1)))
    else:
        th = (2.0, 0.0, 0.5)
    print(f"{theta:4.1f}  {lam[-1]:7.3f}  {th[0]:6.3f}  {overlap:9.3f}  {th[1]:6.3f}"
          f"  {err:10.3f}  {th[2]:6.3f}")

At \theta = 0.8 the classification is a coin toss (error 0.496). At \theta = 2 the error is 0.043 against 0.042 predicted, and at \theta = 3 it is 0.002.

The sparse regime

Real networks are sparse: p_{\mathrm{in}} = a/n and p_{\mathrm{out}} = b/n with fixed a, b. Formally, the dense threshold \theta = 1 becomes (a - b)^2 > 2\, (a + b) , which is the Kesten–Stigum threshold. Decelle, Krzakala, Moore and Zdeborová (2011) conjectured that it is the threshold for any method, and Mossel, Neeman and Sly (2015, 2018) and Massoulié (2014) proved it.

The adjacency spectrum, however, fails before that: vertices of high degree, of order \log n / \log\log n, create eigenvalues larger than the community outlier, with eigenvectors localized on them. The non-backtracking matrix (Krzakala, Moore, Mossel, Neeman, Sly, Zdeborová and Zhang, 2013; Bordenave, Lelarge and Massoulié, 2018) and the Bethe Hessian (Saade, Krzakala and Zdeborová, 2014) repair spectral clustering down to the threshold. For networks with heterogeneous degrees, Ali and Couillet (2018) derive degree-corrected spectral methods with random matrix tools.

12. Kernel random matrices and kernel methods

Kernel methods compute an n \times n matrix K_{ij} = f(\cdot) of pairwise similarities and use its spectrum, for kernel PCA, spectral clustering or kernel ridge regression. In high dimension, all pairwise distances concentrate, and the kernel matrix becomes, to first order, a linear kernel plus a low-rank part and a multiple of the identity.

Linearization of kernel matrices

Let x_1, \dots, x_n \in \mathbb{R}^p be independent \mathcal{N}(0, I_p) vectors, X = [x_1, \dots, x_n], and n/p \to \gamma \in (0, \infty). The key observation is the concentration of the inner products: \frac{x_i^\top x_j}{p} = O\big(p^{-1/2}\big) \ (i \neq j), \qquad \frac{\|x_i\|^2}{p} = 1 + O\big(p^{-1/2}\big) . Every off-diagonal entry of the kernel matrix is therefore evaluated near 0, and a Taylor expansion of f there is accurate.

So the spectrum of K is, up to o(1), one large eigenvalue of order n f(0) along \mathbf{1}, plus a Marchenko–Pastur bulk scaled by f'(0) and shifted by f(1) - f(0) - f'(0).

A nonlinear kernel, in this regime, sees the data almost like a linear one. The same analysis applies to distance kernels f(\|x_i - x_j\|^2 / p), expanded around 2 (El Karoui, 2010).

Kernel spectral clustering

Consider a two-class mixture x_i = j_i\, \mu + z_i, with labels j_i \in \{-1, +1\} (half of each), z_i \sim \mathcal{N}(0, I_p) and \|\mu\|^2 = \ell of order 1. A single coordinate cannot separate the classes; only the aggregation over p coordinates can.

The phase transition of linear PCA survives the kernel: with f'(0) > 0, kernel spectral clustering on this mixture is not better than spectral clustering on X^\top X. Couillet and Benaych-Georges (2016) analyse general mixtures, where classes also differ by their covariances: the information then sits in f' and f'' at the concentration point, and the choice of the kernel controls which statistics of the classes are used.

The code checks both statements: the linearization error decreases like p^{-1/2}, and the clustering accuracy follows Proposition 12.2.

import numpy as np
from scipy.stats import norm

rng = np.random.default_rng(13)

# (a) K = f(X^T X / p) with f = exp is close to its linearization (El Karoui, 2010)
for p in [100, 200, 400, 800]:
    n = 2 * p
    X = rng.standard_normal((p, n))                    # columns x_i ~ N(0, I_p)
    G = X.T @ X / p
    K = np.exp(G)                                      # f(0) = f'(0) = 1, f(1) = e
    K_lin = (1 + 1 / (2 * p)) * np.ones((n, n)) + G + (np.e - 2) * np.eye(n)
    gap = np.abs(np.linalg.eigvalsh(K - K_lin)).max()
    print(f"p={p:4d}, n={n:4d}: ||K - K_lin|| = {gap:.3f}, ||K|| = {np.linalg.eigvalsh(K)[-1]:.0f}")

# (b) kernel spectral clustering of x_i = j_i mu + z_i, with the same kernel
p, n = 400, 800
c = p / n
j = np.repeat([1.0, -1.0], n // 2)
H = np.eye(n) - np.ones((n, n)) / n                    # centring removes the f(0) 1 1^T part
print(" ||mu||^2  alignment  theory  error rate  theory")
for ell in [0.5, 1.0, 1.5, 3.0]:
    mu = np.zeros(p); mu[0] = np.sqrt(ell)
    Xm = np.outer(mu, j) + rng.standard_normal((p, n))
    v = np.linalg.eigh(H @ np.exp(Xm.T @ Xm / p) @ H)[1][:, -1]
    align = (v @ j) ** 2 / n
    err = min(np.mean(np.sign(v) != j), np.mean(np.sign(v) != -j))
    th = (1 - c / ell**2) / (1 + 1 / ell) if ell > np.sqrt(c) else 0.0
    th_err = norm.cdf(-np.sqrt(th / (1 - th)))
    print(f"{ell:8.1f}  {align:9.3f}  {th:6.3f}  {err:10.3f}  {th_err:6.3f}")

The linearization error goes from 0.98 to 0.41 as p goes from 100 to 800, against a matrix norm of 200 to 1600. With c = 0.5, the threshold is \ell = 0.71: at \ell = 0.5 the error is 0.43, close to chance, and at \ell = 3 it is 0.059 against 0.060 predicted by the Gaussian heuristic of ch. 11.

13. Complements for self-study

This chapter collects results that the course states or uses without full proof: the general-entries proof of the Marchenko–Pastur theorem, the moment method, the free probability view of the structured models, the fluctuations of linear statistics, and the shrinkage of covariance estimators. It is not taught in class; read it after ch. 3 and ch. 4.

13.1 The Marchenko–Pastur theorem for general entries

Theorem 3.1 holds for independent and identically distributed entries with mean 0 and variance 1, without any higher moment. The proof of Bai and Silverstein (2010) has three steps.

1. Truncation and centring. Replace X_{ij} by X_{ij} \mathbf{1}_{|X_{ij}| \leq \eta_n \sqrt{n}}, recentred and rescaled, with \eta_n \to 0 slowly. Proposition 1.6 and a bound on the Lévy distance by the Frobenius norm of the difference show that the ESD does not change in the limit.
2. Leave-one-out. With Q_{-j} the resolvent without the column x_j, Sherman–Morrison gives 1 + z\, m_N = \frac{1}{N} \mathrm{tr}(Q \Sigma) = \frac{1}{N} \sum_{j=1}^n \frac{\frac{1}{n} x_j^\top Q_{-j} x_j}{1 + \frac{1}{n} x_j^\top Q_{-j} x_j} .
3. Concentration. The trace lemma (Proposition 0.5, in its fourth-moment form after truncation) gives, uniformly in j, \frac{1}{n} x_j^\top Q_{-j} x_j \approx \frac{1}{n} \mathrm{tr}\, Q_{-j} \approx c_n m_N . The denominators stay away from 0 because |1 + q| \geq v / |z| when q is such a quadratic form.

Inserting the approximation gives 1 + z m_N \approx \frac{m_N}{1 + c_n m_N}, which rearranges into the equation of Theorem 3.1. The almost sure convergence of m_N - \mathbb{E}\, m_N uses a martingale decomposition along the columns instead of the Poincaré–Nash inequality.

13.2 The moment method

For the GOE of Theorem 2.5, expand the moments as sums over closed walks: \mathbb{E}\, \frac{1}{N} \mathrm{tr}\, W^{k} = \frac{1}{N} \sum_{i_1, \dots, i_k} \mathbb{E}\big[W_{i_1 i_2} W_{i_2 i_3} \cdots W_{i_k i_1}\big] .

• Only walks that use each edge at least twice contribute. Those that contribute at order 1 use each edge exactly twice and trace a tree with k/2 edges.
• Such trees are counted by the Catalan number C_{k/2}, which is the number of non-crossing pair partitions of k points. Odd moments vanish.
• So the even moments tend to 1, 2, 5, 14, \dots, the moments of the semicircle law, which is determined by its moments because it has compact support.

For \Sigma = \frac{1}{n} X X^\top, the walks alternate between row and column indices, and the count weights each non-crossing partition by a power of c. This gives the Narayana polynomials of Proposition 3.2.

The same expansion with k of order \log N controls the largest eigenvalue (Füredi and Komlós, 1981, for Wigner matrices; Yin, Bai and Krishnaiah, 1988, for sample covariance matrices). The fourth moment is necessary: Bai, Silverstein and Yin (1988) show that if \mathbb{E}\, X_{11}^4 = \infty, then \limsup \lambda_{\max}(\Sigma) = \infty almost surely, while the ESD still converges to \mu_c.

13.3 Free probability

Two sequences of matrices are asymptotically free when, roughly, their eigenvectors are in generic position with respect to each other, for instance B_N = U B U^* with U Haar-distributed (Voiculescu, 1991). Then the limit spectra of sums and products depend only on the limit spectra of the summands.

• Sums. Let G(z) = -m(z) and define the R-transform by G\big(R(w) + 1/w\big) = w. The R-transforms of free matrices add. The semicircle law of variance s has R(w) = s\, w, so the sum of two free semicircles of variance 1 is a semicircle of variance 2.
• Marchenko–Pastur. Substituting z = R(G) + 1/G into the equation of Theorem 3.1 gives R_{\mu_c}(w) = \frac{1}{1 - c\, w} , the R-transform of a free Poisson law with rate 1/c and jump size c.
• Products. The S-transforms of free positive matrices multiply. The limit ESD of \hat C = C^{1/2} \Sigma C^{1/2} is the free multiplicative convolution H \boxtimes \mu_c, and Silverstein’s equation (Proposition 4.2) is its explicit form.

Mingo and Speicher (2017) give a complete introduction; Potters and Bouchaud (2020) use the R- and S-transforms throughout for applications.

13.4 Fluctuations of linear spectral statistics

A linear spectral statistic is \sum_k f(\lambda_k). The mutual information of ch. 5 (f = \log(1 + \rho\, \cdot)) and \log \det \hat C (f = \log) are examples. Proposition 2.4 already suggests that such sums fluctuate by O(1), without the \sqrt{N} of independent sums.

The centring uses \mu_{c_n}, not \mu_c, and for f = \log it has a closed form: \int \log x\; \mu_c(dx) = \frac{c - 1}{c}\, \log(1 - c) - 1 \qquad (c < 1) .

import numpy as np

rng = np.random.default_rng(14)
print("  N    n   mean (theory)       variance (theory)")
for N, n in [(50, 100), (200, 400)]:
    c = N / n
    centre = N * ((c - 1) / c * np.log(1 - c) - 1)          # N * int log x dMP_c(x)
    stats = []
    for _ in range(2000):
        X = rng.standard_normal((N, n))
        stats.append(np.linalg.slogdet(X @ X.T / n)[1] - centre)
    stats = np.array(stats)
    print(f"{N:4d} {n:4d}   {stats.mean():6.3f} ({0.5 * np.log(1 - c):6.3f})"
          f"     {stats.var():6.3f} ({-2 * np.log(1 - c):6.3f})")

For N = 200 and n = 400 the mean is -0.354 against -0.347 and the variance 1.38 against 1.39.

For MIMO channels, Hachem, Loubaton and Najim (2008) proved the corresponding CLT for the mutual information with a variance profile. It gives the outage probability \mathbb{P}(I_N < R) at a target rate R. With the numbers of ch. 5 (N = n = 8, 10 dB, mean 21.8 bits, standard deviation about 1.27 bits), the Gaussian approximation gives an outage probability of \Phi(-1.41) \approx 8\% at 20 bits per channel use.

13.5 Shrinkage of the sample covariance

Since the sample eigenvalues are biased, a better covariance estimator keeps the sample eigenvectors and changes the eigenvalues. Ledoit and Wolf (2004) shrink linearly towards a multiple of the identity, with an intensity estimated consistently when N/n \to c.

The best rotation-equivariant estimator \sum_k d_k\, \hat u_k \hat u_k^\top in Frobenius norm uses the oracle values d_k = \hat u_k^\top C\, \hat u_k. Ledoit and Péché (2011) show that d_k is asymptotically a function of \hat\lambda_k alone: d(\lambda) = \frac{\lambda}{\big|1 - c - c\, \lambda\, \breve m(\lambda)\big|^2}, \qquad \breve m(\lambda) = \lim_{\eta \downarrow 0} m_F(\lambda + i\eta) . Ledoit and Wolf (2012, 2020) estimate this nonlinear shrinkage from the data. The code compares the oracle with the formula, computing m_F by the deterministic equivalent of ch. 4:

import numpy as np

rng = np.random.default_rng(1)
N, n = 400, 800
c = N / n
pop = np.repeat([1.0, 4.0], N // 2)                  # population eigenvalues
Y = np.sqrt(pop)[:, None] * rng.standard_normal((N, n))
lam, U = np.linalg.eigh(Y @ Y.T / n)
oracle = np.sum(U**2 * pop[:, None], axis=0)         # u_k^T C u_k

def m_F(z, iters=3000):                              # deterministic equivalent of ch. 4
    d = 1.0 + 0j
    for _ in range(iters):
        d = 0.5 * d + 0.5 * np.sum(pop / (pop / (1 + d) - z)) / n
    return np.mean(1 / (pop / (1 + d) - z))

for k in [20, 100, 200, 300, 380]:
    m = m_F(lam[k] + 1e-4j)
    d = lam[k] / abs(1 - c - c * lam[k] * m) ** 2
    print(f"sample eigenvalue {lam[k]:6.3f}: oracle {oracle[k]:.3f}, Ledoit-Peche formula {d:.3f}")

A sample eigenvalue of 0.22 is pulled up to about 1.31, and one of 7.32 pulled down to 3.85: the shrunk values stay between the population eigenvalues 1 and 4. The formula follows the oracle within 10%, and within 0.2% at the two extremes.

13.6 Universality by replacement

Many results of the course were proved for Gaussian entries. The Lindeberg replacement strategy transfers them: swap the entries one at a time for Gaussian ones and control the change of a smooth function of the matrix by a third-order Taylor expansion (Chatterjee, 2006). Matching the first two moments gives the ESD; matching four moments gives local statistics, which is the four moment theorem of Tao and Vu (2011).

14. Research directions and open problems

This chapter is self-study. It points to active research that uses the tools of the course, in a neutral voice, with references to published work.

14.1 High-dimensional regression and double descent

For least squares with n samples, p isotropic Gaussian features, noise variance \sigma^2 and p/n = \gamma < 1, the excess prediction risk tends to \sigma^2\, \frac{\gamma}{1 - \gamma} , a direct consequence of \mathbb{E}\, \mathrm{tr}(X^\top X)^{-1}, as in Proposition 9.2. It explodes at \gamma = 1 and, for the minimum-norm interpolator beyond it, decreases again: the double descent curve. Hastie, Montanari, Rosset and Tibshirani (2022) derive these curves for ridge and ridgeless regression with general covariances, using deterministic equivalents. Mei and Montanari (2022) do it for random features regression, where a nonlinear random matrix appears.

14.2 Neural networks

Pennington and Worah (2017) and Louart, Liao and Couillet (2018) characterize the spectrum of the Gram matrix of random features \sigma(W X), which ch. 12 anticipates: a nonlinearity acts through a few Gaussian moments of its derivatives. Fan and Wang (2020) extend this to the conjugate and neural tangent kernels of networks of linear width.

Pennington and Bahri (2017) model the Hessian of the loss with free probability. Martin and Mahoney (2021) observe heavy-tailed empirical spectra in the weight matrices of trained networks, and propose them as a diagnostic of training quality; a theory that explains these spectra is open.

14.3 Beyond Gaussian and independent data

Real data are neither Gaussian nor made of independent entries. Concentration of measure gives a framework: data that are Lipschitz images of Gaussian vectors, such as the outputs of generative networks, behave like Gaussian mixtures for the spectral quantities of the course (Seddik, Louart, Tamaazousti and Couillet, 2020). Couillet and Liao (2022) develop this program for machine learning.

Heavy tails change the picture: the limit ESD is no longer Marchenko–Pastur when the variance is infinite (Ben Arous and Guionnet, 2008). Dependent samples, such as time series, require new deterministic equivalents.

14.4 Detection and estimation below the spectral threshold

Below the BBP threshold of ch. 6, no eigenvalue separates from the bulk, but the signal is not invisible to all tests. Onatski, Moreira and Hallin (2013) show that likelihood ratio tests based on all eigenvalues keep a non-trivial power below the threshold.

Perry, Wein, Bandeira and Moitra (2018) and El Alaoui, Krzakala and Jordan (2020) determine when detection is possible at all in spiked models with structured priors. The gap between what is possible and what efficient algorithms achieve is an open question in many such models.

14.5 Random tensors

Data with more than two indices lead to spiked tensor models. Montanari and Richard (2014) introduced tensor PCA and found a gap between the statistical threshold and the one reached by known polynomial-time algorithms. Ben Arous, Mei, Montanari and Nica (2019) study the landscape of its likelihood. Goulart, Couillet and Comon (2022) analyse random tensors through random matrices built from their contractions.

14.6 Communications and positioning

Massive MIMO, with many more base-station antennas than users (Marzetta, 2010), is the regime where the deterministic equivalents of ch. 5 are most accurate. Hoydis, ten Brink and Debbah (2013) use them to size such systems.

The positioning problem of ch. 9 with many sources, unknown channel models and few samples remains a field of applications for G-estimation. Estimating the whole population spectrum from the sample one is a related inverse problem (El Karoui, 2008; Ledoit and Wolf, 2020).

Open problems

• Non-asymptotic versions of the deterministic equivalents, with explicit constants, for moderate N and n.
• Spectral methods with optimal thresholds for sparse, heterogeneous or dynamic graphs.
• The choice of a kernel, or a nonlinearity, that extracts the most information from structured data in the proportional regime.
• A predictive theory of the spectra of trained, not random, neural network weights.

Exercises

Levels: ★ apply the course, ★★ combine results or prove a short statement, ★★★ open-ended or research-flavoured. Exercises 7, 15, 16 and 17 are computational; the others can be checked with the code of the chapters.

1. ★ Stieltjes transforms by hand. (a) Compute the Stieltjes transform of the uniform law on [0, 1] and recover its density at x = 1/2 with Theorem 1.4. (b) For A = \mathrm{diag}(0, 2), compute m_{\mu_A}(1 + i) and check properties 1–3 of Proposition 1.3.

   Solution: (a) With the principal branch of the logarithm, m(z) = \int_0^1 \frac{d\lambda}{\lambda - z} = \log \frac{1 - z}{-z} . Let z = x + i\varepsilon with 0 < x < 1. As \varepsilon \downarrow 0, the argument of 1 - z tends to 0 and that of -z to -\pi.

   So \Im\, m tends to \pi, and the density is 1.

   (b) The two atoms give m = \frac{1}{2}\Big( \frac{1}{-1 - i} + \frac{1}{1 - i} \Big) = \frac{1}{2}\Big( \frac{-1 + i}{2} + \frac{1 + i}{2} \Big) = \frac{i}{2} . Its imaginary part is positive, |m| = 0.5 \leq 1/\Im z = 1, and -iy\, m(iy) \to 1 by direct expansion.

2. ★ Centring the data does not change the limit. Let \bar x = \frac{1}{n} X \mathbf{1} be the sample mean and \Sigma_c the sample covariance of the centred vectors x_j - \bar x. Show that \Sigma - \Sigma_c has rank one, and deduce that \Sigma_c has the same Marchenko–Pastur limit and the same limit of its largest eigenvalue as \Sigma.

   Solution: expanding the sum gives \Sigma_c = \Sigma - \bar x \bar x^\top, a rank-one positive semidefinite correction. By Proposition 1.6, \sup_t |F^{\Sigma} - F^{\Sigma_c}| \leq 1/N, so the limits coincide (Theorem 1.5).

   By interlacing (Proposition 0.2), \lambda_1(\Sigma_c) \geq \lambda_2(\Sigma). Since \Sigma_c \preceq \Sigma, Theorem 0.1 also gives \lambda_1(\Sigma_c) \leq \lambda_1(\Sigma).

   Both bounds tend to (1 + \sqrt{c})^2: the first because the second eigenvalue also sticks to the edge (Theorem 3.4 and the ESD near the edge).

3. ★ Moments of the semicircle law. (a) Compute \int x^2 \sigma(dx) and \int x^4 \sigma(dx) from the density. (b) For the GOE of Theorem 2.5, compute \mathbb{E}\, \frac{1}{N} \mathrm{tr}\, W^2 exactly.

   Solution: (a) With x = 2 \sin t, \int x^{2k} \sigma(dx) = \frac{2^{2k+2}}{2\pi} \int_{-\pi/2}^{\pi/2} \sin^{2k} t \cos^2 t\, dt . Since \int \sin^2 \cos^2 = \pi/8 and \int \sin^4 \cos^2 = \pi/16, the moments are 1 and 2, the Catalan numbers C_1 and C_2.

   (b) The N(N - 1) off-diagonal entries have variance 1/N and the N diagonal ones 2/N, so \mathbb{E}\, \tfrac{1}{N} \mathrm{tr}\, W^2 = \frac{1}{N}\Big( N - 1 + 2 \Big) = 1 + \frac{1}{N} .

4. ★ Marchenko–Pastur in numbers. (a) For c = 0.25, give the support, the density at x = 1 and the asymptotic condition number. (b) With N = 100, how many samples are needed for an asymptotic condition number at most 4? (c) For N = 400 and n = 100, describe the spectrum.

   Solution: (a) The support is [0.25, 2.25], the density at 1 is \sqrt{1.25 \times 0.75} / (2\pi \times 0.25) = 0.616, the condition number 9.

   (b) With s = \sqrt{c}, we need (1 + s)/(1 - s) \leq 2, so s \leq 1/3, c \leq 1/9 and n \geq 900.

   (c) c = 4: a fraction 1 - 1/c = 0.75 of the eigenvalues is 0, and the others fill [(1 - 2)^2, (1 + 2)^2] = [1, 9].

5. ★ The capacity of a large MIMO link. (a) For c = 1 and \rho = 10, solve the fixed-point equation of Theorem 5.3 in closed form and compute \mathcal{V}_1(10). (b) Show that \mathcal{V}_1(\rho) = \log \rho - 1 + o(1) as \rho \to \infty, and interpret.

   Solution: (a) With c = 1 the equation reads \delta (11 + \delta) = 10 (1 + \delta), so \delta^2 + \delta - 10 = 0 and \delta = (\sqrt{41} - 1)/2 = 2.7016. Then 1 + \rho/(1 + \delta) = 1 + \delta = 3.7016, and \mathcal{V}_1(10) = 2 \log 3.7016 - \frac{2.7016}{3.7016} = 1.888 \text{ nats} = 2.72 \text{ bits} .

   (b) In general \delta^2 + \delta - \rho = 0, so \delta \sim \sqrt{\rho} and \mathcal{V}_1(\rho) = 2 \log(1 + \delta) - \frac{\delta}{1 + \delta} = \log \rho - 1 + o(1) . The slope 1 per receive antenna is the multiplexing gain; the constant -1 is \int \log x\, \mu_1(dx).

6. ★ Spikes in numbers. Let c = 0.5. (a) Where is the outlier of a spike \ell = 2, and what is the eigenvector alignment? (b) The top sample eigenvalue is 4.5 and the noise level is 1: estimate \ell. (c) With N = 100, how many samples are needed for a spike \ell = 0.5 to produce an outlier?

   Solution: (a) (1 + 2)(1 + 0.25) = 3.75 and \zeta = (1 - 0.125)/1.25 = 0.70. (b) \hat\lambda - 1 - c = 3, so \hat\ell = \frac{3 + \sqrt{9 - 2}}{2} = 2.82 . (c) We need \sqrt{N/n} < 0.5, so n > 400.

7. ★ A Tracy–Widom threshold, in Python. Adapt the code of ch. 7 to compute the 99% quantile of F_1, deduce the 1% threshold for \lambda_{\max}(\Sigma) with N = 100 and n = 400, and check the false-alarm rate by simulation.

   Hints and expected result: the quantile is about 2.02, so the threshold is (898.5 + 2.02 \times 15.93)/400 = 2.327. Over 2000 simulations the false-alarm rate is about 0.8%, slightly below 1%, as at the 5% level in ch. 7.

8. ★ Detectability of communities. A graph has n = 10\,000 vertices, \bar p = 0.05, p_{\mathrm{in}} = 0.055 and p_{\mathrm{out}} = 0.045. Is spectral clustering informative? Predict the overlap and the error rate.

   Solution: the threshold difference is 2\sqrt{0.05 \times 0.95 / 10\,000} = 0.00436, so \theta = 0.01 / 0.00436 = 2.29 > 1. The predicted squared overlap is 1 - 1/\theta^2 = 0.81 and the error rate \Phi\big(-\sqrt{\theta^2 - 1}\big) = \Phi(-2.06) \approx 0.019 .

9. ★★ The Poincaré–Nash bound for the GOE. With W = (G + G^\top)/\sqrt{2N} and m_N = \frac{1}{N}\mathrm{tr}(W - z)^{-1}, show that \mathbb{E}\, |m_N - \mathbb{E}\, m_N|^2 \leq \frac{2}{N^2 v^4} .

   Solution: as a function of the N^2 independent entries of G, using the symmetry of Q, \frac{\partial m_N}{\partial G_{ab}} = -\frac{1}{N\sqrt{2N}}\, \mathrm{tr}\big(Q^2 (E_{ab} + E_{ba})\big) = -\frac{2}{N\sqrt{2N}}\, (Q^2)_{ab} . The sum of the squared moduli is \frac{2}{N^3} \|Q^2\|_F^2, and \|Q^2\|_F^2 \leq N / v^4. Theorem 2.2 concludes.

10. ★★ The second moment of a sample covariance matrix. For independent entries with mean 0, variance 1 and fourth moment \kappa, show that \mathbb{E}\, \frac{1}{N} \mathrm{tr}\, \Sigma^2 = 1 + c_n + \frac{\kappa - 2}{n} , and compare with Proposition 3.2.

    Solution: expand the trace as \mathrm{tr}\, \Sigma^2 = \frac{1}{n^2} \sum_{i, k} \Big(\sum_j X_{ij} X_{kj}\Big)^2 . For i \neq k the expectation of the square is n; for i = k it is n\kappa + n(n - 1). Summing, \mathbb{E}\, \mathrm{tr}\, \Sigma^2 = \frac{1}{n^2}\Big(N(N-1)\, n + N\big(n^2 + (\kappa - 1)\, n\big)\Big) . Dividing by N gives the formula; its limit 1 + c is the second Narayana moment. The fourth moment only enters at order 1/n.

11. ★★ Variance profiles that factorize. (a) Show that when \sigma_{ij}^2 = d_i depends on the row only, the system of Theorem 4.5 reduces to Theorem 4.1 with C = \mathrm{diag}(d_i). (b) What happens when \sigma_{ij}^2 = \tilde d_j depends on the column only?

    Solution: (a) The right-hand side of the equation for \tilde t_j no longer depends on j. So \tilde t_j = \tilde t = -1/(z(1 + \delta)) with \delta = \frac{1}{n}\sum_i d_i t_i. Then t_i = \frac{-1}{z(1 + d_i \tilde t)} = \frac{1}{\frac{d_i}{1 + \delta} - z} , the diagonal of \bar Q, and \delta = \frac{1}{n} \mathrm{tr}\, C \bar Q.

    (b) By symmetry it is the population model for the n \times n matrix Y^* Y, with \mathrm{diag}(\tilde d_j) as population covariance: column variances are a covariance between samples.

12. ★★ An outlier below the bulk. Let c < 1 and C = I + \ell u u^\top with \ell \in (-1, 0). Adapt the proof of Theorem 6.4 to show that an outlier appears below (1 - \sqrt{c})^2 iff \ell < -\sqrt{c}, at the value (1 + \ell)(1 + c/\ell).

    Check it for c = 0.1 and \ell = -0.9.

    Solution: steps 1–4 of the proof are unchanged and give \lambda = (1 + \ell)(1 + c/\ell), now with \lambda to the left of the bulk. The map \ell \mapsto (1 + \ell)(1 + c/\ell) equals (1 - \sqrt{c})^2 at \ell = -\sqrt{c} and is below it on (-1, -\sqrt{c}). The branch condition holds exactly there.

    For c = 0.1 and \ell = -0.9, the outlier is at 0.1 \times (1 - 1/9) = 0.089, below the edge 0.468. A simulation with N = 200 and n = 2000 gives 0.089.

    For \ell = -0.2, there is no outlier: the smallest eigenvalue stays at the edge.

13. ★★ The GLRT statistic is pivotal. (a) Show that T_N is invariant under y_t \mapsto \alpha U y_t for a non-zero scalar \alpha and a unitary U. Deduce that its law under \mathcal{H}_0 depends only on the dimensions.

    (b) For N = 50, n = 100 and \rho = 1, compare the limit of T_N under \mathcal{H}_1 with the 5% threshold of ch. 8.

    Solution: (a) \hat R \mapsto |\alpha|^2 U \hat R U^* multiplies the largest eigenvalue and the normalized trace by |\alpha|^2. Under \mathcal{H}_0, y_t = \sigma w_t, so the law of T_N is that of the case \sigma = 1.

    (b) The limit is (1 + 1)(1 + 0.5) = 3.0, essentially the finite-size threshold 3.001 at N = 50. At this size the test is at the edge of its range: T_N has mean about 2.94, still below its limit, and a simulation as in ch. 8 gives a power of about 0.33.

    Asymptotically the picture is better. With c = 0.5 fixed, the threshold decreases to the edge 2.914 while the limit under \mathcal{H}_1 stays at 3.0, so the power tends to 1: \rho = 1 > \sqrt{c} is detectable. The simulated power is about 0.58 at N = 200 and 0.9 at N = 800.

14. ★★ Two G-estimators. (a) Using the Wishart moments of Proposition 9.2, show that the corrected estimator of \frac{1}{N} \mathrm{tr}\, C^2 has a relative bias \frac{1}{n} - \frac{2}{n^2} for real Gaussian data. (b) For C = I, N = 100 and n = 200, compute the expectations of the naive and corrected estimators of \frac{1}{N}\mathrm{tr}\, C^{-1}.

    Solution: (a) Let t_1 = \frac{1}{N}\mathrm{tr}\, C and t_2 = \frac{1}{N}\mathrm{tr}\, C^2. By Proposition 9.2, \mathbb{E}\, \frac{1}{N} \mathrm{tr}\, \hat C^2 = \frac{n + 1}{n}\, t_2 + c_n t_1^2 . The trace \mathrm{tr}\, \hat C has mean \mathrm{tr}\, C and variance \frac{2}{n} \mathrm{tr}\, C^2, so \mathbb{E}\, c_n \Big(\frac{1}{N}\mathrm{tr}\, \hat C\Big)^2 = c_n t_1^2 + \frac{2}{n^2}\, t_2 . The difference is t_2 (1 + \frac{1}{n} - \frac{2}{n^2}).

    (b) The naive estimator has mean n/(n - N - 1) = 200/99 = 2.02; the corrected one, \frac{99}{200} \times 2.02 = 1.

15. ★★★ The resolution of G-MUSIC, in Python. With the code of ch. 10 (N = 40, two sources of power 3), estimate the probability of resolution of MUSIC and G-MUSIC for separations from 1° to 6° and n \in \{40, 80, 160\}. Relate the results to the spikes \ell_{1,2} = P(1 \pm |a(\theta_1)^* a(\theta_2)|).

    Hints and expected result: for 10° and 12°, |a(\theta_1)^* a(\theta_2)| = 0.39, so the spikes are 4.17 and 1.83, both above \sqrt{c} = 0.71 at n = 80.

    Their alignments are \zeta = 0.87 and 0.67. This imbalance is what merges the two minima of MUSIC (15% resolution against 83% for G-MUSIC).

    As the separation shrinks, the smaller spike approaches \sqrt{c} and both methods fail; larger n lowers \sqrt{c} and moves the failure to smaller separations. A complete answer plots both probabilities against the separation and marks where the smaller spike crosses \sqrt{c}.

16. ★★★ Kernel clustering when the classes differ by their covariance. Let x_i \sim \mathcal{N}(0, (1 + j_i a / \sqrt{p}) I_p) with p = 400, n = 800 and a \in \{1, 2, 3, 4\}.

    Compare spectral clustering (top eigenvector of HKH) with three kernels: the linear kernel X^\top X / p, the Gaussian kernel e^{-t/2} and the kernel (t - 2)^2, where t = \|x_i - x_j\|^2 / p. Explain the ranking with ch. 12.

    Hints and expected result: the linear kernel fails (error about 0.49): the classes have the same mean. The class information is in \psi_i = \|x_i\|^2 / p, whose means differ by 2a/\sqrt{p} while each fluctuates by \sqrt{2/p}. A rule based on \psi_i alone has error about \Phi(-a/\sqrt{2}): 0.24, 0.08, 0.02, 0.002.

    Expand f(\psi_i + \psi_j - 2 x_i^\top x_j / p) around 2 and write \tilde\psi = \psi - \mathbf{1}. The first-order terms in \psi are killed by the centring.

    At second order, the cross term 2 \tilde\psi_i \tilde\psi_j of the square brings the information in as f''(2)\, \tilde\psi \tilde\psi^\top. The bulk -2 f'(2)\, X^\top X / p hides it when f'(2) \neq 0.

    So the kernel with f'(2) = 0 is near the \psi-rule (errors about 0.25, 0.08, 0.02, 0.003). The Gaussian kernel is much worse while the signal is weak (about 0.48, 0.45, 0.07, 0.01): it stays close to chance up to a = 2. This is the mechanism analysed by Couillet and Benaych-Georges (2016).

17. ★★★ Outage of a MIMO channel. For N = n \in \{2, 4, 8\} at 10 dB, estimate by Monte Carlo the mean and standard deviation of I_N, and compare the Gaussian approximation of ch. 13 of the outage probability \mathbb{P}(I_N < R) with the empirical one, at the rates where the approximation predicts 5% and 1%. Discuss the tails.

    Hints and expected result: the means are about 5.55, 10.96 and 21.8 bits and the standard deviations about 1.32, 1.28 and 1.27 bits: the fluctuations do not grow with N. With 20 000 draws, the empirical outage at the predicted 5% point is about 4.6%, 4.8% and 4.8%, and at the predicted 1% point about 0.6%, 0.9% and 0.9%, up to a Monte Carlo standard error of about 0.1 percentage point.

    The Gaussian approximation is already good at N = 4; at N = 2 the left tail is lighter than Gaussian. Far tails, needed for very low outage targets, call for large-deviation estimates rather than the CLT.

Projects and assessment

The official assessment, per master-mva.com, counts the lab reports for 50% and a final exam, a written test or a project, for 50%; the Université Paris-Saclay course detail classifies it as “Evaluation Continue non Intégrale” (partial continuous assessment). The four lab slots of the session plan produce the lab reports. The subjects below are sized for the project option of the final exam, and the written exam proposal shows what a written test on this page can ask. None of these proposals is official material.

Solo projects

1. [Solo · ≈ 20 h] The BBP phase transition: a paper study. Reproduce the eigenvalue and eigenvector transitions of Baik, Ben Arous and Péché (2005) and Paul (2007) by simulation, for real and complex Gaussian data.
   • Implement: the spiked covariance model, the spike-inversion estimator of ch. 6 and the alignment correction.
   • Success criterion: limits matched within 0.02 for N \geq 400 away from the threshold; a fit on four matrix sizes showing that the transition window shrinks like N^{-1/3}.
   • Deliverables: code, a 4–6 page report, a 10-minute talk. Chapters 3, 6, 7.
2. [Solo · ≈ 20 h] G-MUSIC for large arrays: a paper study. Study Mestre and Lagunas (2008) and compare MUSIC, the spiked G-MUSIC of ch. 10 and the G-MUSIC estimator of the paper on synthetic uniform linear arrays.
   • Success criterion: for two sources 2° apart, the SNR at which each method resolves them with probability 0.9, for N \in \{20, 40\} and n \in \{N/2, N, 2N\}, and the paper’s qualitative findings reproduced.
   • Deliverables: code, report, talk. Chapters 6, 9, 10.
3. [Solo · ≈ 15 h] Computing Tracy–Widom laws: a paper study. Implement the Nyström method of Bornemann (2010) for F_2 (Airy kernel) and F_1 (Ferrari–Spohn formula).
   • Success criterion: means and variances to 6 digits against the values tabulated in Bornemann’s 2010 review; a table of the empirical levels of Johnstone’s (2001) approximation for N, n \in \{5, 10, 20, 50\}.
   • Deliverables: code, report, talk. Chapters 7, 8.

Group projects

1. [Group of 2–3 · ≈ 30 h each] Spectrum sensing for cognitive radio. Reproduce the comparison of detectors of Bianchi, Debbah, Maïda and Najim (2011) on synthetic complex Gaussian data.
   • Work split: one member writes the GLRT with Tracy–Widom thresholds and a finite-size calibration; one the condition-number and energy detectors under noise uncertainty; one the estimation of the number of sources (Kritchman and Nadler, 2009).
   • Success criterion: ROC curves for N \in \{10, 50\}, the detectability limit \rho = \sqrt{c} checked within 10%, and the energy detector’s failure under a 0.5 dB noise error shown.
   • Deliverables: code, a 6-page report, a 15-minute talk. Chapters 6, 7, 8.
2. [Group of 2–3 · ≈ 30 h each] Community detection from dense to sparse graphs. Compare spectral methods on stochastic block models and on small real graphs bundled with NetworkX (karate club, Les Misérables co-occurrences).
   • Work split: the dense regime and Theorem 11.2; the sparse regime with the Bethe Hessian of Saade, Krzakala and Zdeborová (2014) and the non-backtracking matrix; the real graphs and the evaluation protocol.
   • Success criterion: overlap curves against (a - b)^2 / (2(a + b)) that locate the Kesten–Stigum threshold within 10% for the Bethe Hessian, and the failure of the adjacency spectrum in the sparse regime.
   • Deliverables: code, report, talk. Chapters 6, 11.
3. [Group of 2–3 · ≈ 30 h each] Locating transmitters with G-estimation. Extend ch. 9 to several sources per sensor, following Couillet, Silverstein, Bai and Debbah (2011).
   • Work split: the estimators (cluster means, Mestre’s estimator, the estimator of the paper); the positioning pipeline with a noisy path-loss model; a study of what happens when the separation condition fails.
   • Success criterion: position error against n for N = 32 and 64, with the G-estimators’ error decreasing while the naive one stays biased.
   • Deliverables: code, report, talk. Chapters 4, 9.

Written exam (proposal)

Format: 3 hours, one handwritten A4 sheet allowed, no electronic device; four independent exercises.

1. Course question (4 points). Define the Stieltjes transform of a probability measure. State the inversion formula and the continuity theorem. State the Marchenko–Pastur theorem, with the support and the density of the limit law.

   Expected: Definition 1.2, Theorems 1.4 and 1.5 with their hypotheses, Theorem 3.1 and Proposition 3.2, including the atom at 0 when c > 1.

2. Computation (6 points). A sensor array has N = 25 antennas and records n = 100 samples of white noise of variance 1, possibly with sources. The two largest sample eigenvalues are 3.0 and 2.0.

   (a) Give the support of the noise eigenvalues. (b) Which sample eigenvalues are outliers? (c) Estimate the corresponding spike and the alignment of its sample eigenvector. (d) How many sources are detected?

   Expected: (a) c = 0.25 and the support is [0.25, 2.25]. (b) Only 3.0 is above the edge. (c) With \hat\lambda - 1 - c = 1.75, \hat\ell = \frac{1.75 + \sqrt{1.75^2 - 1}}{2} = 1.59, \qquad \zeta = \frac{1 - 0.25/1.59^2}{1 + 0.25/1.59} = 0.78 . (d) One source; a second source with spike below \sqrt{c} = 0.5 would be undetectable.

3. Short proof (5 points). Let W be a GOE matrix normalized as in Theorem 2.5, \theta > 1 and v a unit vector. Assuming Proposition 6.2, prove that the largest eigenvalue of W + \theta v v^\top tends to \theta + 1/\theta and that |\langle v, \hat v\rangle|^2 \to 1 - 1/\theta^2.

   Expected: the determinant lemma, the equation \theta\, v^\top R(\lambda) v = 1, the limit -m_{\mathrm{sc}}(\lambda) = 1/\theta, the quadratic m^2 + \lambda m + 1 = 0, then the eigenvector identity |\langle v, \hat v\rangle|^2 = 1/(\theta^2 v^\top R^2 v) and the derivative of m_{\mathrm{sc}}, as in Theorem 6.3.

4. Problem (5 points). (a) From the Marchenko–Pastur equation, compute \int x^{-1}\, \mu_c(dx) for c < 1. (b) Let \hat C = C^{1/2} \Sigma C^{1/2} with Gaussian data. Show that \frac{1}{N}\mathrm{tr}\, \hat C^{-1} - \frac{1}{1 - c}\, \frac{1}{N}\mathrm{tr}\, C^{-1} \to 0 , and propose a G-estimator of \frac{1}{N}\mathrm{tr}\, C^{-1}. (c) Explain the consequence for a plug-in Mahalanobis distance.

   Expected: (a) at z = 0 the equation gives (c - 1)\, m + 1 = 0, so \int x^{-1} d\mu_c = m_c(0) = \frac{1}{1 - c}.

   (b) By cyclicity, \hat C^{-1} and C^{-1} \Sigma^{-1} have the same trace. By orthogonal invariance, \mathbb{E}\, \Sigma^{-1} is a multiple of I, equal to \frac{n}{n - N - 1} I, and the variance vanishes. The estimator is (1 - c_n) \frac{1}{N}\mathrm{tr}\, \hat C^{-1}.

   (c) Plug-in Mahalanobis distances are inflated by the factor 1/(1 - c) on average.

Grading rubric (proposal)

References

Official references (master-mva.com)

• G. W. Anderson, A. Guionnet and O. Zeitouni. An Introduction to Random Matrices. Cambridge Studies in Advanced Mathematics 118, Cambridge University Press, 2010.
• R. Couillet and M. Debbah. Random Matrix Methods for Wireless Communications. Cambridge University Press, 2011.
• T. Tao. Topics in Random Matrix Theory. Graduate Studies in Mathematics 132, American Mathematical Society, 2012.

Textbooks

• Z. Bai and J. W. Silverstein. Spectral Analysis of Large Dimensional Random Matrices, 2nd edition. Springer, 2010.
• R. Couillet and Z. Liao. Random Matrix Methods for Machine Learning. Cambridge University Press, 2022.
• L. Erdős and H.-T. Yau. A Dynamical Approach to Random Matrix Theory. Courant Lecture Notes 28, American Mathematical Society, 2017.
• V. L. Girko. Statistical Analysis of Observations of Increasing Dimension. Kluwer Academic Publishers, 1995.
• J. A. Mingo and R. Speicher. Free Probability and Random Matrices. Fields Institute Monographs 35, Springer, 2017.
• L. Pastur and M. Shcherbina. Eigenvalue Distribution of Large Random Matrices. Mathematical Surveys and Monographs 171, American Mathematical Society, 2011.
• M. Potters and J.-P. Bouchaud. A First Course in Random Matrix Theory: for Physicists, Engineers and Data Scientists. Cambridge University Press, 2020.
• A. M. Tulino and S. Verdú. Random matrix theory and wireless communications. Foundations and Trends in Communications and Information Theory, 1(1), 2004.

Limit laws and their tools (ch. 1–4, 13)

• Z. D. Bai and J. W. Silverstein. No eigenvalues outside the support of the limiting spectral distribution of large-dimensional sample covariance matrices. Annals of Probability, 26(1), 1998.
• Z. D. Bai and J. W. Silverstein. Exact separation of eigenvalues of large dimensional sample covariance matrices. Annals of Probability, 27(3), 1999.
• Z. D. Bai and J. W. Silverstein. CLT for linear spectral statistics of large-dimensional sample covariance matrices. Annals of Probability, 32(1A), 2004.
• Z. D. Bai, J. W. Silverstein and Y. Q. Yin. A note on the largest eigenvalue of a large dimensional sample covariance matrix. Journal of Multivariate Analysis, 26(2), 1988.
• Z. D. Bai and Y. Q. Yin. Limit of the smallest eigenvalue of a large dimensional sample covariance matrix. Annals of Probability, 21(3), 1993.
• S. Chatterjee. A generalization of the Lindeberg principle. Annals of Probability, 34(6), 2006.
• R. B. Dozier and J. W. Silverstein. On the empirical distribution of eigenvalues of large dimensional information-plus-noise-type matrices. Journal of Multivariate Analysis, 98(4), 2007.
• Z. Füredi and J. Komlós. The eigenvalues of random symmetric matrices. Combinatorica, 1(3), 1981.
• S. Geman. A limit theorem for the norm of random matrices. Annals of Probability, 8(2), 1980.
• W. Hachem, P. Loubaton and J. Najim. Deterministic equivalents for certain functionals of large random matrices. Annals of Applied Probability, 17(3), 2007.
• V. A. Marčenko and L. A. Pastur. Distribution of eigenvalues for some sets of random matrices. Mathematics of the USSR-Sbornik, 1(4), 1967.
• J. W. Silverstein and Z. D. Bai. On the empirical distribution of eigenvalues of a class of large dimensional random matrices. Journal of Multivariate Analysis, 54(2), 1995.
• J. W. Silverstein and S.-I. Choi. Analysis of the limiting spectral distribution of large dimensional random matrices. Journal of Multivariate Analysis, 54(2), 1995.
• T. Tao and V. Vu. Random matrices: universality of local eigenvalue statistics. Acta Mathematica, 206(1), 2011.
• D. Voiculescu. Limit laws for random matrices and free products. Inventiones Mathematicae, 104(1), 1991.
• E. P. Wigner. Characteristic vectors of bordered matrices with infinite dimensions. Annals of Mathematics, 62(3), 1955.
• E. P. Wigner. On the distribution of the roots of certain symmetric matrices. Annals of Mathematics, 67(2), 1958.
• Y. Q. Yin, Z. D. Bai and P. R. Krishnaiah. On the limit of the largest eigenvalue of the large dimensional sample covariance matrix. Probability Theory and Related Fields, 78(4), 1988.

Spikes and the edge (ch. 6–7)

• J. Baik, G. Ben Arous and S. Péché. Phase transition of the largest eigenvalue for nonnull complex sample covariance matrices. Annals of Probability, 33(5), 2005.
• J. Baik and J. W. Silverstein. Eigenvalues of large sample covariance matrices of spiked population models. Journal of Multivariate Analysis, 97(6), 2006.
• F. Benaych-Georges and R. R. Nadakuditi. The eigenvalues and eigenvectors of finite, low rank perturbations of large random matrices. Advances in Mathematics, 227(1), 2011.
• F. Benaych-Georges and R. R. Nadakuditi. The singular values and vectors of low rank perturbations of large rectangular random matrices. Journal of Multivariate Analysis, 111, 2012.
• F. Bornemann. On the numerical evaluation of Fredholm determinants. Mathematics of Computation, 79(270), 2010.
• F. Bornemann. On the numerical evaluation of distributions in random matrix theory: a review. Markov Processes and Related Fields, 16(4), 2010.
• M. Capitaine, C. Donati-Martin and D. Féral. The largest eigenvalues of finite rank deformation of large Wigner matrices: convergence and nonuniversality of the fluctuations. Annals of Probability, 37(1), 2009.
• N. El Karoui. Tracy–Widom limit for the largest eigenvalue of a large class of complex sample covariance matrices. Annals of Probability, 35(2), 2007.
• P. L. Ferrari and H. Spohn. A determinantal formula for the GOE Tracy–Widom distribution. Journal of Physics A: Mathematical and General, 38(33), 2005.
• S. P. Hastings and J. B. McLeod. A boundary value problem associated with the second Painlevé transcendent and the Korteweg–de Vries equation. Archive for Rational Mechanics and Analysis, 73(1), 1980.
• K. Johansson. Shape fluctuations and random matrices. Communications in Mathematical Physics, 209(2), 2000.
• I. M. Johnstone. On the distribution of the largest eigenvalue in principal components analysis. Annals of Statistics, 29(2), 2001.
• D. Paul. Asymptotics of sample eigenstructure for a large dimensional spiked covariance model. Statistica Sinica, 17(4), 2007.
• S. Péché. Universality results for the largest eigenvalues of some sample covariance matrix ensembles. Probability Theory and Related Fields, 143(3–4), 2009.
• N. S. Pillai and J. Yin. Universality of covariance matrices. Annals of Applied Probability, 24(3), 2014.
• A. Soshnikov. A note on universality of the distribution of the largest eigenvalues in certain sample covariance matrices. Journal of Statistical Physics, 108(5–6), 2002.
• C. A. Tracy and H. Widom. Level-spacing distributions and the Airy kernel. Communications in Mathematical Physics, 159(1), 1994.
• C. A. Tracy and H. Widom. On orthogonal and symplectic matrix ensembles. Communications in Mathematical Physics, 177(3), 1996.

Signal processing and communications (ch. 5, 8–10)

• P. Bianchi, M. Debbah, M. Maïda and J. Najim. Performance of statistical tests for single-source detection using random matrix theory. IEEE Transactions on Information Theory, 57(4), 2011.
• L. S. Cardoso, M. Debbah, P. Bianchi and J. Najim. Cooperative spectrum sensing using random matrix theory. Proceedings of the 3rd International Symposium on Wireless Pervasive Computing (ISWPC), 2008.
• R. Couillet, J. W. Silverstein, Z. Bai and M. Debbah. Eigen-inference for energy estimation of multiple sources. IEEE Transactions on Information Theory, 57(4), 2011.
• G. J. Foschini and M. J. Gans. On limits of wireless communications in a fading environment when using multiple antennas. Wireless Personal Communications, 6(3), 1998.
• W. Hachem, O. Khorunzhiy, P. Loubaton, J. Najim and L. Pastur. A new approach for mutual information analysis of large dimensional multi-antenna channels. IEEE Transactions on Information Theory, 54(9), 2008.
• W. Hachem, P. Loubaton and J. Najim. A CLT for information-theoretic statistics of Gram random matrices with a given variance profile. Annals of Applied Probability, 18(6), 2008.
• W. Hachem, P. Loubaton, X. Mestre, J. Najim and P. Vallet. A subspace estimator for fixed rank perturbations of large random matrices. Journal of Multivariate Analysis, 114, 2013.
• J. Hoydis, S. ten Brink and M. Debbah. Massive MIMO in the UL/DL of cellular networks: how many antennas do we need? IEEE Journal on Selected Areas in Communications, 31(2), 2013.
• S. Kritchman and B. Nadler. Non-parametric detection of the number of signals: hypothesis testing and random matrix theory. IEEE Transactions on Signal Processing, 57(10), 2009.
• T. L. Marzetta. Noncooperative cellular wireless with unlimited numbers of base station antennas. IEEE Transactions on Wireless Communications, 9(11), 2010.
• X. Mestre. Improved estimation of eigenvalues and eigenvectors of covariance matrices using their sample estimates. IEEE Transactions on Information Theory, 54(11), 2008.
• X. Mestre and M. A. Lagunas. Modified subspace algorithms for DoA estimation with large arrays. IEEE Transactions on Signal Processing, 56(2), 2008.
• R. R. Nadakuditi and A. Edelman. Sample eigenvalue based detection of high-dimensional signals in white noise using relatively few samples. IEEE Transactions on Signal Processing, 56(7), 2008.
• R. O. Schmidt. Multiple emitter location and signal parameter estimation. IEEE Transactions on Antennas and Propagation, 34(3), 1986.
• R. Tandra and A. Sahai. SNR walls for signal detection. IEEE Journal of Selected Topics in Signal Processing, 2(1), 2008.
• I. E. Telatar. Capacity of multi-antenna Gaussian channels. European Transactions on Telecommunications, 10(6), 1999.
• P. Vallet, P. Loubaton and X. Mestre. Improved subspace estimation for multivariate observations of high dimension: the deterministic signals case. IEEE Transactions on Information Theory, 58(2), 2012.
• S. Verdú and S. Shamai. Spectral efficiency of CDMA with random spreading. IEEE Transactions on Information Theory, 45(2), 1999.
• Y. Zeng and Y.-C. Liang. Eigenvalue-based spectrum sensing algorithms for cognitive radio. IEEE Transactions on Communications, 57(6), 2009.

Graphs and kernels (ch. 11–12)

• H. T. Ali and R. Couillet. Improved spectral community detection in large heterogeneous networks. Journal of Machine Learning Research, 18, 2018.
• C. Bordenave, M. Lelarge and L. Massoulié. Non-backtracking spectrum of random graphs: community detection and non-regular Ramanujan graphs. Annals of Probability, 46(1), 2018.
• R. Couillet and F. Benaych-Georges. Kernel spectral clustering of large dimensional data. Electronic Journal of Statistics, 10(1), 2016.
• A. Decelle, F. Krzakala, C. Moore and L. Zdeborová. Asymptotic analysis of the stochastic block model for modular networks and its algorithmic applications. Physical Review E, 84(6), 2011.
• N. El Karoui. The spectrum of kernel random matrices. Annals of Statistics, 38(1), 2010.
• F. Krzakala, C. Moore, E. Mossel, J. Neeman, A. Sly, L. Zdeborová and P. Zhang. Spectral redemption in clustering sparse networks. Proceedings of the National Academy of Sciences, 110(52), 2013.
• L. Massoulié. Community detection thresholds and the weak Ramanujan property. Proceedings of the 46th ACM Symposium on Theory of Computing (STOC), 2014.
• E. Mossel, J. Neeman and A. Sly. Reconstruction and estimation in the planted partition model. Probability Theory and Related Fields, 162(3–4), 2015.
• E. Mossel, J. Neeman and A. Sly. A proof of the block model threshold conjecture. Combinatorica, 38(3), 2018.
• R. R. Nadakuditi and M. E. J. Newman. Graph spectra and the detectability of community structure in networks. Physical Review Letters, 108(18), 2012.
• A. Saade, F. Krzakala and L. Zdeborová. Spectral clustering of graphs with the Bethe Hessian. Advances in Neural Information Processing Systems 27, 2014.

Statistics and learning (ch. 13–14)

• G. Ben Arous and A. Guionnet. The spectrum of heavy tailed random matrices. Communications in Mathematical Physics, 278(3), 2008.
• G. Ben Arous, S. Mei, A. Montanari and M. Nica. The landscape of the spiked tensor model. Communications on Pure and Applied Mathematics, 72(11), 2019.
• A. El Alaoui, F. Krzakala and M. I. Jordan. Fundamental limits of detection in the spiked Wigner model. Annals of Statistics, 48(2), 2020.
• N. El Karoui. Spectrum estimation for large dimensional covariance matrices using random matrix theory. Annals of Statistics, 36(6), 2008.
• Z. Fan and Z. Wang. Spectra of the conjugate kernel and neural tangent kernel for linear-width neural networks. Advances in Neural Information Processing Systems 33, 2020.
• J. H. de M. Goulart, R. Couillet and P. Comon. A random matrix perspective on random tensors. Journal of Machine Learning Research, 23, 2022.
• T. Hastie, A. Montanari, S. Rosset and R. J. Tibshirani. Surprises in high-dimensional ridgeless least squares interpolation. Annals of Statistics, 50(2), 2022.
• O. Ledoit and S. Péché. Eigenvectors of some large sample covariance matrix ensembles. Probability Theory and Related Fields, 151(1–2), 2011.
• O. Ledoit and M. Wolf. A well-conditioned estimator for large-dimensional covariance matrices. Journal of Multivariate Analysis, 88(2), 2004.
• O. Ledoit and M. Wolf. Nonlinear shrinkage estimation of large-dimensional covariance matrices. Annals of Statistics, 40(2), 2012.
• O. Ledoit and M. Wolf. Analytical nonlinear shrinkage of large-dimensional covariance matrices. Annals of Statistics, 48(5), 2020.
• C. Louart, Z. Liao and R. Couillet. A random matrix approach to neural networks. Annals of Applied Probability, 28(2), 2018.
• C. H. Martin and M. W. Mahoney. Implicit self-regularization in deep neural networks: evidence from random matrix theory and implications for learning. Journal of Machine Learning Research, 22, 2021.
• S. Mei and A. Montanari. The generalization error of random features regression: precise asymptotics and the double descent curve. Communications on Pure and Applied Mathematics, 75(4), 2022.
• A. Montanari and E. Richard. A statistical model for tensor PCA. Advances in Neural Information Processing Systems 27, 2014.
• A. Onatski, M. J. Moreira and M. Hallin. Asymptotic power of sphericity tests for high-dimensional data. Annals of Statistics, 41(3), 2013.
• J. Pennington and Y. Bahri. Geometry of neural network loss surfaces via random matrix theory. Proceedings of the 34th International Conference on Machine Learning (ICML), 2017.
• J. Pennington and P. Worah. Nonlinear random matrix theory for deep learning. Advances in Neural Information Processing Systems 30, 2017.
• A. Perry, A. S. Wein, A. S. Bandeira and A. Moitra. Optimality and sub-optimality of PCA I: spiked random matrix models. Annals of Statistics, 46(5), 2018.
• M. E. A. Seddik, C. Louart, M. Tamaazousti and R. Couillet. Random matrix theory proves that deep learning representations of GAN-data behave as Gaussian mixtures. Proceedings of the 37th International Conference on Machine Learning (ICML), 2020.