Implementation of Bornemann's method for Tracy-Widom distributions (#44)

* Added Bornemann method

* Added tests

* tests passing

* Adjust docstring

* Updated REQUIRE

* Removed dependency on OrdinaryDiffEq

Author: AustenLamacraft

REQUIRE

@@ -1,6 +1,7 @@
 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
+FastGaussQuadrature 0.3.0

src/RandomMatrices.jl

@@ -2,8 +2,7 @@ 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
 import Distributions: ContinuousUnivariateDistribution,

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 Bornemann's method of evaluating
+a finite dimensional approximation to the Fredholm determinant using quadrature.
 
-Computes the Tracy-Widom distribution by directly solving the
-Painlevé II equation using the Vern8 numerical integrator
+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{Float64,1},Array{Float64,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
 """

test/densities/TracyWidom.jl

@@ -1,17 +1,17 @@
 using RandomMatrices
 using Compat.Test
 
+# CDF lies in correct range
+@test all(i->(0<i<1), cdf(TracyWidom, randn(5)))
+
 #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.1) ≈ 0 atol=1e-14
+@test cdf(TracyWidom, 10.1) ≈ 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 cdf(TracyWidom,0,beta=1) ≈ 0.83190806620295 atol=1e-14
+@test cdf(TracyWidom,0,beta=2) ≈ 0.96937282835526 atol=1e-14
 
 @test isfinite(rand(TracyWidom, 10))
 @test isfinite(rand(TracyWidom, 100))