Skip to content

Commit dcea45b

Browse files
committed
chapter 6 of rmt (5.1 to 5.5)
1 parent d74fbab commit dcea45b

File tree

8 files changed

+151
-1
lines changed

8 files changed

+151
-1
lines changed

demos/book/1/semicircle.jl

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -4,6 +4,7 @@
44
#Experiment: Sample random symmetric Gaussian matrices
55
#Plot: Histogram of the eigenvalues
66
#Theory: Semicircle as n->infinity
7+
#See also: RandomMatrices.SemiCircle Distribution
78

89
## Parameters
910
n=1000 # matrix size

demos/book/1/tracywidom.jl

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -2,6 +2,7 @@
22
#Algorithm 1.3 of Random Eigenvalues by Alan Edelman
33

44
#Theory: Compute and plot the Tracy-Widom distribution
5+
#See also: RandomMatrices.TracyWidom(t) for tn=-8
56

67
##Parameters
78
t0=5. # right endpoint

demos/book/4/cauchymp.jl

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -1,4 +1,4 @@
1-
#achchymp.jl
1+
#cauchymp.jl
22
#Algorithm 4.1 of Random Eigenvalues by Alan Edelman
33

44
#Experiment: Calculate Cauchy transform as trace of resolvent of Wishart matrix

demos/book/5/compute_hermite.jl

Lines changed: 26 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,26 @@
1+
include("orthopoly_evaluate.jl")
2+
function compute_hermite(n, x, all=true);
3+
# Evaluate Orthonormalized Hermite polynomials at x
4+
# Physicists Hermite:
5+
# w(x)=exp(-x^2) with integral sqrt(pi)
6+
# Equivalent Mathematica: HermiteH[n,x]/Sqrt[Sqrt[Pi]*n!*2^n]
7+
# Probabilists Hermite:
8+
# w(x)=exp(-x^2/2) with integral sqrt(2*pi)
9+
# Equivalent Mathematica: HermiteH[n,x/Sqrt[2]]/Sqrt[Sqrt[Pi]*n!*2^(1/2+n)]
10+
11+
## Physicists Hermite Tridiagonal matrix (w=exp(-x^2))
12+
T = diagm(sqrt((1:n)/2),1) # n+1 by n+1 matrix
13+
T = T+T'
14+
c = sqrt(pi) #Integral of weight function
15+
## Probabilists Hermite Tridiagonal matrix (w=exp(-x^2/2))
16+
# T = diagm(sqrt((1:n)),1) # n+1 by n+1 matrix
17+
# T = T+T'
18+
# c = sqrt(2*pi)
19+
## Default to 0 through n computation
20+
if all
21+
phi = orthopoly_evaluate_all(T, x)
22+
else
23+
phi = orthopoly_evaluate_n(T, x)
24+
end
25+
phi/sqrt(c)
26+
end

demos/book/5/compute_laguerre.jl

Lines changed: 20 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,20 @@
1+
## This function computes the Generalized Laguerre polynomials using the Tridiagonal matrix
2+
include("orthopoly_evaluate.jl")
3+
function compute_laguerre(n, a, x, all=true)
4+
# Evaluate Orthonormalized Laguerre polynomials at x
5+
# Laguerre:
6+
# w(x)=x^a exp(-x) with integral gamma(a+1) (a>-1)
7+
# Equivalent Mathematica:LaguerreL[n,a,x]/Sqrt[Gamma[n+a+1]/Gamma[n+1]]
8+
## Laguerre Bidiagonal matrix
9+
B = diagm(sqrt((a + 1 +(0:n)))) # diagonal
10+
B = B - diagm(sqrt((1:n)), 1) # superdiagonal
11+
T = B'*B
12+
c = gamma(a+1) #Integral of weight function
13+
## Default to 0 through n computation
14+
if all
15+
phi = orthopoly_evaluate_all(T, x)
16+
else
17+
phi = orthopoly_evaluate_n(T, x)
18+
end
19+
phi/sqrt(c)
20+
end

demos/book/5/hermitekernel.jl

Lines changed: 18 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,18 @@
1+
include("compute_hermite.jl")
2+
function hermitekernel(n,x)
3+
#Compute the Hermite kernel matrix at points x
4+
#including degrees from 0 to n inclusive
5+
6+
# # Method 1: Direct Computation
7+
phi=compute_hermite(n,x);
8+
K=phi'*phi;
9+
10+
# # Method 2: Christoffel-Darboux
11+
# phinp1 = compute_hermite(n+1,x,false)
12+
# phin = compute_hermite(n,x,false)
13+
# K = phinp1'*phin*sqrt((n+1)/2)
14+
# x = x[:]*ones(1,length(x))
15+
# K = (K-K')./(x-x')
16+
17+
K
18+
end

demos/book/5/orthopoly_evaluate.jl

Lines changed: 55 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,55 @@
1+
function orthopoly_evaluate_all(T, x)
2+
## Evaluate 0th through nth orthogonal polynomial at x
3+
# Parameters
4+
# T: Tridiagonal matrix encoding the three term recurrence
5+
# size is (n+1)x(n+1)
6+
# x: evaluation points
7+
#
8+
# Returns
9+
# phi: Polynomials 0 through n at x
10+
# Method: Spectrally solve (T-XI)u = constant * (last column of I)
11+
# WARNING: potential for dividing by 0
12+
# 1st row of phi is 1 -- for w not a probability measure
13+
# divide by 1/sqrt(c), where c is the total weight
14+
# of the real line
15+
n = size(T,1)
16+
Lambda, H = eig(T)
17+
Hn = H[end, :]'
18+
phi = Array(Float64, length(Lambda), length(x))
19+
for m = 1:length(x)
20+
xi = x[m]
21+
# generate the u vector as in (2.5)
22+
u = H * (Hn./(Lambda - xi))
23+
phi[:,m] = u/u[1]
24+
end
25+
phi
26+
end
27+
28+
function orthopoly_evaluate_n(T, x)
29+
## Evaluate nth orthogonal polynomial at x
30+
# Parameters
31+
# T: Tridiagonal matrix encoding the three term recurrence
32+
# size is (n+1)x(n+1)
33+
# x: evaluation points
34+
#
35+
# Returns
36+
# phi: nth at x
37+
# normalized so that 0'th is 1 (w has unit weight)
38+
# Method: characteristic polynomial of n x n T
39+
# with leading coefficient taken from the
40+
# product of the superdiagonal of the
41+
# (n+1)x(n+1) T
42+
43+
n = size(T,1)-1
44+
## compute the characteristic polynomial of T[1:n,1:n]
45+
Lambda = eigvals(T[1:n, 1:n])
46+
phi = ones(size(x))
47+
for i = 1:n
48+
phi = phi.*(x - Lambda[i])
49+
end
50+
## normalization
51+
if n > 0
52+
phi = phi / prod(diag(T,1))
53+
end
54+
phi
55+
end

demos/book/5/sinekernel.jl

Lines changed: 29 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,29 @@
1+
#sinekernel.jl
2+
3+
n=200
4+
x=(0:.025:1)*3*pi
5+
include("hermitekernel.jl")
6+
K=hermitekernel(n-1,x/sqrt(2*n))*pi/sqrt(2*n)
7+
8+
## Plot
9+
using Winston
10+
module Plot3D
11+
using Winston
12+
# Warning: the 3D interface is very new at the time
13+
# of writing this and is likely to not be available
14+
# and/or have changed
15+
if Winston.output_surface == :gtk
16+
include(Pkg.dir("Winston","src","canvas3d_gtk.jl"))
17+
else
18+
include(Pkg.dir("Winston","src","canvas3d.jl"))
19+
end
20+
end
21+
xx = Float64[xx for xx = x, yy = x]
22+
yy = Float64[yy for xx = x, yy = x]
23+
Plot3D.plot3d(Plot3D.surf(xx,K,yy))
24+
25+
Plot3D.plot3d(Plot3D.surf(
26+
(u,v)->u,
27+
(u,v)->sin(u-v)./(u-v),
28+
(u,v)->v,
29+
x, x))

0 commit comments

Comments
 (0)