chapter 7 (7.1 and 7.2) eigenvalue sums

Author: vtjnash

demos/book/6/.gitignore

@@ -0,0 +1,38 @@
+#centrallimit.jl
+#Algorithm 7.2 of Random Eigenvalues by Alan Edelman
+
+using Winston
+
+function centrallimit(m, t)
+    # sums m random matrices and histograms their eigenvalues
+    # normalized by 1/sqrt(m)
+    # Input :
+    #  m : number of matrices to sum
+    #  t : number of trials to perform
+    n = 100     # matrix size
+    dx = .1     # bin size
+    v = Float64[]
+    # we choose a random diagonal {-1,1} matrx as our prototype
+    # A = diagm(sign(randn(1,n)))
+    # proto = kron([1 0; 0 -1], eye(n/2))
+    for i = 1:t                         # for each trial
+        B = zeros(n,n)                  # initialize the accumulator
+        for z = 1:m                     # sum up the random matrices
+            Q = full(QR(randn(n,n) + im*randn(n,n))[:Q])  # (piecewise) Haar orthogonal
+            A = diagm(sign(randn(n,1)))
+            B = B + Q'*A*Q              # A and Q'AQ are asymptotically free
+        end
+        append!(v, real(eigvals(B)/sqrt(m)))  # normalize eigenvalues
+    end
+    grid, count = hist(v, -2-dx/2:dx:2)
+
+    p = FramedPlot()
+    h = Histogram(count*pi/2/(dx*n*t), step(grid))
+    h.x0 = first(grid)
+    add(p, h)
+    if isinteractive()
+        Winston.display(p)
+    else
+        file(p, "centrallimit_m=$(m)_t=$(t).png")
+    end
+end

demos/book/7/centrallimit.jl

@@ -0,0 +1,49 @@
+#eigfreesum.jl
+#Algorithm 7.1 of Random Eigenvalues by Alan Edelman
+
+#Experiment:    n by n diagonal matrices with random entries +1 or -1
+#Plot:          classical and free sums
+#Theory:        arcsine distribution
+
+## Parameters
+t = 1           # number of trials
+n = 1000        # size of matrix
+dx = .1
+
+## Experiment
+function eigfreesum_experiment(t, n, dx)
+    v = Float64[]
+    w = Float64[]
+    for i = 1:t
+        a = sign(randn(n))
+        b = sign(randn(n))
+        Q = full(QR(randn(n,n))[:Q])                    # random orthogonal matrix
+        append!(v, a+b)                                 # classical sum
+        append!(w, eigvals(diagm(a) + Q'*diagm(b)*Q))   # free sum
+    end
+    x = -2-dx/2:dx:2+dx/2
+    _,countw = hist(w, x)
+    _,countv = hist(v, x)
+    x, countw, countv
+end
+grid, countw, countv = eigfreesum_experiment(t, n, dx)
+
+## Theory
+x = -2+dx/20:dx/10:2-dx/20
+y = (1/pi) ./ sqrt(4 - x.^2)
+
+## Plot
+using Winston
+p = FramedPlot()
+hw = Histogram(countw*pi/(n*t*dx), step(grid), "color", "blue")
+hw.x0 = first(grid)
+add(p, hw)
+hv = Histogram(countv*pi/(n*t), step(grid), "color", "green")
+hv.x0 = first(grid)
+add(p, hv)
+add(p, Curve(x, y, "color", "red", "linewidth", 2))
+if isinteractive()
+    Winston.display(p)
+else
+    file(p, "eigfreesum.png")
+end