Inverse Problems and Imaging: Statistical and Stochastic Approaches

Contents

• Introduction and overview
• 0. Review: Fourier analysis, singular values and Gaussian vectors
• 1. Imaging problems, waves and noise
• 2. Least squares, ill-posedness and Bayesian regularisation
• 3. Time-reversal imaging
• 4. Random matrix theory
• 5. Detection tests and what follows them
• 6. Spectral properties of stationary random processes
• 7. Imaging with ambient noise cross-correlations
• 8. Research directions and open problems
• Exercises
• Projects and assessment
• References

Introduction and overview

Audience: an M2 student with a good M1 in applied mathematics, physics or signal processing who wants to image a medium from the waves recorded by an array of sensors — an ultrasound probe, a network of seismometers, an antenna — and to know what noise allows and forbids.

Empirical first: in Code 1.10, 41 sensors record the echoes of two small reflectors with a measurement noise as strong as the signal in a typical entry of the data matrix. Yet the two largest singular values of that matrix stand clearly above the noise edge predicted by random matrix theory, and migration puts the first reflector at the right point to within a tenth of a wavelength (Code 3.12). In seismology, correlating weeks of ambient noise at two stations reveals the travel time between them without any controlled source; ch. 7 explains why.

The module studies multi-sensor imaging, as in ultrasound echography or seismic imaging, as a statistical problem. Waves carry information about a medium to an array; noise corrupts it; an imaging method must decide whether a target is present, where it is and what it is, with error rates that can be computed. Four mathematical tools organise the course: least squares and Bayesian regularisation, time reversal and the Helmholtz–Kirchhoff identity, random matrix theory, and the spectral theory of stationary processes.

It shares random matrices with Large Random Matrices and Applications to Machine Learning, hypothesis testing with Detection Theory and Industrial Applications, Bayesian inverse problems with Learning Physical Phenomena: Metamodelling and Uncertainty Quantification, and learned reconstruction with Deep Learning for Image Restoration and Synthesis.

What you will be able to do

• Formulate an imaging problem as detection, localisation, characterisation or identification, and write the corresponding statistical model.
• Compute Green’s functions of the wave equation and the multistatic response matrix of point reflectors in the Born approximation.
• Solve a linear inverse problem by least squares, diagnose its ill-posedness from the singular values, and regularise it by truncation or Tikhonov’s method.
• Read Tikhonov regularisation as Bayesian inference, and compute posterior means and covariances.
• Derive the resolution of time-reversal and migration imaging from the Helmholtz–Kirchhoff identity.
• Use the Marchenko–Pastur law and the spiked-model phase transition to predict the singular values of noisy response matrices.
• Build detection tests with a prescribed false-alarm rate, compute their power, estimate the position and strength of a detected reflector, and identify it within a dictionary of known targets.
• Compute the power spectral density of a stationary process and the statistical fluctuations of empirical correlations.
• Retrieve travel times and Green’s functions from cross-correlations of ambient noise, and image with them.
• Implement every method above in a few lines of NumPy and check it against the theory.

Official objectives and syllabus (2026–2027)

The course description on master-mva.com sets the objective of presenting several techniques of stochastic and statistical analysis that are useful to solve inverse problems and multi-sensor imaging problems, of the type of ultrasound echography or seismic imaging. It announces eight sessions of three hours and lists four themes, each with an application; the Université Paris-Saclay course detail adds no programme items. The table below follows the themes in their official order.

Beyond the syllabus, this page adds a review chapter (ch. 0), complete proofs of the central results, runnable Python for every chapter, and a self-study chapter on research directions (ch. 8).

Where this course fits in the programme

Official prerequisites: probability, analysis (Fourier) and linear algebra. The module is a second-semester Modelling course of the Santé track.

• Refreshers — Probability and Statistics (Gaussian vectors, likelihood, tests), Linear Algebra, Analysis and Optimization (SVD, least squares, Fourier) and Scientific Python, NumPy and PyTorch.
• Semester 1 courses that prepare it — Introduction to Statistical Learning (ridge regression is Tikhonov regularisation), Convex Optimization, Computational Statistics (Bayesian computation), Machine Learning for Time Series (stationary processes), Statistical Learning for Scientific Imaging and Physical Measurements: Geosciences, Medical Image Analysis.
• Related courses of the same semester — Large Random Matrices, Detection Theory and Industrial Applications, Learning Physical Phenomena, Bayesian Machine Learning, Deep Learning for Image Restoration and Synthesis, Deep Learning for Medical Imaging, Audio Signal Processing: Time–Frequency Analysis, Remote Sensing Data.
• Afterwards — the Research Internship, for instance in a wave-imaging, ultrasound or seismology laboratory.

Course outline

• Review — Fourier transform and correlation, SVD, Weyl and Eckart–Young, complex Gaussian vectors and conditioning.
• Imaging problems, waves and noise — the four imaging tasks, the wave and Helmholtz equations, Green’s functions, the Born approximation, the multistatic response matrix, noise models.
• Least squares, ill-posedness and Bayesian regularisation — pseudo-inverse and Gauss–Markov, noise amplification, spectral filters, the discrepancy principle, Gaussian posteriors.
• Time-reversal imaging — migration as the adjoint, the time-reversal experiment, the Helmholtz–Kirchhoff identity, resolution formulas, signal-to-noise ratio, DORT and MUSIC.
• Random matrix theory — the Marchenko–Pastur law by the moment method, the edge and Tracy–Widom, the spiked model and its phase transition.
• Detection tests — Neyman–Pearson, generalised likelihood ratios, thresholds and power, then localisation, characterisation and identification of what was detected.
• Stationary random processes — Wiener–Khinchin, filtering, the variance of empirical correlations, the periodogram.
• Ambient noise cross-correlations — Green’s function retrieval, directivity of the sources, travel times, passive imaging.
• Research directions and open problems — imaging in complex media, high-dimensional detection, learned reconstruction, noise-based monitoring (self-study).

Session plan

27 h in class = 9 sessions of 3 h. The volume is the 27 h of lectures of the Université Paris-Saclay programme table; master-mva.com announces eight sessions of three hours (24 h). The plan therefore uses nine sessions of three hours, and the split of the chapters into sessions is a proposal. With eight sessions, ch. 6 (session 8) is the natural candidate for guided reading before session 9.

Level 1 = foundations, 2 = core methods and results, 3 = advanced and research-level material; the level never decreases.

Self-study, outside the 27 h (about 3 h in total): ch. 8 Research directions and open problems (≈ 3 h, with Ex. 15–16 and the paper-study projects).

How to study this module

• Read ch. 0 before session 1: the SVD, Eckart–Young and complex Gaussian vectors return in every later chapter.
• Keep one geometry in mind throughout: a linear array of sensors facing a few point reflectors. Codes 1.10 and 3.12, and Exercise 12, all use it.
• After each session, run the chapter’s code and change one parameter (noise level, number of sensors, bandwidth); predict the effect before running it.
• Do the exercises of the Practice column within the week; the ★★ exercises are close to the written-exam questions.
• The Python environment is NumPy and SciPy on a CPU; every code block runs in a few seconds on a laptop.

0. Review: Fourier analysis, singular values and Gaussian vectors

Three M1 tools carry the whole module. The Fourier transform turns the wave equation into the Helmholtz equation and correlations into products; the singular value decomposition organises least squares, imaging and detection; complex Gaussian vectors model the noise on array data.

This chapter fixes the conventions and proves the facts used later. The refreshers on linear algebra and probability give more background.

Fourier transform and correlations

With this sign, \partial_t f has transform -i\omega \hat f, so a time derivative becomes a multiplication.

The correlation formula is the heart of ch. 6 and ch. 7: correlating two recordings multiplies one spectrum by the conjugate of the other, which cancels the common phase of the source.

Singular values

Weyl’s bound says that noise of norm \|E\| cannot move a singular value by more than \|E\|. Ch. 4 shows that for random noise this bound is pessimistic, and computes the exact displacement.

Complex Gaussian vectors

At a fixed frequency, sensor data are complex numbers, and their noise is modelled by circularly symmetric complex Gaussian variables.

Unitary invariance (fact 3) is why the singular values of pure noise carry no information about any direction: ch. 4 and ch. 5 rest on it.

1. Imaging problems, waves and noise

Before computing any image, one must say what is measured, how waves carry information about the medium to the sensors, and what corrupts that information. This chapter sets up the forward model used in the whole module: point reflectors probed by an array, in the Born approximation, with additive noise.

The four imaging tasks

An array of N sensors at positions x_1, \dots, x_N emits waves and records the echoes, or simply listens. The data D are noisy, and the questions asked about the medium form a hierarchy.

Each task needs more information than the previous one, hence a higher signal-to-noise ratio. A reflector can be detected and not located, or located and not characterised; ch. 5 makes these thresholds quantitative.

Waves and Green’s functions

The time-domain formula says that an impulse emitted at y reaches x after the travel time |x - y|/c_0, with amplitude decaying like the inverse distance. Green’s functions are also symmetric, \hat G(\omega, x, y) = \hat G(\omega, y, x), even in heterogeneous media: this is the reciprocity principle.

Reflectors and the multistatic response matrix

The Born approximation neglects multiple scattering. It is accurate when the reflectivity is weak and the reflectors are small compared with the wavelength, which is the regime of this module.

In the time domain, in three dimensions, the entry A_{nm}(t) for a pulse f is an echo arriving at the round-trip travel time T_{nm}: A_{nm}(t) = -\frac{\rho_1}{16\pi^2 c_0^2\, |x_n - z_1|\, |z_1 - x_m|}\; f''\big(t - T_{nm}\big), \qquad T_{nm} = \frac{|x_n - z_1| + |z_1 - x_m|}{c_0} .

Noise models

This module mostly uses the first model, which allows exact computations; the research chapter returns to clutter.

import numpy as np
from scipy.special import hankel1

rng = np.random.default_rng(0)
lam = 1.0; k = 2 * np.pi / lam                 # lengths in wavelengths
N = 41
xs = np.stack([(np.arange(N) - 20) * lam / 2, np.zeros(N)], axis=1)   # array, aperture 20 lambda

def G(x, z):
    """2D outgoing Green's function (i/4) H0^(1)(k|x - z|)."""
    return 0.25j * hankel1(0, k * np.linalg.norm(x - z, axis=-1))

refl = [(np.array([-2.0, 30.0]), 1.0), (np.array([3.0, 32.0]), 0.6)]
A = sum(rho * k**2 * np.outer(G(xs, z), G(xs, z)) for z, rho in refl)   # Born MSR matrix
sigma = 0.01
W = sigma * (rng.standard_normal((N, N)) + 1j * rng.standard_normal((N, N))) / np.sqrt(2)
s_clean = np.linalg.svd(A, compute_uv=False)
s_noisy = np.linalg.svd(A + W, compute_uv=False)
print("typical |A_nm|        :", np.round(np.median(np.abs(A)), 4))
print("clean singular values :", np.round(s_clean[:4], 4))
print("noisy singular values :", np.round(s_noisy[:4], 4))
print("noise edge 2 sigma sqrt(N):", np.round(2 * sigma * np.sqrt(N), 4))

2. Least squares, ill-posedness and Bayesian regularisation

Most imaging problems, once linearised, read y = Kx + w: unknown parameters x, a known forward operator K, noisy data y. This chapter asks when least squares can invert K, why it fails for imaging operators, and how regularisation and Bayesian priors restore stability.

Throughout, K \in \mathbb{C}^{M \times P} has the singular value decomposition K = \sum_j s_j u_j v_j^*, with non-increasing s_j. The noise w has mean zero and covariance \sigma^2 I.

We write s_j for the singular values of K to keep \sigma for the noise.

Least squares

The mean squared error of least squares follows from the covariance: \mathbb{E}\, \|\hat x - x\|^2 = \sigma^2\, \mathrm{tr}\, (K^* K)^{-1} = \sigma^2 \sum_j \frac{1}{s_j^2} .

The best unbiased estimator can thus be useless: one tiny singular value suffices to ruin it.

Ill-posedness

In infinite dimension, a compact operator, such as a smoothing integral operator, has singular values accumulating at zero, so its inverse is unbounded. After discretisation the inverse exists but its norm is huge, and the problem stays ill posed in a quantitative sense.

The data coefficients u_j^* y split into a signal part, s_j v_j^* x, which decays with j, and a noise part of constant size \sigma. Beyond the index where s_j falls below the noise, the coefficients are pure noise; Picard’s condition says that a meaningful solution needs the signal coefficients to decay faster than s_j.

Regularisation by spectral filtering

Increasing \alpha trades variance for bias; the best \alpha depends on the unknown x, so it must be chosen from the data.

The rationale is that \mathbb{E}\|w\|^2 = M\sigma^2: fitting the data more closely than the noise level means fitting the noise. Other rules, such as the L-curve and generalised cross-validation, do not need \sigma.

Bayesian analysis

The Bayesian approach treats x as random, with a prior law that encodes what is known before the measurement. Regularisation then appears as a consequence rather than a recipe. We state the real case; the circular complex case is identical with conjugate transposes.

The posterior covariance quantifies the uncertainty. Along v_j, the posterior standard deviation is close to \sigma / s_j when s_j^2 \gg \alpha, and close to \gamma when s_j^2 \ll \alpha: the data inform the first directions, the prior alone the others.

Non-Gaussian priors change the picture. A Laplace prior gives an \ell_1-penalised MAP, which favours sparse reflectivities; the posterior is then no longer Gaussian and is explored by Markov chain Monte Carlo, as in Computational Statistics. The hyperparameters \sigma^2 and \gamma^2 can be estimated by maximising the marginal likelihood of y, which is \mathcal{N}(0, K \Gamma K^\top + \sigma^2 I).

import numpy as np

rng = np.random.default_rng(1)
n = 200
t = np.linspace(0, 1, n)
K = np.exp(-(t[:, None] - t[None, :])**2 / (2 * 0.03**2))       # Gaussian blur
K /= K.sum(axis=1, keepdims=True)
x_true = (np.abs(t - 0.3) < 0.08) + 0.6 * np.exp(-((t - 0.7) / 0.05)**2)
sigma = 1e-3
y = K @ x_true + sigma * rng.standard_normal(n)

U, s, Vt = np.linalg.svd(K)
c = U.T @ y
err = lambda x: np.linalg.norm(x - x_true) / np.linalg.norm(x_true)
print(f"condition number {s[0] / s[-1]:.1e}; least squares rel. error {err(Vt.T @ (c / s)):.1e}")

def tikhonov(alpha):
    return Vt.T @ (s / (s**2 + alpha) * c)

# Morozov discrepancy principle: residual = sqrt(n) * sigma
alphas = np.logspace(-10, 0, 201)
res = np.array([np.linalg.norm(K @ tikhonov(a) - y) for a in alphas])
a_mor = alphas[np.argmin(np.abs(res - np.sqrt(n) * sigma))]
errs = np.array([err(tikhonov(a)) for a in alphas])
print(f"Tikhonov, Morozov alpha={a_mor:.1e}: rel. error {err(tikhonov(a_mor)):.3f}")
print(f"Tikhonov, oracle  alpha={alphas[errs.argmin()]:.1e}: rel. error {errs.min():.3f}")
for r in [10, 20, 40]:
    print(f"TSVD r={r:2d}: rel. error {err(Vt[:r].T @ (c[:r] / s[:r])):.3f}")

# Bayesian reading: prior N(0, gamma^2 I) with gamma^2 = sigma^2 / alpha
gamma = sigma / np.sqrt(a_mor)
post_sd = 1 / np.sqrt(s**2 / sigma**2 + 1 / gamma**2)   # posterior s.d. along each v_j
print(f"directions informed by the data (s_j^2 > alpha): {np.sum(s**2 > a_mor)} of {n}")
print(f"posterior s.d. along v_1, v_20, v_100: {post_sd[0]:.1e}, {post_sd[19]:.1e}, {post_sd[99]:.1e}")

3. Time-reversal imaging

The first application of ch. 2 is imaging itself. Least squares on a fine grid of pixels is hopelessly ill posed; its strongly regularised limit, back-propagation of the time-reversed data, is stable, and its resolution can be computed exactly. This chapter derives migration from ch. 2, explains physical time reversal through the Helmholtz–Kirchhoff identity, and ends with the subspace methods DORT and MUSIC.

Migration as strongly regularised least squares

Fix a search grid z_1, \dots, z_P of candidate positions with unknown reflectivities \rho \in \mathbb{C}^P. By Proposition 1.7 the data depend linearly on \rho: \mathcal{K} \rho = k^2 \sum_{p=1}^P \rho_p\, g(z_p)\, g(z_p)^\top \in \mathbb{C}^{N \times N} .

With the Frobenius inner product on matrices, the adjoint is (\mathcal{K}^* A)_p = k^2\, g(z_p)^*\, A\, \overline{g(z_p)} .

Migration is therefore least squares with so much regularisation that only the adjoint survives. By Theorem 2.8 its variance is under control; its bias is a blur, since the image of the exact data is \mathcal{K}^* \mathcal{K} \rho. The blur kernel is the point spread function computed below.

This delay-and-sum formula is the beamforming used in ultrasound scanners and the Kirchhoff migration of seismology.

The time-reversal experiment

Experiments by Fink and coworkers show that the re-emitted wave refocuses on the source. The next theorem explains why, and how sharp the focus is.

The focal spot cannot be smaller than half a wavelength, whatever the mirror: this is the diffraction limit. In the time domain, \mathrm{Im}\, \hat G corresponds to G(t) - G(-t), a wave converging on the source and then diverging from it.

Stability: the signal-to-noise ratio of migration

Iterative time reversal, DORT and MUSIC

Re-emitting the conjugate of the recorded field twice maps an input v to \bar A A v = A^* A v, since A is symmetric. Iterating time reversal is thus the power method on the time-reversal operator A^* A: it converges to the top singular vector, and the array focuses on the strongest reflector. The DORT method of Prada and Fink uses all the singular vectors at once.

Orthogonality holds approximately when the reflectors are farther apart than the Rayleigh distance of Proposition 3.7. Closer reflectors mix their singular vectors, but not their signal subspace, which MUSIC exploits.

With noisy data, the range is replaced by the span of the R leading left singular vectors. This needs an estimate of R, which ch. 5 obtains from random matrix theory.

MUSIC comes from array processing (Schmidt, 1986); Cheney (2001) relates it to the linear sampling method of inverse scattering.

import numpy as np
from scipy.special import hankel1

rng = np.random.default_rng(0)
N = 41
xs = np.stack([(np.arange(N) - 20) * 0.5, np.zeros(N)], axis=1)   # aperture a = 20 (wavelengths)
refl = [(np.array([-2.0, 30.0]), 1.0), (np.array([3.0, 32.0]), 0.6)]
ks = 2 * np.pi * np.linspace(0.8, 1.2, 21)                         # band B = 0.4, centre lambda = 1
sigma = 0.01

def g(k, z):                     # illumination vectors, z of shape (..., 2) -> (..., N)
    r = np.linalg.norm(z[..., None, :] - xs, axis=-1)
    return 0.25j * hankel1(0, k * r)

def msr(k):                      # Born response matrix + measurement noise
    A = sum(rho * k**2 * np.outer(g(k, z), g(k, z)) for z, rho in refl)
    return A + sigma * (rng.standard_normal((N, N)) + 1j * rng.standard_normal((N, N))) / np.sqrt(2)

data = [msr(k) for k in ks]
X, Y = np.meshgrid(np.arange(-6, 7.01, 0.1), np.arange(26, 36.01, 0.1))
Z = np.stack([X, Y], axis=-1)

def kirchhoff(freqs):            # I(z) = sum_k g(z)^* A_k conj(g(z))
    I = 0
    for i in freqs:
        G = g(ks[i], Z)
        I = I + np.einsum("...n,nm,...m->...", G.conj(), data[i], G.conj())
    return np.abs(I)

I_multi = kirchhoff(range(len(ks)))
iy, ix = np.unravel_index(I_multi.argmax(), I_multi.shape)
print(f"migration peak at ({X[iy, ix]:.1f}, {Y[iy, ix]:.1f})")

# cross-range and range profiles through the first reflector
I_one = kirchhoff([10])                                   # centre frequency only
row = I_one[np.argmin(np.abs(Y[:, 0] - 30.0))]
xr = X[0]; j0 = np.argmin(np.abs(xr + 2.0))
first_null = xr[j0 + np.argmax(np.diff(row[j0:]) > 0)] + 2.0
print(f"cross-range first null at {first_null:.1f} (Rayleigh lambda L / a = {30 / 20:.1f})")
col = I_multi[:, np.argmin(np.abs(xr + 2.0))]
y0 = np.argmin(np.abs(Y[:, 0] - 30.0))
print(f"range first null at {Y[y0 + np.argmax(np.diff(col[y0:]) > 0), 0] - 30:.1f} (c / 2B = {1 / (2 * 0.4):.2f})")

# MUSIC at the centre frequency, signal space of dimension 2
U = np.linalg.svd(data[10])[0][:, :2]
G = g(ks[10], Z); G /= np.linalg.norm(G, axis=-1, keepdims=True)
proj = np.linalg.norm(G - (G @ U.conj()) @ U.T, axis=-1)
music = 1 / proj**2
for half in [X < 0.5, X >= 0.5]:                    # one search window per half-plane
    m = np.where(half, music, 0)
    iy, ix = np.unravel_index(m.argmax(), m.shape)
    print(f"MUSIC peak at ({X[iy, ix]:.1f}, {Y[iy, ix]:.1f})")

In a randomly scattering medium, physical time reversal refocuses even better than in a homogeneous one: multiple scattering widens the effective aperture, and the refocused field is statistically stable (Fink, 1997; Blomgren, Papanicolaou and Zhao, 2002). Migration of the recorded data, by contrast, is blurred by clutter; ch. 8 returns to this contrast.

4. Random matrix theory

Ch. 1 and ch. 3 showed that a few singular values of a noisy response matrix stand out of a noise bulk. To turn this observation into a detection test, one needs the law of the singular values of a pure-noise matrix, and the exact way a weak signal emerges from it. Random matrix theory provides both, in the limit of many sensors.

Throughout, W is an N \times N matrix with independent and identically distributed entries of mean zero and variance \mathbb{E}|W_{ij}|^2 = 1, and X = W / \sqrt N.

The Marchenko–Pastur law

The bulk of noise singular values thus fills [0, 2\sqrt N] for W itself, with the quarter-circle profile. No noise direction is privileged, by the unitary invariance of Proposition 0.7.

The convergence of the largest singular value is due to Geman (1980) and Yin, Bai and Krishnaiah (1988). The law \mathrm{TW}_2 was found by Tracy and Widom (1994) for Gaussian Hermitian matrices. Its appearance for the largest eigenvalue of W W^* is due to Johansson (2000) in the complex case and Johnstone (2001) in the real case; the formula above follows by taking a square root.

The spiked model and the BBP transition

A reflector adds a rank-one matrix to the noise. The question is whether its singular value separates from the bulk.

Below the threshold \theta = 1, the signal is invisible in the spectrum: the top singular value sticks to the edge and its vector is asymptotically orthogonal to u. Above it, the singular value separates and is biased upwards by 1/\theta.

For a reflector, s_A = |\rho| k^2 \|g\|^2 grows like N while the threshold \sigma\sqrt N grows like \sqrt N: adding sensors always ends up revealing a fixed reflector.

import numpy as np

rng = np.random.default_rng(2)
def noise(N):                    # i.i.d. CN(0, 1) entries
    return (rng.standard_normal((N, N)) + 1j * rng.standard_normal((N, N))) / np.sqrt(2)

# 1. Quarter-circle law for the singular values of W / sqrt(N)
N = 400
s = np.linalg.svd(noise(N) / np.sqrt(N), compute_uv=False)
grid = np.linspace(0, 2, 2001)
cdf = np.cumsum(np.sqrt(4 - grid**2) / np.pi) * (grid[1] - grid[0])
ks = np.max(np.abs(np.searchsorted(np.sort(s), grid) / N - cdf))
print(f"N={N}: largest singular value {s[0]:.3f} (limit 2), Kolmogorov distance {ks:.3f}")

# 2. Edge fluctuations: sigma_1 ~ 2 sqrt(N) + 2^(-2/3) N^(-1/6) TW2, E[TW2] = -1.771
for N in [50, 200]:
    top = np.array([np.linalg.norm(noise(N), 2) for _ in range(300)])
    pred = 2 * np.sqrt(N) + 2**(-2/3) * N**(-1/6) * (-1.771)
    print(f"N={N}: mean sigma_1 {top.mean():.3f}, Tracy-Widom prediction {pred:.3f}, sd {top.std():.3f}")

# 3. Spiked model X = theta u v^* + W / sqrt(N): BBP transition
N = 400
u = np.ones(N) / np.sqrt(N); v = np.exp(2j * np.pi * rng.random(N)) / np.sqrt(N)
for theta in [0.5, 1.0, 1.5, 2.0, 3.0]:
    Uh, sh, Vh = np.linalg.svd(theta * np.outer(u, v.conj()) + noise(N) / np.sqrt(N))
    s_lim = theta + 1 / theta if theta > 1 else 2.0
    ov_lim = max(1 - 1 / theta**2, 0.0)
    print(f"theta={theta:.1f}: sigma_1={sh[0]:.3f} (limit {s_lim:.3f}), "
          f"|<u,u_hat>|^2={abs(u @ Uh[:, 0].conj())**2:.2f} (limit {ov_lim:.2f})")

5. Detection tests and what follows them

The first imaging question is whether there is anything to image. This chapter builds tests with a prescribed false-alarm rate from the response matrix, computes their power with ch. 4, and then localises, characterises and identifies what was detected.

Throughout, H_0 says A^{\mathrm{meas}} = W and H_1 says A^{\mathrm{meas}} = A + W, with W as in Definition 1.8 and A the MSR matrix of one reflector.

Tests and the Neyman–Pearson lemma

In imaging, H_1 contains unknown parameters: the position and the reflectivity. The generalised likelihood-ratio test (GLRT) replaces them by their maximum-likelihood estimates. Under the Gaussian model, minus the log-likelihood is \|A^{\mathrm{meas}} - A\|_F^2 / \sigma^2 plus a constant.

Known position: migration is the matched filter

Unknown position: the top singular value

The test ignores the physics of the illumination vectors, which makes it robust: it also detects targets that are not points, or a model with an inaccurate speed of sound.

When \sigma is unknown, it is estimated from the bulk, for instance by \hat\sigma^2 = \frac{1}{N^2} \sum_{j \geq 2} \sigma_j(A^{\mathrm{meas}})^2 , or by fitting the quarter-circle law of Theorem 4.2 to the empirical singular values. Removing the top value makes the estimate consistent under both hypotheses.

Below the threshold, the top singular value is blind, but the data are not. The energy detector, which uses the whole spectrum, keeps a nontrivial power there.

The two tests are complementary: the energy detector never reaches power one at fixed \theta, while the top singular value is asymptotically perfect above the threshold. Onatski, Moreira and Hallin (2013) study the best achievable power below it.

After detection: localisation and characterisation

Once a target is declared, the top singular triplet carries its position and strength.

• Localisation. Image with the top left singular vector, \mathcal{I}_{\mathrm{SV}}(z) = |\hat g(z)^* \hat u_1|^2, or with MUSIC (Theorem 3.11). By Theorem 4.4 the peak value is about 1 - 1/\theta^2 instead of 1, so localisation degrades smoothly as \theta decreases to 1.
• Number of targets. The number of singular values above \tau_\alpha estimates the rank R needed by MUSIC.
• Identification. Comparing estimated features, such as reflectivity or, for extended targets, polarisation tensors, with a dictionary of known targets turns identification into classification; Proposition 5.11 below computes its error rate.

import numpy as np

rng = np.random.default_rng(3)
N, sigma, n_mc = 100, 1.0, 400
def noise():
    return sigma * (rng.standard_normal((N, N)) + 1j * rng.standard_normal((N, N))) / np.sqrt(2)

u = np.exp(2j * np.pi * rng.random(N)) / np.sqrt(N)      # unit "illumination" vector
v = u.conj()                                              # MSR matrices are symmetric: v = conj(u)
stats = {
    "top singular value": lambda M: np.linalg.norm(M, 2),
    "Frobenius energy":   lambda M: np.linalg.norm(M),
    "matched filter":     lambda M: abs(u.conj() @ M @ v),   # known position
}
h0 = {name: np.array([f(noise()) for _ in range(n_mc)]) for name, f in stats.items()}
thr = {name: np.quantile(x, 0.95) for name, x in h0.items()}          # false-alarm rate 5 %
edge = sigma * (2 * np.sqrt(N) + 2**(-2/3) * N**(-1/6) * (-0.232))    # TW2 95 % quantile ~ -0.232
print(f"threshold on sigma_1: simulated {thr['top singular value']:.2f}, Tracy-Widom {edge:.2f}")
for theta in [0.2, 0.8, 1.2, 1.6]:
    s_ref = theta * sigma * np.sqrt(N)                    # singular value of the clean matrix
    A = s_ref * np.outer(u, v.conj())
    pd = {name: np.mean([f(A + noise()) > thr[name] for _ in range(n_mc)]) for name, f in stats.items()}
    print(f"theta={theta}: " + ", ".join(f"{k} {p:.2f}" for k, p in pd.items()))

# characterisation: invert sigma_1 = s + sigma^2 N / s
s_ref = 2.0 * sigma * np.sqrt(N)
s1 = np.linalg.norm(s_ref * np.outer(u, v.conj()) + noise(), 2)
s_hat = (s1 + np.sqrt(s1**2 - 4 * sigma**2 * N)) / 2
print(f"true {s_ref:.1f}, raw sigma_1 {s1:.2f}, debiased {s_hat:.2f}")

Identification within a dictionary

Identification asks which of several known targets is present, for instance which of two materials fills a small inclusion at a known position. Under the Gaussian noise model it is a classification problem with an exact error rate.

6. Spectral properties of stationary random processes

Sensor noise and ambient noise are random signals whose statistics do not change in time. Their spectral description answers the two questions of ch. 7: what does the correlation of two recordings estimate, and how long must one record for the estimate to be reliable?

Stationary processes

The power spectral density

The spectral density says how the variance C(0) = \frac{1}{2\pi}\int S is distributed over frequencies. Conversely, Bochner’s theorem states that every continuous positive-definite function is the Fourier transform of a non-negative measure.

Filtering and cross-spectra

In ch. 7 the medium is the filter: the signal recorded at a sensor is the noise of each source filtered by the Green’s function from the source to the sensor.

Fluctuations of empirical correlations

In practice the covariance is unknown and is estimated from a single record of duration T. For Gaussian processes the error of this estimate can be computed exactly, thanks to the following moment formula.

A record of duration T and bandwidth B contains about BT independent pieces of information. One day of seismic noise in a band of 1 Hz gives a relative error of about 1/\sqrt{86\,400} \approx 0.3\,\% on a correlation of order C(0).

The periodogram

Averaging the periodograms of K independent segments, as in Bartlett’s and Welch’s methods, divides the variance by K at the price of a frequency resolution of order 2\pi K / T. The same computation shows that Fourier coefficients at frequencies more than 2\pi/T apart are asymptotically independent, which justifies the noise model of Definition 1.8.

import numpy as np

rng = np.random.default_rng(4)
dt = 1e-3                                         # sampling step (s)
def bandlimited_noise(n_samples, f_lo=20.0, f_hi=40.0):
    """Gaussian noise with a flat spectrum on f_lo < |f| < f_hi (Hz), unit variance."""
    w = rng.standard_normal(n_samples)
    f = np.fft.rfftfreq(n_samples, dt)
    x = np.fft.irfft(np.fft.rfft(w) * ((f > f_lo) & (f < f_hi)), n_samples)
    return x / np.sqrt(2 * (f_hi - f_lo) * dt)  # variance = 2 B dt * (1 / (2 B dt)) = 1

B = 20.0
for T in [1.0, 4.0, 16.0]:
    n = int(T / dt)
    c0 = np.array([np.mean(bandlimited_noise(n)**2) for _ in range(200)])   # C_T(0)
    print(f"T={T:4.0f} s: sd of C_T(0) = {c0.std():.3f}, theory 1/sqrt(BT) = {1 / np.sqrt(B * T):.3f}")

# the periodogram does not converge; averaging over segments (Welch) does
x = bandlimited_noise(2**16)
for n_seg in [1, 64]:
    segs = x.reshape(n_seg, -1)
    P = np.mean(np.abs(np.fft.rfft(segs, axis=1))**2, axis=0) * dt / segs.shape[1]
    f = np.fft.rfftfreq(segs.shape[1], dt)
    band = (f > 22) & (f < 38)
    print(f"{n_seg:2d} segment(s): mean of P in band {P[band].mean():.4f}, "
          f"relative sd {P[band].std() / P[band].mean():.2f}")
print("true S(omega) in band:", 1 / (2 * B))

7. Imaging with ambient noise cross-correlations

Ambient noise, generated by ocean waves, wind or traffic, is recorded by every seismometer, all the time. This chapter shows that correlating the noise recorded at two sensors produces, on average, the signal that one sensor would record if the other were a source. Passive sensors can then measure travel times and image reflectors without any controlled source.

The noise-source model

Green’s function retrieval

The same proof works in a heterogeneous medium, as long as the speed is constant outside a ball enclosed by the source sphere: Green’s identity only uses that the Helmholtz operator is real. The cross-correlation then contains the full Green’s function, with all its scattered arrivals.

Directivity of the sources

Real noise sources are not uniformly distributed: ocean microseisms come mainly from the coasts. The far-field computation below shows which sources matter.

For a band-pass noise, F oscillates and the last integral is small when \bar K is smooth. The pulse at \tau = +T is weighted by the sources in the direction e, beyond x_1 as seen from x_2.

The pulse at \tau = -T is weighted by the sources beyond x_2.

With \bar K \equiv 1 one recovers Theorem 7.3. This is the stationary-phase argument of Snieder (2004), in exact form for a homogeneous medium.

Travel times and tomography

The peak of the cross-correlation envelope gives the travel time between two sensors. Its reliability follows from ch. 6: by Theorem 6.8, the empirical cross-correlation fluctuates around its mean with a standard deviation of order C(0)/\sqrt{BT}, where C(0) is the power of the recordings. A coherent peak of relative amplitude \eta therefore emerges after a recording time of order 1/(\eta^2 B).

With travel times T_{ij} between many pairs of stations, the slowness 1/c(x) can be reconstructed. Linearised around a reference model, each travel time is a line integral of the slowness perturbation along a ray, so the problem is a linear inverse problem of ch. 2, regularised by Tikhonov’s method. Shapiro, Campillo, Stehly and Ritzwoller (2005) applied this to surface waves across California.

Passive imaging

Theorem 7.3 says that a pair of passive sensors behaves, on average, like an active source–receiver pair. An array of passive sensors therefore yields a virtual response matrix, and the imaging methods of ch. 3 apply.

Migrating the positive-lag cross-correlations with the delay-and-sum formula of Proposition 3.3 images the reflector, with the resolution of Proposition 3.7. Garnier and Papanicolaou (2009) analyse this passive imaging, including its statistical stability and the role of the source distribution. When the illumination is partial, the missing directions produce artefacts, a topic taken up in ch. 8.

import numpy as np

c0, R, S = 1.0, 50.0, 6000                       # speed, radius of the source sphere, sources
x1, x2 = np.array([-2.6, 0, 0]), np.array([2.6, 0, 0])
i = np.arange(S) + 0.5                            # Fibonacci points on the sphere
phi, th = np.pi * (1 + 5**0.5) * i, np.arccos(1 - 2 * i / S)
Y = R * np.stack([np.cos(phi) * np.sin(th), np.sin(phi) * np.sin(th), np.cos(th)], axis=1)
dA = 4 * np.pi * R**2 / S                         # area element of each source
r1, r2 = np.linalg.norm(Y - x1, axis=1), np.linalg.norm(Y - x2, axis=1)
G = lambda k, r: np.exp(1j * k * r) / (4 * np.pi * r)
r12 = np.linalg.norm(x1 - x2)

# Helmholtz-Kirchhoff: sum_y conj G(x1,y) G(x2,y) dA  ~  Im G(x1,x2) / k
for k in [2 * np.pi, 3 * np.pi]:
    lhs = np.sum(np.conj(G(k, r1)) * G(k, r2)) * dA
    print(f"k={k:.2f}: source sum {lhs.real:+.6f}{lhs.imag:+.6f}i, "
          f"Im G / k = {np.sin(k * r12) / (4 * np.pi * r12) / k:+.6f}")

# envelope of the cross-correlation for a noise spectrum centred on 1 Hz
w = np.linspace(0.05, 4 * np.pi, 800); dw = w[1] - w[0]
F = np.exp(-((w - 2 * np.pi) / 2.0)**2)
tau = np.linspace(-10, 10, 2001)
for name, K in [("uniform", np.ones(S)), ("mostly x < 0", (1 - Y[:, 0] / R)**2)]:
    Chat = F * ((K * np.conj(G(w[:, None] / c0, r1)) * G(w[:, None] / c0, r2)).sum(1) * dA)
    env = np.abs(np.exp(-1j * np.outer(tau, w)) @ Chat) * dw / np.pi   # analytic-signal envelope
    pos, neg = tau > 0, tau < 0
    print(f"sources {name:12s}: envelope max at tau={tau[pos][env[pos].argmax()]:+.2f} "
          f"({env[pos].max():.1e}) and tau={tau[neg][env[neg].argmax()]:+.2f} ({env[neg].max():.1e})")
print(f"travel time |x1 - x2| / c0 = {r12 / c0:.2f}")

8. Research directions and open problems

This chapter is for self-study. It sketches where the four themes of the module lead in current research, with entry points into the literature; it is a map, not a survey.

Imaging in random and cluttered media

The additive noise model of Definition 1.8 allows exact computations, but in echography and seismology the dominant randomness is clutter: echoes from many small heterogeneities of the medium. Migration of cluttered data produces speckle, images whose fine details depend on the unknown realisation of the medium.

Coherent interferometry, developed by Borcea, Papanicolaou and Tsogka (2005), migrates local cross-correlations of the data instead of the data themselves. The resulting images are statistically stable, at the price of some resolution, which is the trade-off of Theorem 6.8 transposed to imaging.

In ultrasound, Aubry and Derode (2009) used random matrix theory to separate single scattering, which carries the image, from multiple scattering, which behaves like a random matrix. Lambert, Cobus, Frappart, Fink and Aubry (2020) introduced a distortion matrix that corrects aberrations caused by an unknown speed of sound.

Time reversal behaves differently from migration in random media: multiple scattering enlarges the effective aperture and the refocusing is statistically stable (Blomgren, Papanicolaou and Zhao, 2002). The asymptotic theory of waves in randomly layered media is developed in the book by Fouque, Garnier, Papanicolaou and Sølna (2007).

High-dimensional detection and estimation

Ch. 4 and ch. 5 used the simplest spiked model. Several directions refine it.

• Denoising before imaging. Gavish and Donoho (2014) showed that, for a square matrix with known noise level, the optimal hard threshold for singular values is 4/\sqrt 3 times \sqrt N \sigma. Optimal shrinkage of the singular values themselves (Gavish and Donoho, 2017) improves on hard thresholding, and Donoho, Gavish and Johnstone (2018) treat the analogous shrinkage of eigenvalues in spiked covariance models.
• Estimating the number of targets. Information criteria, introduced for array processing by Wax and Kailath (1985), compete with thresholds based on the BBP transition.
• Below the threshold. Onatski, Moreira and Hallin (2013) compute the asymptotic power of tests that use the whole spectrum, the regime where the top singular value is blind.
• Realistic noise. MSR noise is often symmetric by reciprocity, correlated between neighbouring sensors, or non-Gaussian. Each changes the bulk, the edge or the transition, and finite-N corrections matter for arrays of a few dozen sensors.

The official reference by Ammari and coauthors (2013) develops multistatic imaging of small and extended targets along these lines.

Noise-based imaging and monitoring

The founding experiments of correlation-based imaging include ultrasonic thermal noise in an aluminium block (Weaver and Lobkis, 2001) and the seismic coda (Campillo and Paul, 2003). Wapenaar (2004) proved Green’s function retrieval for arbitrary inhomogeneous elastic media, and Snieder (2004) gave the stationary-phase explanation of Proposition 7.4. Bakulin and Calvert (2006) turned the idea into the virtual source method of exploration geophysics.

Repeated cross-correlations also monitor slow changes of a medium. Brenguier and coauthors (2008) measured small relative changes of seismic velocity along the San Andreas fault after earthquakes, using years of continuous noise records. The theory of these measurements combines ch. 6 and ch. 7 with the analysis of small perturbations of travel times.

The book by Garnier and Papanicolaou (2016), an official reference of the module, gives the mathematical theory of passive imaging, including scattering media, where multiple scattering can help by diversifying the illumination.

Learned reconstruction and generative priors

The Bayesian framework of ch. 2 extends to priors learned from data. Arridge, Maass, Öktem and Schönlieb (2019) survey data-driven models for inverse problems. Jin, McCann, Froustey and Unser (2017) post-process a simple reconstruction with a convolutional network, and Adler and Öktem (2018) unroll an iterative scheme into a trained network.

Diffusion models provide generative priors, and Chung, Kim, McCann, Klasky and Ye (2023) sample approximate posteriors with them. Whether such posteriors quantify uncertainty reliably, in the sense of Proposition 2.12, is largely open. The MVA courses on image restoration and synthesis and on diffusion and multiscale models develop these models.

Full-waveform inversion and super-resolution

Migration linearises the problem; full-waveform inversion fits the complete wave equation to the data by nonlinear least squares (Tarantola, 1984; Virieux and Operto, 2009). The misfit is not convex, and a poor initial model converges to a wrong one whose arrivals are shifted by a period, the cycle-skipping problem. Engquist and Froese (2014) proposed optimal-transport distances between seismic signals to convexify the misfit.

The diffraction limit of Corollary 3.6 holds for linear imaging. With a sparsity prior, Candès and Fernandez-Granda (2014) showed that point sources separated by a few times the Rayleigh distance can be recovered exactly by convex optimisation from low-frequency data.

A list of open problems

1. Detection thresholds and power for structured noise (symmetric, correlated, heavy-tailed) at moderate N.
2. Statistical stability of passive imaging with non-uniform, non-stationary noise sources.
3. Calibrated uncertainty for learned reconstructions and generative priors.
4. What information about the medium survives strong multiple scattering, and how to extract it.
5. Joint estimation of the noise level, the number of targets and their positions, with finite-sample guarantees.

Exercises

Levels: ★ apply the course, ★★ combine or prove a short result, ★★★ open-ended or research-flavoured. Exercises 12 and 14 are computational.

1. ★ Resolution of a seismic survey. Waves travel at c_0 = 3000 m/s with central frequency 30 Hz and bandwidth 20 Hz. A linear array of aperture a = 1 km images a reflector at depth L = 3 km. Compute the wavelength, the cross-range resolution and the range resolution of echoes. (ch. 1, 3)

   Solution: \lambda = 3000/30 = 100 m. By Proposition 3.7, the cross-range resolution is \lambda L / a = 100 \times 3 = 300 \text{ m}, and the range resolution of echoes is c_0/(2B) = 3000/40 = 75 m. Range is resolved four times better than cross-range: with a narrow aperture, bandwidth is the main source of resolution.

2. ★ Eigenvalues are not singular values. Three sensors see one reflector with illumination vector g = (1, i, 1)^\top and \rho k^2 = 2. Compute the MSR matrix, its singular values and vectors, the migration image at the reflector, and the nonzero eigenvalue of A. (ch. 1, 3)

   Solution: A = 2 g g^\top has rank one, and \|g\|^2 = 3, so by Proposition 1.7 \sigma_1 = 2 \times 3 = 6, \qquad u_1 = g/\sqrt 3, \qquad v_1 = \bar g/\sqrt 3 = (1, -i, 1)^\top/\sqrt 3 .

   The migration image is g^* A \bar g = 2 \|g\|^4 = 18. However A g = 2 g\, (g^\top g) and g^\top g = 1 - 1 + 1 = 1, so the nonzero eigenvalue is 2, not 6. A complex symmetric matrix is not Hermitian: imaging uses its singular values, or the eigenvalues of the Hermitian time-reversal operator A^* A.

3. ★ Filter factors. Let s_j = 10^{-j} for j = 0, \dots, 4 and \alpha = 10^{-4}. Compute the Tikhonov factors \varphi_j and the noise amplification \varphi_j / s_j, and compare with TSVD at the same \alpha. (ch. 2)

   Solution: \varphi_j = s_j^2/(s_j^2 + \alpha) gives 0.9999, 0.990, 0.5, 0.0099 and 0.0001. The amplifications \varphi_j / s_j are 1.0, 9.9, 50, 9.9 and 1.0, never above \frac{1}{2\sqrt\alpha} = 50 , as in Theorem 2.8. TSVD keeps j = 0, 1, 2 with amplifications 1, 10 and 100, and discards the rest: it has twice the worst-case amplification and a sharper cut.

4. ★ Exponential covariance. Let C(\tau) = s^2 e^{-|\tau|/\tau_c}. Compute the spectral density, and the relative standard deviation of C_T(0) for a Gaussian process with this covariance. (ch. 6)

   Solution: by evenness of C, S(\omega) = 2\, \mathrm{Re} \int_0^\infty s^2 e^{-\tau(1/\tau_c - i\omega)}\, d\tau , that is S(\omega) = \frac{2 s^2 \tau_c}{1 + \omega^2 \tau_c^2} , a Lorentzian of width 1/\tau_c. By Theorem 6.8, T\, \mathrm{Var}\, C_T(0) \to 2\int C^2 = 2 s^4 \tau_c, so the relative standard deviation is \sqrt{2\tau_c / T}: the record must be much longer than the correlation time.

5. ★ Detection by numbers. An array of N = 64 sensors has noise level \sigma = 0.02 per entry. Give the edge of the noise bulk and the detectability threshold. For a reflector with s_A = 0.4, predict the top singular value and the overlap of its vector. If the observed top singular value is 0.47, estimate s_A. (ch. 4, 5)

   Solution: the edge is 2\sigma\sqrt N = 0.32 and the threshold \sigma\sqrt N = 0.16. With \theta = 0.4/0.16 = 2.5, Corollary 4.6 and Theorem 4.4 give \sigma_1 \approx 0.4 + \frac{0.0004 \times 64}{0.4} = 0.464, \qquad |u^* \hat u_1|^2 \approx 1 - \frac{1}{6.25} = 0.84 .

   Proposition 5.8 gives \hat s_A = (0.47 + \sqrt{0.2209 - 0.1024})/2 = 0.407.

6. ★ Correlations of white noise. Let n_t be independent \mathcal{N}(0,1) and C_T(\tau) = T^{-1}\sum_{t=1}^T n_t\, n_{t+\tau}. Compute the mean and variance of C_T(\tau). (ch. 6)

   Solution: the mean is \mathbb{1}[\tau = 0]. For \tau \neq 0 the products n_t n_{t+\tau} are uncorrelated with variance 1, so the variance is 1/T.

   For \tau = 0, \mathrm{Var}(n_t^2) = 2 gives 2/T.

   This is the discrete form of Theorem 6.8, where the integral becomes \sum_s C(s)^2 + C(s+\tau) C(s-\tau).

7. ★★ Uniqueness in the discrepancy principle. Show that \alpha \mapsto \|x_\alpha\| is non-increasing for Tikhonov regularisation, and that the discrepancy equation of Proposition 2.9 has a unique solution when some u_j^* y with s_j > 0 is nonzero. (ch. 2)

   Solution: in the singular bases, \|x_\alpha\|^2 = \sum_j \frac{s_j^2\, |u_j^* y|^2}{(s_j^2 + \alpha)^2} , a sum of non-increasing functions of \alpha. In the residual of Proposition 2.9, each term with s_j > 0 and u_j^* y \neq 0 is strictly increasing in \alpha, so r(\alpha) is strictly increasing and takes the target value once.

8. ★★ Bayes risk. With the prior x \sim \mathcal{N}(0, \gamma^2 I), compute the Bayes risk of the posterior mean and of least squares in terms of the s_j. Evaluate them for s_j = 10^{-j}, j = 0, \dots, 4, \sigma = 0.01 and \gamma = 1. (ch. 2)

   Solution: by Proposition 2.12 the posterior mean has risk \mathrm{tr}\, C, and least squares has risk \sigma^2 \sum_j s_j^{-2} for every x: \mathrm{tr}\, C = \sum_j \frac{\sigma^2 \gamma^2}{\sigma^2 + \gamma^2 s_j^2}, \qquad R_{\mathrm{LS}} = \sigma^2 \sum_j s_j^{-2} .

   Numerically, the terms of \mathrm{tr}\, C are 0.0001, 0.0099, 0.5, 0.990 and 0.9999, so \mathrm{tr}\, C \approx 2.50. The prior mean alone has risk 5\gamma^2 = 5, and least squares about 1.01 \cdot 10^4.

9. ★★ Power of the known-position test. In Proposition 5.3, let s = |\rho| k^2 \|g\|^2 and t = \sqrt{\ln(1/\alpha)}. Show that, if s/\sigma > t, P_{\mathrm{D}} \geq 1 - \exp\big(-(s/\sigma - t)^2\big) , and evaluate the bound for \alpha = 5\,\% and s/\sigma = 4. Compare with the threshold of blind detection. (ch. 3, 5)

   Solution: under H_1 the normalised statistic is |s e^{i\phi}/\sigma + Z'| with Z' \sim \mathcal{CN}(0,1). It falls below t only if |Z'| \geq s/\sigma - t.

   Since |Z'|^2 is exponential, this event has probability \exp(-(s/\sigma - t)^2).

   With t = 1.731 and s/\sigma = 4, the bound gives P_{\mathrm{D}} \geq 1 - e^{-5.15} = 0.994.

   Blind detection needs s > \sigma\sqrt N, that is s/\sigma > 10 for N = 100: knowing the position saves a factor of order \sqrt N.

10. ★★ Focal spot in two dimensions. Show that a circular time-reversal mirror surrounding a source in two dimensions refocuses a field proportional to J_0(k|x - y|), and compute the radius of the first zero. Compare with three dimensions. (ch. 3)

    Solution: the proof of Theorem 3.5 holds in two dimensions, since \hat G = O(r^{-1/2}) and \partial_r \hat G - ik\hat G = O(r^{-3/2}) make the boundary error O(R^{-1}). The field is proportional to \mathrm{Im}\, \hat G = \mathrm{Im}\, \frac{i}{4} H_0^{(1)}(kr) = \frac{1}{4} J_0(kr) .

    The first zero of J_0 is at 2.405, hence r = 2.405\, \lambda / 2\pi \approx 0.38\, \lambda, about three quarters of the \lambda/2 of three dimensions.

11. ★★ The second moment of the Marchenko–Pastur law. Show directly that N^{-1}\, \mathbb{E}\, \mathrm{tr}\, (X X^*)^2 \to 2 for X = W/\sqrt N, and check that C_2 = 2. (ch. 4)

    Solution: expand \mathbb{E}\, \mathrm{tr}\, (W W^*)^2 = \sum_{i, j, k, l} \mathbb{E}\big[W_{ij}\, \overline{W_{kj}}\, W_{kl}\, \overline{W_{il}}\big] .

    The expectation vanishes unless the entries pair up. Either k = i, which gives \sum_{i,j,l} \mathbb{E}\, |W_{ij}|^2 |W_{il}|^2 = N^3 + O(N^2) , or l = j with k \neq i, which gives another N^3 + O(N^2).

    Dividing by N \cdot N^2 gives the limit 2, the Catalan number C_2 = \binom42/3 = 2.

12. ★★ Localisation by the top singular vector (computational). For the array of Code 1.10 and one reflector, choose the noise level so that \theta = 0.7, 1.5, 3, 6. Compare the top singular value with Corollary 4.6, the overlap |\hat g(z_1)^* \hat u_1|^2 with 1 - 1/\theta^2, and the cross-range position of the maximum of |\hat g(z)^* \hat u_1|. (ch. 4, 5)

    Hints and expected result: the code below gives ratios of \sigma_1 to the prediction between 0.97 and 1.04. The overlaps are 0.02, 0.71, 0.87 and 0.96, against 0, 0.56, 0.89 and 0.97: at N = 41, finite-size effects are visible near the threshold. Localisation fails at \theta = 0.7, is off by 0.15\lambda at \theta = 1.5, and is exact on the grid for \theta \geq 3.

import numpy as np
from scipy.special import hankel1

rng = np.random.default_rng(5)
N, k = 41, 2 * np.pi
xs = np.stack([(np.arange(N) - 20) * 0.5, np.zeros(N)], axis=1)
g = lambda z: 0.25j * hankel1(0, k * np.linalg.norm(z[..., None, :] - xs, axis=-1))
z1 = np.array([-2.0, 30.0])
A = k**2 * np.outer(g(z1), g(z1))                       # one reflector, rho = 1
s_A = np.linalg.norm(A, 2)
line = np.stack([np.linspace(-8, 4, 241), np.full(241, 30.0)], axis=1)   # cross-range search line
G_line = g(line); G_line /= np.linalg.norm(G_line, axis=1, keepdims=True)
gh = g(z1) / np.linalg.norm(g(z1))
for theta in [0.7, 1.5, 3.0, 6.0]:
    sigma = s_A / (theta * np.sqrt(N))                  # noise level giving this theta
    W = sigma * (rng.standard_normal((N, N)) + 1j * rng.standard_normal((N, N))) / np.sqrt(2)
    U, s, Vh = np.linalg.svd(A + W)
    pred = s_A + sigma**2 * N / s_A if theta > 1 else 2 * sigma * np.sqrt(N)
    overlap = abs(gh.conj() @ U[:, 0])**2
    x_hat = line[np.argmax(np.abs(G_line.conj() @ U[:, 0])), 0]
    print(f"theta={theta}: sigma_1/pred={s[0] / pred:.3f}, overlap={overlap:.2f} "
          f"(limit {max(1 - 1 / theta**2, 0):.2f}), estimated cross-range {x_hat:+.2f} (true -2.00)")

13. ★★ One-sided noise. In Proposition 7.4, let \bar K(\mu) = 1 + \mu. Which pulses does \partial_\tau C contain, and where are the sources that produce them? Show that the remaining integral is small for a band-pass noise. (ch. 7)

    Solution: \bar K(1) = 2 and \bar K(-1) = 0, so only the pulse at \tau = +T survives, with twice the uniform amplitude. It comes from the sources in the direction e, beyond x_1 as seen from x_2.

    The remaining term is, with F_1 a primitive of F, \frac{1}{8\pi T} \int_{-1}^1 F(\tau - \mu T)\, d\mu = \frac{F_1(\tau + T) - F_1(\tau - T)}{8\pi T^2} .

    For a noise of central frequency \omega_0, F_1 is of order F(0)/\omega_0, so this term is smaller than the pulse by a factor of order 1/(\omega_0 T).

14. ★★ How long to record (computational). In two dimensions, put 1000 noise sources on a circle of radius 50 around two sensors 5.2 wavelengths apart. Simulate records of duration T = 100, 400, 1600 and 6400 s in the band 0.5–1.5 Hz, compute the empirical cross-correlation and its relative error, and estimate the travel time from the symmetrised envelope. (ch. 6, 7)

    Hints and expected result: Fourier coefficients of a record of length T at the frequencies m/T are nearly independent circular Gaussians (Proposition 6.10), which the code exploits. The relative error times \sqrt{T/100} stays between 7.2 and 7.7: the error decays like T^{-1/2}, as Theorem 6.8 predicts. The travel time is missed in five of the six records of 100 and 400 s, and found within 0.15 s in every record of 1600 and 6400 s.

import numpy as np

rng = np.random.default_rng(6)
c0, R, S = 1.0, 50.0, 1000                       # 2D: sources on a circle of radius R
x1, x2 = np.array([-2.6, 0.0]), np.array([2.6, 0.0])
ang = 2 * np.pi * np.arange(S) / S
Y = R * np.stack([np.cos(ang), np.sin(ang)], axis=1)
r1, r2 = np.linalg.norm(Y - x1, axis=1), np.linalg.norm(Y - x2, axis=1)
T12 = np.linalg.norm(x1 - x2) / c0
tau = np.linspace(-30, 30, 1201)

def correlate(T):
    w = 2 * np.pi * np.arange(int(0.5 * T), int(1.5 * T)) / T   # Fourier bins in the band 0.5-1.5 Hz
    u1 = np.zeros(len(w), complex); u2 = np.zeros(len(w), complex); mean = np.zeros(len(w), complex)
    for chunk in np.array_split(np.arange(S), 20):
        G1 = np.exp(1j * np.outer(r1[chunk], w) / c0) / np.sqrt(r1[chunk, None])   # far-field 2D Green
        G2 = np.exp(1j * np.outer(r2[chunk], w) / c0) / np.sqrt(r2[chunk, None])
        xi = rng.standard_normal(G1.shape) + 1j * rng.standard_normal(G1.shape)    # source spectra
        u1 += (G1 * xi).sum(0); u2 += (G2 * xi).sum(0)
        mean += 2 * (np.conj(G1) * G2).sum(0)                                     # E[conj(u1) u2]
    E = np.exp(-1j * np.outer(tau, w))
    return E @ (np.conj(u1) * u2), E @ mean       # empirical and expected correlations (analytic)

for T in [100, 400, 1600, 6400]:                 # record length (s)
    errs, est = [], []
    for _ in range(3):                           # three independent records
        C_T, C = correlate(T)
        errs.append(np.linalg.norm(C_T - C) / np.linalg.norm(C))
        sym = (np.abs(C_T) + np.abs(C_T)[::-1])[tau > 0]               # symmetric envelope
        est.append(tau[tau > 0][sym.argmax()])
    err = np.sqrt(np.mean(np.square(errs)))
    print(f"T={T:5d} s: relative error {err:.2f} (x sqrt(T/100): {err * np.sqrt(T / 100):.1f}), "
          f"travel-time estimates {np.round(est, 2)} (true {T12:.2f})")

15. ★★★ How many reflectors? Design an estimator of the number R of reflectors from the singular values of a noisy MSR matrix, by counting the values above the threshold of Proposition 5.5. Compare it with the information criteria of Wax and Kailath (1985), and study how both fail when two reflectors come closer than the Rayleigh distance. (ch. 3, 4, 5)

    Hints and expected result: with R reflectors well separated, Proposition 3.10 gives R singular values s_r, and Theorem 4.4 extends to finite rank: each s_r > \sigma\sqrt N produces an outlier. The counting estimator is consistent when all s_r exceed the threshold by a margin.

    When two reflectors approach, their illumination vectors become nearly parallel and the second singular value of A drops; expect a curve of the detection probability of the second reflector against the separation, falling near \lambda L / a.

16. ★★★ Time reversal against migration in clutter. Add to the geometry of Code 3.12 a few hundred weak random scatterers, simulated with the Born approximation or with the Foldy–Lax multiple-scattering equations. Compare, over many realisations, the time-reversal refocusing on the target and the migration image of the target. (ch. 3, 8)

    Hints and expected result: in the Born simulation the clutter echoes add a random matrix of rank equal to the number of scatterers, far from the i.i.d. model of ch. 4. With multiple scattering, time reversal should refocus with a spot no wider than in the homogeneous case and with small fluctuations between realisations, while the migration image of the target fluctuates; quantify both by the mean and variance of the peak over realisations, and relate them to the discussion of ch. 8.

Projects and assessment

The official assessment is a project, with a numerical implementation or purely theoretical, ending with a written report and an oral presentation (master-mva.com); the Université Paris-Saclay course detail describes it as partial continuous assessment, and no weights are displayed. The subjects below are sized for that project. Each names its paper or model, its deliverables and a success criterion; the theoretical ones replace code by complete proofs. The written-exam proposal is a self-test on the whole module.

Solo projects

1. [Solo · ≈ 20 h] Rectangular spiked models for two arrays. Paper study of Benaych-Georges and Nadakuditi (2012). Derive the singular value and vector limits for an N \times M matrix with N/M \to c, following the sketch of Theorem 4.4, and apply them to an MSR matrix whose emitter and receiver arrays differ. Success criterion: simulated limits within 2 % of the formulas for N = 400 and c \in \{0.25, 0.5, 1\}, and a figure of detection probability against \theta. Deliverables: code, 4–6 page report, 10-minute talk. Chapters 0, 4 and 5.
2. [Solo · ≈ 20 h] Denoising the response matrix before imaging. Paper study of Gavish and Donoho (2014). Apply the optimal hard threshold 4/\sqrt 3 \cdot \sqrt N \sigma to the noisy MSR matrices of Code 3.12, then image with migration and MUSIC, and compare with truncation at the noise edge 2\sigma\sqrt N. Success criterion: for \theta \in \{1, 1.5, 2, 3, 5\} and 100 realisations each, a figure of the median localisation error of raw, edge-truncated and threshold-truncated data, which locates the range of \theta where denoising halves the error, or shows that it never does. Deliverables: code, 4–6 page report, 10-minute talk. Chapters 3, 4 and 5.
3. [Solo · ≈ 20 h] The point-spread function of migration (theoretical). Derive the full two-dimensional point-spread function of multi-frequency Kirchhoff migration for a linear array in the paraxial regime, including an off-broadside reflector, and prove the resolution formulas of Proposition 3.7 with explicit error terms. Success criterion: complete proofs, and a numerical check against Code 3.12 within 10 %. Deliverables: 4–6 page report with proofs, short script, 10-minute talk. Chapters 1 and 3.
4. [Solo · ≈ 20 h] Selective focusing with DORT. Paper study of Prada and Fink (1994). Reproduce selective focusing on two reflectors of different strengths with the singular vectors of the time-reversal operator, then measure how the vectors mix as the reflectors approach. Success criterion: a figure of selective focusing, and a curve of the mixing angle against separation that turns over near \lambda L / a. Deliverables: code, 4–6 page report, 10-minute talk. Chapter 3.

Group projects

1. [Group of 2–3 · ≈ 30 h each] Passive travel-time tomography. In the spirit of Shapiro, Campillo, Stehly and Ritzwoller (2005), on synthetic data. Build a 2D medium with a slow anomaly and 20 stations, simulate ambient noise from sources on a surrounding circle with straight-ray travel times, extract inter-station travel times from cross-correlations, and invert them by Tikhonov regularisation. Chapters 2, 6 and 7.
   • Split: noise simulation and correlations; travel-time picking with error bars from the \sqrt{BT} rule; tomography with the discrepancy principle.
   • Success criterion: the anomaly recovered within one station spacing and with relative L^2 error below 30 %. Deliverables: code, 8-page report, 15-minute talk.
2. [Group of 2–3 · ≈ 30 h each] A detection–localisation–characterisation pipeline. Simulate multi-frequency MSR matrices of a 64-element probe facing point reflectors, with i.i.d. noise. Implement the thresholds of ch. 5, rank estimation, MUSIC and migration, and debiased reflectivity estimates. Split: simulation and detection; imaging; characterisation and ROC study. Success criterion: empirical false-alarm rate within one percentage point of the nominal 5 %, a detection curve showing the transition of Proposition 5.6, and reflectivity errors below 10 % for \theta \geq 2. Deliverables: code, 8-page report, 15-minute talk. Chapters 1, 3, 4 and 5.
3. [Group of 2–3 · ≈ 30 h each] Bayesian deconvolution with uncertainty. On the blur of Code 2.13 and on the imaging operator of ch. 3, compare three ways of setting the hyperparameters: a Gibbs sampler with conjugate inverse-gamma priors on \sigma^2 and \gamma^2, empirical Bayes by marginal likelihood, and the discrepancy principle. Split: Gibbs sampler; empirical Bayes and discrepancy; imaging operator. Success criterion: pointwise 95 % credible bands whose empirical coverage lies between 90 % and 99 % over 100 synthetic signals. Deliverables: code, 8-page report, 15-minute talk. Chapters 2 and 3.

Written exam (proposal)

Two hours, one handwritten sheet of notes allowed, four exercises covering ch. 2–6.

1. Question (5 points). State and prove the Helmholtz–Kirchhoff identity, deduce the time-reversal focal spot of a full-aperture mirror in three dimensions, and give its first zero.

   Expected: Green’s second identity applied to \hat G(\cdot, y) and \overline{\hat G(\cdot, x)}, reciprocity, the radiation condition on a large sphere (Theorem 3.5); the field proportional to \sin(k r)/(4\pi r), first zero at \lambda / 2 (Corollary 3.6).

2. Question (5 points). Let K = \mathrm{diag}(1, 0.1), \sigma = 0.1 and x \sim \mathcal{N}(0, I). Compute the posterior mean and covariance, the equivalent Tikhonov parameter, and compare the Bayes risk of the posterior mean with that of least squares.

   Expected: \alpha = \sigma^2/\gamma^2 = 0.01; posterior mean (y_1 / 1.01,\ 5\, y_2); posterior variances 0.0099 and 0.5, so \mathrm{tr}\, C = 0.51; least squares has variances 0.01 and 1, total 1.01 (Theorem 2.10, Proposition 2.12).

3. Question (5 points). An array of N = 100 sensors has noise level \sigma = 0.05. Give the threshold of the top-singular-value test at 5 %. For s_A = 1, predict the top singular value and the overlap; for s_A = 0.4, say which test has power and how much.

   Expected: threshold 0.05 \times 19.93 \approx 0.997; \theta = 2, \sigma_1 \approx 1.25, overlap 0.75; for \theta = 0.8 the top singular value is asymptotically powerless, and the energy detector has power about 1 - \Phi(1.645 - 0.64) \approx 0.16 (Propositions 5.5–5.7).

4. Question (5 points). Prove the Wiener–Khinchin theorem for an integrable covariance, and deduce that the spectral density is non-negative. Compute the spectral density of C(\tau) = e^{-|\tau|}.

   Expected: the double integral reduced to \int (T - |\tau|) C(\tau) e^{i\omega\tau} d\tau, dominated convergence (Theorem 6.3); S(\omega) = 2/(1 + \omega^2).

Grading rubric (proposal)

References

Official references (master-mva.com, “Références”)

1. Course lecture notes (“polycopié”).
2. H. Ammari et al. Mathematical and Statistical Methods for Multistatic Imaging. Lecture Notes in Mathematics 2098, Springer, 2013.
3. J. Garnier and G. Papanicolaou. Passive Imaging with Ambient Noise. Cambridge University Press, 2016.

Textbooks

4. G. W. Anderson, A. Guionnet and O. Zeitouni. An Introduction to Random Matrices. Cambridge University Press, 2010.
5. Z. Bai and J. W. Silverstein. Spectral Analysis of Large Dimensional Random Matrices, 2nd edition. Springer, 2010.
6. N. Bleistein, J. K. Cohen and J. W. Stockwell Jr. Mathematics of Multidimensional Seismic Imaging, Migration, and Inversion. Springer, 2001.
7. P. J. Brockwell and R. A. Davis. Time Series: Theory and Methods, 2nd edition. Springer, 1991.
8. D. Colton and R. Kress. Inverse Acoustic and Electromagnetic Scattering Theory, 3rd edition. Springer, 2013.
9. H. W. Engl, M. Hanke and A. Neubauer. Regularization of Inverse Problems. Kluwer, 1996.
10. J.-P. Fouque, J. Garnier, G. Papanicolaou and K. Sølna. Wave Propagation and Time Reversal in Randomly Layered Media. Springer, 2007.
11. P. C. Hansen. Rank-Deficient and Discrete Ill-Posed Problems. SIAM, 1998.
12. R. A. Horn and C. R. Johnson. Matrix Analysis, 2nd edition. Cambridge University Press, 2013.
13. J. Kaipio and E. Somersalo. Statistical and Computational Inverse Problems. Springer, 2005.
14. S. M. Kay. Fundamentals of Statistical Signal Processing, Volume II: Detection Theory. Prentice Hall, 1998.
15. E. L. Lehmann and J. P. Romano. Testing Statistical Hypotheses, 3rd edition. Springer, 2005.
16. A. M. Stuart. Inverse problems: a Bayesian perspective. Acta Numerica, 19:451–559, 2010.
17. A. Tarantola. Inverse Problem Theory and Methods for Model Parameter Estimation. SIAM, 2005.

Articles

18. J. Adler and O. Öktem. Learned primal-dual reconstruction. IEEE Transactions on Medical Imaging, 37(6):1322–1332, 2018.
19. S. Arridge, P. Maass, O. Öktem and C.-B. Schönlieb. Solving inverse problems using data-driven models. Acta Numerica, 28:1–174, 2019.
20. A. Aubry and A. Derode. Random matrix theory applied to acoustic backscattering and imaging in complex media. Physical Review Letters, 102(8):084301, 2009.
21. 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):1643–1697, 2005.
22. A. Bakulin and R. Calvert. The virtual source method: theory and case study. Geophysics, 71(4):SI139–SI150, 2006.
23. 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:120–135, 2012.
24. P. Blomgren, G. Papanicolaou and H. Zhao. Super-resolution in time-reversal acoustics. Journal of the Acoustical Society of America, 111(1):230–248, 2002.
25. L. Borcea, G. Papanicolaou and C. Tsogka. Interferometric array imaging in clutter. Inverse Problems, 21(4):1419–1460, 2005.
26. F. Brenguier, M. Campillo, C. Hadziioannou, N. M. Shapiro, R. M. Nadeau and E. Larose. Postseismic relaxation along the San Andreas fault at Parkfield from continuous seismological observations. Science, 321(5895):1478–1481, 2008.
27. M. Campillo and A. Paul. Long-range correlations in the diffuse seismic coda. Science, 299(5606):547–549, 2003.
28. E. J. Candès and C. Fernandez-Granda. Towards a mathematical theory of super-resolution. Communications on Pure and Applied Mathematics, 67(6):906–956, 2014.
29. M. Cheney. The linear sampling method and the MUSIC algorithm. Inverse Problems, 17(4):591–595, 2001.
30. H. Chung, J. Kim, M. T. McCann, M. L. Klasky and J. C. Ye. Diffusion posterior sampling for general noisy inverse problems. International Conference on Learning Representations, 2023.
31. D. L. Donoho, M. Gavish and I. M. Johnstone. Optimal shrinkage of eigenvalues in the spiked covariance model. Annals of Statistics, 46(4):1742–1778, 2018.
32. B. Engquist and B. D. Froese. Application of the Wasserstein metric to seismic signals. Communications in Mathematical Sciences, 12(5):979–988, 2014.
33. M. Fink. Time reversed acoustics. Physics Today, 50(3):34–40, 1997.
34. J. Garnier and G. Papanicolaou. Passive sensor imaging using cross correlations of noisy signals in a scattering medium. SIAM Journal on Imaging Sciences, 2(2):396–437, 2009.
35. M. Gavish and D. L. Donoho. The optimal hard threshold for singular values is 4/\sqrt 3. IEEE Transactions on Information Theory, 60(8):5040–5053, 2014.
36. M. Gavish and D. L. Donoho. Optimal shrinkage of singular values. IEEE Transactions on Information Theory, 63(4):2137–2152, 2017.
37. S. Geman. A limit theorem for the norm of random matrices. Annals of Probability, 8(2):252–261, 1980.
38. K. H. Jin, M. T. McCann, E. Froustey and M. Unser. Deep convolutional neural network for inverse problems in imaging. IEEE Transactions on Image Processing, 26(9):4509–4522, 2017.
39. K. Johansson. Shape fluctuations and random matrices. Communications in Mathematical Physics, 209(2):437–476, 2000.
40. I. M. Johnstone. On the distribution of the largest eigenvalue in principal components analysis. Annals of Statistics, 29(2):295–327, 2001.
41. W. Lambert, L. A. Cobus, T. Frappart, M. Fink and A. Aubry. Distortion matrix approach for ultrasound imaging of random scattering media. Proceedings of the National Academy of Sciences, 117(26):14645–14656, 2020.
42. V. A. Marchenko and L. A. Pastur. Distribution of eigenvalues for some sets of random matrices. Mathematics of the USSR-Sbornik, 1(4):457–483, 1967.
43. A. Onatski, M. J. Moreira and M. Hallin. Asymptotic power of sphericity tests for high-dimensional data. Annals of Statistics, 41(3):1204–1231, 2013.
44. C. Prada and M. Fink. Eigenmodes of the time reversal operator: a solution to selective focusing in multiple-target media. Wave Motion, 20(2):151–163, 1994.
45. R. O. Schmidt. Multiple emitter location and signal parameter estimation. IEEE Transactions on Antennas and Propagation, 34(3):276–280, 1986.
46. N. M. Shapiro, M. Campillo, L. Stehly and M. H. Ritzwoller. High-resolution surface-wave tomography from ambient seismic noise. Science, 307(5715):1615–1618, 2005.
47. R. Snieder. Extracting the Green’s function from the correlation of coda waves: a derivation based on stationary phase. Physical Review E, 69(4):046610, 2004.
48. A. Tarantola. Inversion of seismic reflection data in the acoustic approximation. Geophysics, 49(8):1259–1266, 1984.
49. C. A. Tracy and H. Widom. Level-spacing distributions and the Airy kernel. Communications in Mathematical Physics, 159(1):151–174, 1994.
50. J. Virieux and S. Operto. An overview of full-waveform inversion in exploration geophysics. Geophysics, 74(6):WCC1–WCC26, 2009.
51. K. Wapenaar. Retrieving the elastodynamic Green’s function of an arbitrary inhomogeneous medium by cross correlation. Physical Review Letters, 93(25):254301, 2004.
52. M. Wax and T. Kailath. Detection of signals by information theoretic criteria. IEEE Transactions on Acoustics, Speech, and Signal Processing, 33(2):387–392, 1985.
53. R. L. Weaver and O. I. Lobkis. Ultrasonics without a source: thermal fluctuation correlations at MHz frequencies. Physical Review Letters, 87(13):134301, 2001.
54. P. D. Welch. The use of fast Fourier transform for the estimation of power spectra: a method based on time averaging over short, modified periodograms. IEEE Transactions on Audio and Electroacoustics, 15(2):70–73, 1967.
55. 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):509–521, 1988.