chapter 6 of rmt (5.6 to 5.8)

Author: vtjnash

demos/book/5/finitesemi.jl

@@ -0,0 +1,44 @@
+#finitesemi.jl
+#Algorithm 5.6 of Random Eigenvalues by Alan Edelman
+
+#Experiment:    Eigenvalues of GUE matrices
+#Plot:          Histogram of eigenvalues
+#Theory:        Semicircle and finite semicircle
+
+## Parameters
+n = 3           # size of matrix
+s = 10000       # number of samples
+d = 0.1         # bin size
+
+## Experiment
+function finitesemi_experiment(n,s,d)
+    e = Float64[]
+    for i = 1:s
+        a = randn(n,n)+im*randn(n,n)
+        a = (a+a')/(2*sqrt(4n))
+        v = eigvals(a)
+        append!(e, v)
+    end
+    return hist(e, -1.5:d:1.5)
+end
+grid, count = finitesemi_experiment(n,s,d)
+
+## Theory
+x = -1:0.01:1
+y = sqrt(1-x.^2)
+
+## Plot
+using Winston
+p = FramedPlot()
+h = Histogram(count*pi/(2*d*n*s), step(grid))
+h.x0 = first(grid)
+add(p, h)
+add(p, Curve(x, y, "color", "red", "linewidth", 2))
+include("levels.jl")
+add(p, Curve(levels(n)..., "color", "blue", "linewidth", 2))
+if isinteractive()
+    Winston.display(p)
+else
+    file(p, "spherecomponent.png")
+end
+

demos/book/5/levels.jl

@@ -0,0 +1,27 @@
+#levels.jl
+#Algorithm 5.7 of Random Eigenvalues by Alan Edelman
+
+function levels(n)
+    # Plot exact semicircle formula for GUE
+    # This is given in formula (5.2.16) in Mehta as the sum of phi_j^2
+    # but Christoffel-Darboux seems smarter
+    ## Parameter
+    # n: Dimension of matrix for which to determine exact semicircle
+    x = (-1:0.001:1) * sqrt(2n) * 1.3
+    pold = zeros(size(x))   # -1st Hermite Polynomial
+    p = ones(size(x))       #  0th Hermite Polynomial
+    k = p
+    for j = 1:n # Compute the three-term recurrence
+        pnew = (sqrt(2) * x .* p - sqrt(j - 1) * pold) / sqrt(j)
+        pold = p
+        p = pnew
+    end
+    pnew = (sqrt(2) * x .* p - sqrt(n) * pold)/sqrt(n + 1)
+
+    k = n * p .^ 2 - sqrt(n * (n+1)) * pnew .* pold # Use p.420 of Mehta
+    k = k .* exp(-x .^ 2) / sqrt(pi) # Normalize
+    
+    # Rescale so that "semicircle" is on [-1, 1] and area is pi/2
+    (x / sqrt(2n), k * pi / sqrt(2n))
+end
+

demos/book/5/levels2.jl

@@ -0,0 +1,30 @@
+#levels2.jl
+#Algorithm 5.8 of Random Eigenvalues by Alan Edelman
+
+function levels2(n)
+    # Plot exact semicircle formula for GUE
+    xfull = (-1:0.001:1) * sqrt(2n) * 1.3
+
+    # Form the T matrix
+    T = diagm(sqrt(1:n-1),1)
+    T = T + T'
+    
+    # Do the eigen-decomposition of T, T = UVU'
+    V, U = eig(T)
+
+    # Precompute U'*e_n
+    # tmp_en = U' * ((0:n-1) == n-1)'
+    tmp_en = U[end, :]'
+
+    y = Array(Float64, length(xfull))
+    for i = 1:length(xfull)
+        x = xfull[i]
+        # generate the v vector as in (2.5)
+        v = U * (tmp_en ./ (V - sqrt(2) * x))
+        # multiply the normalization term
+        y[i] = norm((sqrt(pi)) ^ (-1/2) * exp(-x^2/2) * v/v[1])^2
+    end
+
+    # Rescale so that "semicircle" is on [-1, 1] and area is pi/2
+    (xfull/sqrt(2n), y * pi / sqrt(2n))
+end