From d30ca8379f15d87b1f716a344145080d2f099083 Mon Sep 17 00:00:00 2001 From: Austen Lamacraft Date: Thu, 24 May 2018 12:41:03 +0200 Subject: [PATCH 1/6] Added Bornemann method --- src/densities/TracyWidom.jl | 106 ++++++++++++++++++++++-------------- 1 file changed, 66 insertions(+), 40 deletions(-) diff --git a/src/densities/TracyWidom.jl b/src/densities/TracyWidom.jl index 46f7f90..601ae31 100644 --- a/src/densities/TracyWidom.jl +++ b/src/densities/TracyWidom.jl @@ -4,18 +4,21 @@ export TracyWidom Tracy-Widom distribution The probability distribution of the normalized largest eigenvalue of a random -Hermitian matrix. +matrix with iid Gaussian matrix elements of variance 1/2 and mean 0. The cdf of Tracy-Widom is given by `` -F_2 (s) = lim_{n→∞} Pr(√2 n^{1/6} (λₙ - √(2n) ≤ s) +F_beta (s) = lim_{n→∞} Pr(√2 n^{1/6} (λ_max - √(2n) ≤ s) `` +where beta = 1, 2, or 4 for the orthogonal, unitary, or symplectic ensembles. + References: 1. doi:10.1016/0370-2693(93)91114-3 2. doi:10.1007/BF02100489 +3. doi.org/10.1007/BF02099545 Numerical routines adapted from Alan Edelman's course notes for MIT 18.338, Random Matrix Theory, 2016. @@ -24,59 +27,82 @@ struct TracyWidom <: ContinuousUnivariateDistribution end """ -Probability density function of the Tracy-Widom distribution +Cumulative density function of the Tracy-Widom distribution -Computes the Tracy-Widom distribution by directly solving the -Painlevé II equation using the Vern8 numerical integrator +Computes the Tracy-Widom distribution by Bornemann's method of evaluating +a finite dimensional approximation to the Fredholm determinant using quadrature. + +doi.org/10.1090/S0025-5718-09-02280-7 # Arguments * `d::TracyWidom` or `Type{TracyWidom}`: an instance of `TracyWidom` or the type itself -* `t::Real`: The point at which to evaluate the pdf -* `t0::Real = -8.0`: The point at which to start integrating +* `s::Real`: The point at which to evaluate the cdf +* `beta::Integer = 2`: The Dyson index defining the distribution. Takes values 1, 2, or 4. +* `num_points::Integer = 25`: The number of points in the quadrature """ -function pdf(d::TracyWidom, t::S, t0::S = convert(S, -8.0)) where {S<:Real} - t≤t0 && return 0.0 - t≥5 && return 0.0 +function cdf{T<:Real}(d::TracyWidom, s::T; beta::Integer=2, num_points::Integer=25) + beta ∈ (1,2,4) || throw(ArgumentError("Beta must be 1, 2, or 4")) + quad = gausslegendre(num_points) + _TWcdf(s, beta, quad) +end - sol = _solve_painleve_ii(t0, t) +function cdf{T<: Real}(d::Type{TracyWidom}, s::T; beta::Integer=2, num_points::Integer=25) + cdf(d(), s, beta=beta, num_points=num_points) +end - ΔF2=exp(-sol[end-1][3]) - exp(-sol[end][3]) # the cumulative distribution - f2=ΔF2/(sol.t[end]-sol.t[end-1]) # the density at t +function cdf{T<:Real}(d::TracyWidom, s_vals::AbstractArray{T}; beta::Integer=2, num_points::Integer=25) + beta ∈ (1,2,4) || throw(ArgumentError("Beta must be 1, 2, or 4")) + quad = gausslegendre(num_points) + [_TWcdf(s, beta, quad) for s in s_vals] end -pdf(d::Type{TracyWidom}, t::Real) = pdf(d(), t) -""" -Cumulative density function of the Tracy-Widom distribution +function cdf{T<:Real}(d::Type{TracyWidom}, s_vals::AbstractArray{T}; beta::Integer=2, num_points::Integer=25) + cdf(d(), s_vals, beta=beta, num_points=num_points) +end -Computes the Tracy-Widom distribution by directly solving the -Painlevé II equation using the Vern8 numerical integrator +function _TWcdf{T<:Real}(s::T, beta::Integer, quad::Tuple{Array{T,1},Array{T,1}}) + if beta == 2 + kernel = ((ξ,η) -> _K2tilde(ξ,η,s)) + return _fredholm_det(kernel, quad) + elseif beta == 1 + kernel = ((ξ,η) -> _K1tilde(ξ,η,s)) + return _fredholm_det(kernel, quad) + elseif beta == 4 + kernel2 = ((ξ,η) -> _K2tilde(ξ,η,s*sqrt(2))) + kernel1 = ((ξ,η) -> _K1tilde(ξ,η,s*sqrt(2))) + F2 = _fredholm_det(kernel2, quad) + F1 = _fredholm_det(kernel1, quad) + return (F1 + F2/F1) / 2 + end +end -See `pdf(::TracyWidom)` for a description of the arguments. -""" -function cdf{S<:Real}(d::TracyWidom, t::S, t0::S = convert(S, -8.0)) - t≤t0 && return 0.0 - t≥5 && return 1.0 - sol = _solve_painleve_ii(t0, t) - F2=exp(-sol[end][3]) +function _fredholm_det{T<:Real}(kernel::Function, quad::Tuple{Array{T,1},Array{T,1}}) + nodes, weights = quad + N = length(nodes) + sqrt_weights = sqrt.(weights) + weights_matrix = kron(transpose(sqrt_weights),sqrt_weights) + K_matrix = [kernel(ξ,η) for ξ in nodes, η in nodes] + det(eye(N) - weights_matrix .* K_matrix) end -cdf(d::Type{TracyWidom}, t::Real) = cdf(d(), t) - -# An internal function which sets up the Painleve II differential equation and -# runs it through the Vern8 numerical integrator -function _solve_painleve_ii{S<:Real}(t0::S, t::S) - function deq(dy, y, p, t) - dy[1] = y[2] - dy[2] = t*y[1]+2y[1]^3 - dy[3] = y[4] - dy[4] = y[1]^2 + +_ϕ(ξ, s) = s + 10*tan(π*(ξ+1)/4) +_ϕprime(ξ) = (5π/2)*(sec(π*(ξ+1)/4))^2 + +# For beta = 2 +function _airy_kernel(x, y) + if x==y + return (airyaiprime(x))^2 - x * (airyai(x))^2 + else + return (airyai(x) * airyaiprime(y) - airyai(y) * airyaiprime(x)) / (x - y) end - a0 = airyai(t0) - T = typeof(big(a0)) - y0=T[a0, airyaiprime(t0), 0, airyai(t0)^2] # Initial conditions - prob = ODEProblem(deq,y0,(t0,t)) - solve(prob, Vern8(), abstol=1e-12, reltol=1e-12) # Solve the ODE end +_K2tilde(ξ,η,s) = sqrt(_ϕprime(ξ) * _ϕprime(η)) * _airy_kernel(_ϕ(ξ,s), _ϕ(η,s)) + +# For beta = 1 +_A_kernel(x,y) = airyai((x+y)/2) / 2 +_K1tilde(ξ,η,s) = sqrt(_ϕprime(ξ) * _ϕprime(η)) * _A_kernel(_ϕ(ξ,s), _ϕ(η,s)) + """ Samples the largest eigenvalue of the n × n GUE matrix """ From 449d826064832848747a53eb8ccead46fd52dcd0 Mon Sep 17 00:00:00 2001 From: Austen Lamacraft Date: Thu, 24 May 2018 16:18:05 +0200 Subject: [PATCH 2/6] Added tests --- src/RandomMatrices.jl | 1 + src/densities/TracyWidom.jl | 2 +- test/densities/TracyWidom.jl | 17 ++++++++--------- 3 files changed, 10 insertions(+), 10 deletions(-) diff --git a/src/RandomMatrices.jl b/src/RandomMatrices.jl index 179ff59..0fc1bcb 100644 --- a/src/RandomMatrices.jl +++ b/src/RandomMatrices.jl @@ -4,6 +4,7 @@ using Combinatorics using GSL using OrdinaryDiffEq using DiffEqBase # This line is only needed on v0.5, and comes free from OrdinaryDiffEq on v0.6 +using SpecialFunctions, FastGaussQuadrature import Base: isinf, rand, convert import Distributions: ContinuousUnivariateDistribution, diff --git a/src/densities/TracyWidom.jl b/src/densities/TracyWidom.jl index 601ae31..dcfa0cd 100644 --- a/src/densities/TracyWidom.jl +++ b/src/densities/TracyWidom.jl @@ -37,7 +37,7 @@ doi.org/10.1090/S0025-5718-09-02280-7 # Arguments * `d::TracyWidom` or `Type{TracyWidom}`: an instance of `TracyWidom` or the type itself * `s::Real`: The point at which to evaluate the cdf -* `beta::Integer = 2`: The Dyson index defining the distribution. Takes values 1, 2, or 4. +* `beta::Integer = 2`: The Dyson index defining the distribution. Takes values 1, 2, or 4 * `num_points::Integer = 25`: The number of points in the quadrature """ function cdf{T<:Real}(d::TracyWidom, s::T; beta::Integer=2, num_points::Integer=25) diff --git a/test/densities/TracyWidom.jl b/test/densities/TracyWidom.jl index a7b4799..8688b6e 100644 --- a/test/densities/TracyWidom.jl +++ b/test/densities/TracyWidom.jl @@ -2,16 +2,15 @@ using RandomMatrices using Compat.Test #Test far outside support -#Tracy-Widom has support on all x>0, but the integration won't pick it up -@test pdf(TracyWidom, -10) == pdf(TracyWidom, 10) == 0 -@test cdf(TracyWidom, -10) == 0 -@test cdf(TracyWidom, 10) == 1 +@test cdf(TracyWidom, -10) ≈ 0 atol=1e-14 +@test cdf(TracyWidom, 10) ≈ 1 atol=1e-14 -if isdefined(:OrdinaryDiffEq) && isa(OrdinaryDiffEq, Module) - t = rand() - @test pdf(TracyWidom, t) > 0 - @test 0 < cdf(TracyWidom, t) < 1 -end +# Test exact values +# See https://arxiv.org/abs/0904.1581 +@test TWcdf(0,beta=1) ≈ 0.83190806620295 atol=1e-14 +@test TWcdf(0,beta=2) ≈ 0.96937282835526 atol=1e-14 + +@test 0 < cdf(TracyWidom, rand()) < 1 @test isfinite(rand(TracyWidom, 10)) @test isfinite(rand(TracyWidom, 100)) From 58215f494323c34de399b4dbafd7200d6938c084 Mon Sep 17 00:00:00 2001 From: Austen Lamacraft Date: Thu, 24 May 2018 16:40:13 +0200 Subject: [PATCH 3/6] tests passing --- src/densities/TracyWidom.jl | 2 +- test/densities/TracyWidom.jl | 13 +++++++------ 2 files changed, 8 insertions(+), 7 deletions(-) diff --git a/src/densities/TracyWidom.jl b/src/densities/TracyWidom.jl index dcfa0cd..f6029f6 100644 --- a/src/densities/TracyWidom.jl +++ b/src/densities/TracyWidom.jl @@ -60,7 +60,7 @@ function cdf{T<:Real}(d::Type{TracyWidom}, s_vals::AbstractArray{T}; beta::Integ cdf(d(), s_vals, beta=beta, num_points=num_points) end -function _TWcdf{T<:Real}(s::T, beta::Integer, quad::Tuple{Array{T,1},Array{T,1}}) +function _TWcdf{T<:Real}(s::T, beta::Integer, quad::Tuple{Array{Float64,1},Array{Float64,1}}) if beta == 2 kernel = ((ξ,η) -> _K2tilde(ξ,η,s)) return _fredholm_det(kernel, quad) diff --git a/test/densities/TracyWidom.jl b/test/densities/TracyWidom.jl index 8688b6e..64884ba 100644 --- a/test/densities/TracyWidom.jl +++ b/test/densities/TracyWidom.jl @@ -1,16 +1,17 @@ using RandomMatrices using Compat.Test +# CDF lies in correct range +@test all(i->(0 Date: Thu, 24 May 2018 16:42:20 +0200 Subject: [PATCH 4/6] Adjust docstring --- src/densities/TracyWidom.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/densities/TracyWidom.jl b/src/densities/TracyWidom.jl index f6029f6..5338dbb 100644 --- a/src/densities/TracyWidom.jl +++ b/src/densities/TracyWidom.jl @@ -27,7 +27,7 @@ struct TracyWidom <: ContinuousUnivariateDistribution end """ -Cumulative density function of the Tracy-Widom distribution +Cumulative density function of the Tracy-Widom distribution. Computes the Tracy-Widom distribution by Bornemann's method of evaluating a finite dimensional approximation to the Fredholm determinant using quadrature. From e67d9235dd120a2779774a7fd43f570830f69389 Mon Sep 17 00:00:00 2001 From: Austen Lamacraft Date: Thu, 24 May 2018 17:12:34 +0200 Subject: [PATCH 5/6] Updated REQUIRE --- REQUIRE | 2 ++ 1 file changed, 2 insertions(+) diff --git a/REQUIRE b/REQUIRE index 54fc222..1112328 100644 --- a/REQUIRE +++ b/REQUIRE @@ -4,3 +4,5 @@ Distributions 0.8.10 OrdinaryDiffEq 3.0.0 GSL 0.3.1 Compat 0.60 +SpecialFunctions 0.4.0 +FastGaussQuadrature 0.3.0 From 8f2ab443b6b1c89ee40be7a3e6e6b11d4ffbd178 Mon Sep 17 00:00:00 2001 From: Austen Lamacraft Date: Fri, 25 May 2018 23:23:12 +0100 Subject: [PATCH 6/6] Removed dependency on OrdinaryDiffEq --- REQUIRE | 1 - src/RandomMatrices.jl | 2 -- 2 files changed, 3 deletions(-) diff --git a/REQUIRE b/REQUIRE index 1112328..989e22c 100644 --- a/REQUIRE +++ b/REQUIRE @@ -1,7 +1,6 @@ julia 0.6 Combinatorics 0.2 Distributions 0.8.10 -OrdinaryDiffEq 3.0.0 GSL 0.3.1 Compat 0.60 SpecialFunctions 0.4.0 diff --git a/src/RandomMatrices.jl b/src/RandomMatrices.jl index 0fc1bcb..fd022ba 100644 --- a/src/RandomMatrices.jl +++ b/src/RandomMatrices.jl @@ -2,8 +2,6 @@ module RandomMatrices using Combinatorics using GSL -using OrdinaryDiffEq -using DiffEqBase # This line is only needed on v0.5, and comes free from OrdinaryDiffEq on v0.6 using SpecialFunctions, FastGaussQuadrature import Base: isinf, rand, convert