PaperlayerView exact revision on arXiv

Source-linked AI summary

On the Numerical Evaluation of Distributions in Random Matrix Theory: A Review

Folkmar Bornemann

arXiv:0904.1581v5math.PRcs.NAmath.NA

TL;DR

Accurate numerical evaluation is needed for random-matrix distributions represented through Painlevé transcendents or Fredholm determinants. The paper compares these approaches and reports new determinant formulas for extreme levels, arguing that determinant evaluation is simpler and more readily extended to joint probabilities.

  • Problem

    Numerical evaluation is needed to obtain specific values such as moments, quantiles, or correlations from random-matrix distributions.

    Source: §1 · p. 1

  • Method

    The paper reviews and compares numerical evaluation of random-matrix distributions represented by Painlevé transcendents or Fredholm determinants.

    Source: §1 · p. 1

  • Results

    Extensive numerical exploration led to new determinantal formulas for the k-th largest or smallest levels in orthogonal and symplectic edge scaling limits.

    Source: §1 · p. 1

  • Takeaways & Limitations

    The paper argues that Fredholm determinant approximation is conceptually simpler and more efficient, and readily generalizes to joint probabilities and correlations.

    Source: §1 · p. 1

    For GOE formulas with discontinuous matrix-kernel terms, theoretical justification is difficult and unmodified numerical methods are rather inefficient.

    Source: §2.2.3 · p. 7

Abstract

from arXiv · show

In this paper we review and compare the numerical evaluation of those probability distributions in random matrix theory that are analytically represented in terms of Painlev\'e transcendents or Fredholm determinants. Concrete examples for the Gaussian and Laguerre (Wishart) beta-ensembles and their various scaling limits are discussed. We argue that the numerical approximation of Fredholm determinants is the conceptually more simple and efficient of the two approaches, easily generalized to the computation of joint probabilities and correlations. Having the means for extensive numerical explorations at hand, we discovered new and surprising determinantal formulae for the k-th largest (or smallest) level in the edge scaling limits of the Orthogonal and Symplectic Ensembles; formulae that in turn led to improved numerical evaluations. The paper comes with a toolbox of Matlab functions that facilitates further mathematical experiments by the reader.

1. Introduction

The paper addresses accurate numerical evaluation of structured random-matrix distributions, comparing Painlevé representations with Fredholm determinants. It challenges the view that Painlevé methods are numerically preferable and explores determinant methods’ scope, including joint probabilities and new formulas for edge levels.

1. Introduction

  • The paper reviews accurate numerical evaluation of RMT distributions, which matters when applications require specific values such as moments, quantiles, or correlations.The introduction frames these distributions as special functions arising from structured matrix models.Source: §1 · p. 1

1.1. The Common Point of View.

  • A prevailing view held that numerical evaluation of Fredholm determinants was quite involved, despite determinant formulas for the GUE bulk-spacing and soft-edge distributions.Earlier work evaluated the bulk-spacing determinant using prolate spheroidal wave functions, while the Tracy–Widom distribution was evaluated through a Painlevé II initial-value problem.Source: §1.1 · pp. 1–2, §1.1.1 · p. 2, §1.1.2 · p. 2

1.2. Challenging the Common Point of View.

  • The paper challenges Painlevé methods’ presumed numerical advantage by comparing simplicity, efficiency, accuracy, and extendibility.It notes that stable Painlevé computation requires asymptotic expansions of connection formulas.Source: §1.2 · p. 3
  • Fredholm determinant representations offer a simple, fast, accurate, and general numerical method, while nonlinear PDE representations have limited numerical use for multivariate functions.The paper notes that joint probabilities often have determinant representations.Source: §1.2 · p. 3
  • Extensive numerical experiments led to new determinantal formulas for k-th extreme levels and examples of joint probabilities, and the paper supplies a Matlab toolbox for further experiments.The formulas concern soft-edge GOE and GSE levels and extend by structural analogy to hard-edge LOE and LSE levels.Source: §1.3 · p. 3

2. Some Distribution Functions of RMT and Their Representation

This section surveys distribution functions for Gaussian and Laguerre ensembles and their scaling limits, alongside determinant and Painlevé representations. It also identifies a numerical complication in GOE determinant formulas and a kernel singularity for non-integer Laguerre parameters.

2. Some Distribution Functions of RMT and Their Representation

  • The section surveys fundamental RMT functions without aiming for completeness, focusing on Gaussian and Laguerre ensembles and examples useful for judging numerical methods.The ensembles considered have real spectra for Gaussian models and nonnegative spectra for Laguerre models, with β taking values 1, 2, or 4.Source: §2 · p. 4

2.1. Gaussian Ensembles.

  • Gaussian ensemble functions include bulk level-spacing probabilities and soft-edge distributions of the k-th largest eigenvalue.The bulk scaling sets the expected neighboring-eigenvalue distance to one; the soft edge concerns the largest eigenvalues.Source: §2.1.2 · p. 5, §2.1.3 · p. 5

2.2. Determinantal Representations for Gaussian Ensembles.

  • GUE determinant formulas use Hermite kernels, with sine- and Airy-kernel limits at the bulk and soft edges, respectively.The GSE soft-edge formulas include finite-rank matrix kernels, with alternative formulas for specified cases described as more efficient numerically.Source: §2.2.1 · pp. 5–6
  • GOE determinant formulas use matrix kernels with a discontinuous term, complicating theoretical justification and making unmodified numerical methods rather inefficient.For the specified special cases, the section points to alternative formulas amenable to efficient numerical evaluation.Source: §2.2.3 · p. 7

2.3. Painlevé Representations for Gaussian Ensembles.

  • For integrable kernels, determinants can be represented by PDEs in two bounds; fixing one bound gives an ODE that is a Painlevé equation in the RMT cases discussed.The section considers intervals such as (0, s) for bulk and (s, ∞) for soft-edge distributions.Source: §2.3 · pp. 7–8
  • The cited Gaussian Painlevé representations include GUE cases involving Painlevé IV, V, and II, and GOE/GSE bulk and soft-edge cases.The section attributes these representations to Tracy and Widom, Jimbo et al., and Dieng, among others.Source: §2.3.1 · p. 8, §2.3.3 · p. 9

2.4. Laguerre Ensembles.

  • The Laguerre ensembles use weights with parameter α > −1; their hard-edge limits and bulk limits are discussed, while the section confines its discussion to LUE.The Laguerre soft-edge limit matches the Gaussian result independently of α; LOE and LSE hard-edge determinant formulas are deferred to Section 7.Source: §2.4 · p. 9, §2.4.1 · p. 9

2.5–2.6. LUE Representations

  • LUE determinant formulas use Laguerre kernels and have hard-edge limits, while Painlevé representations use Painlevé V at finite size and Painlevé III at the hard edge.For non-integer α, the Laguerre and Bessel kernels have algebraic singularities at zero, a fact relevant to choosing numerical methods.Source: §2.5 · pp. 9–10, §2.6 · p. 10

3. Numerics of Painlevé Equations: the Need for Connection Formulae

Painlevé initial-value calculations can become unstable because small errors in their asymptotic starting conditions grow into large errors or blow-ups. Connection formulae expose the distinct endpoint behavior and support a more stable two-point boundary-value approach, though that approach still requires substantial numerical work.

3.1. The Straightforward Approach: Solving the Initial Value Problem.

  • The five Painlevé examples are asymptotic IVPs initialized from a one-sided endpoint expansion, but this straightforward approach unavoidably encounters instabilities.Higher terms in the asymptotic expansion can be obtained symbolically to improve the starting approximation.Source: §3.1 · p. 11

3.1.1. An example: the Tracy–Widom distribution F2.

  • The Hastings–McLeod IVP loses more than 10 digits of absolute precision and ceases to exist near x = −5.56626 when integrated leftward in double precision.The calculation starts at b− = 12; initial values are accurate to machine precision for b− ≥ 8.Source: §3.1.1 · pp. 11–12

3.1.2. Explaining the instability of the IVP.

  • In hardware arithmetic, perturbing the Hastings–McLeod asymptotic condition makes an order-one solution error unavoidable before x ≈ −12.For Painlevé II, θ < 1 and θ > 1 lead respectively to oscillatory and blow-up behavior as x → −∞.Source: §3.1.2 · pp. 12–14

3.1.3. Explaining the separation of asymptotic regimes.

  • For Painlevé II, the Fredholm determinant vanishes exactly at a pole of u(x; θ), linking blow-up locations to eigenvalues of the Airy-kernel operator.For θ > 1, the largest s satisfying θ = ρ(s)^−1/2 gives the pole location; analogous finite poles occur for z > 1 in the other listed Painlevé representations.Source: §3.1.3 · pp. 14–15

3.2. The Stable Approach: Solving the Boundary Value Problem.

  • The stable alternative solves a two-point BVP, using a connection formula to infer the asymptotic condition at one endpoint from the other.Finite boundary points and asymptotic truncation orders must be balanced for accuracy and efficiency.Source: §3.2 · p. 15

3.2.1. An Example: the Tracy–Widom distribution F2(x).

  • For F2, spectral collocation offers near-machine-precision accuracy through exponentially fast convergence, while the fixed-order bvp4c example had limited accuracy and starting-iterate difficulties.The spectral-collocation example uses a+ = −30 and b− = 8; Dieng’s bvp4c example uses a+ = −10 and b− = 6.Source: §3.2.1 · pp. 15–16

3.3. A List of Connection Formulae.

  • The connection-formula catalogue covers additional Painlevé III, IV, and V representations, including finite poles for the GUE bulk Painlevé V representation when z > 1.The listed cases include LUE hard-edge and n-dimensional LUE representations.Source: §3.3 · p. 16

3.4. Summary.

  • The BVP workflow requires locating a connection formula, balancing endpoint truncations and boundary points, and addressing solver setup; it is not yet a black-box approach.The authors say successful parameter balancing requires considerable experimentation and solver work considerable numerical-analysis experience.Source: §3.4 · pp. 16–17

4. Numerics of Fredholm Determinants and Their Derivatives

The section presents quadrature-based finite-dimensional approximations for evaluating Fredholm determinants and their derivatives. Analyticity supports fast convergence, and examples apply the method to LUE and GUE distributions.

4.1. The Basic Method.

  • An m-point quadrature rule turns the Fredholm determinant into a finite-dimensional determinant, requiring kernel evaluations and only one method parameter, m.Clenshaw–Curtis quadrature is used, with a meromorphic transformation for semi-infinite intervals.Source: §4.1 · p. 17, §4.2 · p. 17
  • For analytic kernels, the determinant approximation converges exponentially; m between 10 and 100 typically yields machine precision, and evaluating F2(s) takes a few milliseconds.Source: §4.1 · p. 17

4.2. Numerical Evaluation of Finite-Dimensional Determinants.

  • A determinant at one z is computed by LU decomposition in O(m^3); eigenvalues instead allow subsequent evaluations of dm(z) in O(m) each.Source: §4.2 · pp. 17–18

4.3. Numerical Evaluation of Higher Derivatives.

  • Cauchy integrals and the trapezoidal rule provide exponentially convergent higher derivatives, while the integration radius must be chosen carefully for numerical stability.The authors recommend κ ≈ 1, with radii r = 1 for bulk and r = 0.1 for edge limits.Source: §4.3 · p. 18

4.4. Error Control.

  • The difference between successive contracted approximations estimates the error; for F2(−2), Table 3 shows 15 correct digits at m = 64 after machine-precision saturation.Source: §4.4 · p. 19, Table 3 · p. 19

4.5. Numerical Evaluation of Densities.

  • Chebyshev interpolation approximates real-analytic CDFs stably, enabling calculation of low-order derivatives such as densities and integrals such as moments.Source: §4.5 · p. 19

4.6.1. Distribution of smallest and largest level in a specific LUE.

  • For LUE with n = 80 and α = 40, Table 4 reports moments of distributions (4.15)–(4.18), showing ten correctly truncated digits that passed error control; computations used hardware arithmetic.Source: Table 4 · p. 21

4.6.2. The distribution of k-level spacings in the bulk of GUE.

  • GUE bulk spacing functions for k = 0–14 and soft-edge densities for k = 1–6 are computed to about 5 · 10^-15 absolute error, with each plot taking about 30 seconds.Source: Figure 4 · p. 21

5. The Distribution of k-Level Spacings in the Bulk: GOE and GSE

The section uses a factorization of the sine-kernel determinant to express bulk spacing functions through the functions E+ and E−. These representations extend to GOE and GSE and support numerical plots of their spacing functions.

5. The Distribution of k-Level Spacings in the Bulk: GOE and GSE

  • Mehta’s formulae express the GOE and GSE bulk k-level spacing functions through E+ and E− and are directly amenable to the section’s numerical approach.Source: §5 · pp. 22–23
  • Shift invariance of the sine-kernel convolution operator yields an even–odd block factorization, which is used with the Leibniz formula to represent GUE spacing functions through E+ and E−.Source: §5 · pp. 22–23
  • The numerical methods reproduce plots of the GOE and GSE bulk spacing functions using determinantal formulae expressed through E+ and E−.Source: §5 · p. 23

6. The k-th Largest Level at the Soft Edge: GOE and GSE

The section derives new determinantal formulas for the soft-edge cumulative distributions of the k-th largest GOE and GSE levels. It also gives recursive GOE formulas and uses these results for numerical evaluation.

6. The k-th Largest Level at the Soft Edge: GOE and GSE

  • New determinantal formulas express the soft-edge cumulative distributions F1(k; s) and F4(k; s) through Theorem 6.1’s identity, valid for all real s and z in an analytic domain containing 1.The authors first checked the identity numerically at 100,000 randomly chosen arguments, then proved it; the resulting formulas are amenable to efficient numerical evaluation.Source: §6 · pp. 23–25
  • F4(k; s) = F1(2k; s) for k = 1, 2, 3, …, so the GSE soft-edge distributions follow from the even-indexed GOE distributions.The relation follows from the interrelationships between the ensembles in the soft-edge scaling limit.Source: §6 · p. 26
  • The GOE functions ˜E1(k; s) have recursively defined formulas obtained from a functional-equation system, with separate even- and odd-indexed cases.Theorem 6.2 establishes the system’s unique recursively defined solution.Source: §6 · pp. 26–27
  • Figure 6 plots GOE soft-edge densities for k = 1, …, 6 and the absolute errors of Painlevé II evaluations for k = 1, …, 4.The plotted calculations use hardware arithmetic and have absolute error about 5 · 10−15.Source: Figure 6 · p. 28

Remark

The determinantal formula provides a route to a Painlevé representation, and the section reports numerical evaluations of GOE soft-edge distributions alongside error comparisons.

Remark

  • The determinantal formula gives an elementary route to the Painlevé representation of D1(z; (s, ∞)).The representation is stated for s > 0 and z in an analytic domain containing 1.Source: Remark · p. 28
  • GOE soft-edge distributions for k = 1, …, 6 were evaluated from the recursion formulas with absolute tolerance 5 · 10−15.Because of the relation in the preceding section, separate GSE calculations were not needed.Source: Remark · pp. 28–29
  • Figure 6.b compares absolute errors for Dieng’s Painlevé II evaluations at k = 1, …, 4 and shows switches at a+ = −10 and b− = 6.The authors use their numerical solutions as the reference.Source: Remark · p. 29

7. The k-th Smallest Level at the Hard Edge: LOE and LSE

The section derives determinantal formulas for the hard-edge distributions of the k-th smallest levels in the Laguerre ensembles. Their algebraic structure parallels the soft-edge formulas, and a hard-to-soft transition is established.

7. The k-th Smallest Level at the Hard Edge: LOE and LSE

  • Determinantal formulas are derived for the k-th smallest level distributions of LOE, LUE, and LSE with parameter α, with the same algebraic structure as the soft-edge formulas.The derivation proceeds in parallel to the preceding section.Source: §7 · p. 29
  • The hard-edge functional-equation system has the same algebraic structure as the one solved for the soft edge, yielding a recursive solution.The section gives relations between the even- and odd-indexed LOE distributions and E+,α and E−,α.Source: §7 · pp. 30–31
  • As α tends to infinity under the stated scaling, E±,α(k; s) converge to the corresponding soft-edge functions, yielding a hard-to-soft transition.The paper identifies this transition as a means of validating its numerical methods.Source: §7 · p. 31

8. Matrix Kernels and Examples of Joint Probability Distributions

The section extends quadrature-based Fredholm determinant evaluation to matrix kernels and joint probabilities. Its GUE examples show statistically dependent top levels and a visible difference between their summed-level CDF and the convolution of individual CDFs.

8.1. Matrix Kernels.

  • An m-point quadrature rule on each of N intervals approximates the matrix-kernel determinant by an N·m-dimensional determinant.The intervals are open and the smooth matrix kernel generates a trace-class operator on the direct sum of their L2 spaces.Source: §8.1 · p. 32
  • Analytic matrix-kernel entries yield exponential convergence for the quadrature approximation, and the methods of Section 4.2–4.5 apply locally uniformly in z.The approach was used for numerical checks of equation (6.10).Source: §8.1 · p. 32

8.1.1. An example: the joint probability distribution of GUE matrix diffusion.

  • The two-time GUE matrix-diffusion covariance is computed from a joint probability distribution and checked against Monte Carlo data and small- and large-time asymptotic expansions.The figure reports 10-digit accuracy with absolute tolerance 5 · 10^-11.Source: Figure 7 · p. 33, §8.1.1 · p. 33

8.2. Operators Acting on Unions of Intervals.

  • Operators on a union of disjoint intervals can be transformed into matrix-kernel determinants, including weighted versions, and then evaluated by the Section 8.1 method.The operators in the blocks share one kernel function.Source: §8.2 · pp. 34–35

8.3. Generalized Spacing Functions.

  • Generalized spacing functions give the probability of specified level counts across disjoint intervals, but differentiating in multiple variables can make their numerical evaluation expensive.The resulting computation requires a multidimensional Cauchy integral of a function calculated from an approximated Fredholm determinant.Source: §8.3 · p. 35

8.3.1. Efficient numerical evaluation of the finite-dimensional determinants.

  • For the special multi-index α = (k, 0, . . . , 0), a factorization reduces the finite-dimensional determinant to a form handled by the Section 4.2 results.The example of the top two GUE levels pushes the numerical methods close to their limit; the reported results are out of reach of PDE representations.Source: §8.3.1 · p. 35, §8.3.2 · p. 36

8.3.2. An example: the joint distribution of the largest two eigenvalues in GUE.

  • The CDF of the sum of the two largest GUE soft-edge levels visibly differs from the convolution of their individual CDFs, reflecting statistical dependence.The joint-probability calculation used about 16 hours with hardware arithmetic, and its 11 displayed correlation digits were estimated correct.Source: Figure 7 · p. 33, §8.3.2 · p. 36

9. Software

The software section presents a Matlab toolbox for evaluating the paper’s probability distributions, controlling numerical error, and computing densities, moments, and quantiles. Examples check its outputs against known values and distribution constraints.

9.1.1. Quadrature rule.

  • The low-level interface offers Clenshaw–Curtis and Gauss–Jacobi quadrature, while medium-level commands conceal point counts and return determinant values with error estimates.Kernel construction supports grid and outer vectorization modes.Source: §9.1.1 · p. 37, §9.1.2 · p. 38, §9.2 · p. 39

9.2.1. Example 9.1.4 revisited.

  • Automatic error control reports 14 correct digits for both GOE and GUE Tracy–Widom values at zero, matching Prähofer’s tables.The earlier low-level example with 64 quadrature points was correct to one unit of the last decimal place.Source: §9.1.4 · p. 38, §9.2.1 · p. 39

9.2.2. Example: an instance of equation (6.10).

  • An instance of the determinantal equation (6.10) agrees within the default tolerance of 5 · 10^-15.The example compares two computed values for one real s and one complex z.Source: §9.2.2 · p. 40

9.3. High Level Commands.

  • The toolbox provides high-level commands for the paper’s probability distributions, returning values and conservative error estimates.The software also supplies low- and medium-level commands for quadrature, kernel construction, and determinant approximation.Source: §9 · p. 37, §9.3 · p. 40, Table 6 · p. 41

9.3.2. Example: checking the k-level spacing functions against a constraint.

  • The toolbox’s k-level spacing probabilities reproduce the mass and mean constraints in the reported checks.The examples also check generalized probabilities against their stated constraints, with results reported as correctly rounded.Source: §9.3.2 · pp. 40, 42, §9.3.3 · pp. 42–43

9.3.4. Example: calculating quantiles.

  • The computed 5% and 95% GOE Tracy–Widom quantiles are -3.18037997693773 and 0.979316053469556.They are obtained by solving F(1,s) = p for the specified probabilities.Source: §9.3.4 · p. 43

9.4. Densities and Moments.

  • The toolbox computes probability densities and the first four moments by barycentric interpolation, and can automatically choose the interpolation-point count.The returned outputs include error estimates, numerical support, and interpolants for both the PDF and CDF.Source: §9.4 · p. 43

9.4.1. Example: the first four moments of the Tracy–Widom distributions for GOE, GUE, and GSE.

  • For GOE, GUE, and GSE, the plotted Tracy–Widom PDFs and CDFs place larger β distributions farther left.The example fills in missing high-precision GSE moment digits in Tracy and Widom (2008, Table 1).Source: §9.4.1 · pp. 43–44

A.1. Algebraic Kernel Singularities and Gauss–Jacobi Quadrature.

The standard Clenshaw–Curtis approach does not generally converge exponentially for kernels lacking analyticity near the integration interval.

A.1. Algebraic Kernel Singularities and Gauss–Jacobi Quadrature.

  • For nonanalytic kernels, straightforward Clenshaw–Curtis quadrature does not generally yield exponential convergence.The stated examples include Laguerre and Bessel kernels with non-integer parameter, which have algebraic singularities at the boundary.Source: §A.1 · p. 45
  • With the software’s default maximum of 256 quadrature points, this loss of convergence can produce inaccurate results.The automatic error control detects these inaccuracies.Source: §A.1 · p. 45
  • The stated limitation applies in particular to Laguerre and Bessel kernels with non-integer parameter and algebraic endpoint singularities.Source: §A.1 · p. 45

A.1.1. Example: level spacing at the hard edge for non-integer parameter.

For hard-edge Bessel-kernel calculations with non-integer parameters, Gauss–Jacobi quadrature addresses algebraic endpoint singularities and improves the reported accuracy over Clenshaw–Curtis quadrature.

A.1.1. Example: level spacing at the hard edge for non-integer parameter.

  • For the stronger singularity α = −1, Clenshaw–Curtis quadrature gives just 1 digit in the hard-edge example.Source: §A.1.1 · p. 45

A.1.2. Gauss–Jacobi quadrature.

  • Gauss–Jacobi quadrature addresses algebraic endpoint singularities; the example reports 14 digits for both non-integer parameter cases.The singular exponent for the determinant approximation is governed by the singularity of Kα(x, y)uλ(y) as y approaches zero.Source: §A.1.2 · pp. 45–46

A.2. Tables of Some Statistical Properties.

  • The appendix tables report statistical properties of k-level spacing densities and soft-edge k-th-largest-level distributions, with correctly truncated digits passing automatic error control.Computing times are reported in seconds; GOE–GSE interrelations remove the need to tabulate β = 4 separately.Source: §A.2 · p. 47
Loading 0904.1581v5…