In [1]:
using Color
using Gadfly
using Interact
using ProgressMeter

Wigner's semicircle law

Quite possibly the most famous result in random matrix theory.

Take a $N\times N$ real symmetric matrix with Gaussian entries and find its eigenvalues. The eigenvalues are distributed according to the semicircle law.

In [2]:
n=500
M=randn(n, n)
M=(M+M')/√2n

xg = linspace(-2, 2, 100)
plot(Coord.Cartesian,
    Guide.xlabel("Eigenvalue"), Guide.ylabel("Spectral density"),
    layer(x=xg, y=√abs(4 - xg.^2)/(2π), Geom.line,
        Theme(default_color=color("black"))
    ),
layer(x=eigvals(M), Geom.histogram(bincount=round(Int,√n), density=true))
)
Out[2]:
Eigenvalue -8 -7 -6 -5 -4 -3 -2 -1 0 1 2 3 4 5 6 7 8 9 -7.0 -6.8 -6.6 -6.4 -6.2 -6.0 -5.8 -5.6 -5.4 -5.2 -5.0 -4.8 -4.6 -4.4 -4.2 -4.0 -3.8 -3.6 -3.4 -3.2 -3.0 -2.8 -2.6 -2.4 -2.2 -2.0 -1.8 -1.6 -1.4 -1.2 -1.0 -0.8 -0.6 -0.4 -0.2 0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6 1.8 2.0 2.2 2.4 2.6 2.8 3.0 3.2 3.4 3.6 3.8 4.0 4.2 4.4 4.6 4.8 5.0 5.2 5.4 5.6 5.8 6.0 6.2 6.4 6.6 6.8 7.0 7.2 7.4 7.6 7.8 8.0 -10 -5 0 5 10 -7.0 -6.5 -6.0 -5.5 -5.0 -4.5 -4.0 -3.5 -3.0 -2.5 -2.0 -1.5 -1.0 -0.5 0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 4.5 5.0 5.5 6.0 6.5 7.0 7.5 8.0 -0.5 -0.4 -0.3 -0.2 -0.1 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 -0.40 -0.38 -0.36 -0.34 -0.32 -0.30 -0.28 -0.26 -0.24 -0.22 -0.20 -0.18 -0.16 -0.14 -0.12 -0.10 -0.08 -0.06 -0.04 -0.02 0.00 0.02 0.04 0.06 0.08 0.10 0.12 0.14 0.16 0.18 0.20 0.22 0.24 0.26 0.28 0.30 0.32 0.34 0.36 0.38 0.40 0.42 0.44 0.46 0.48 0.50 0.52 0.54 0.56 0.58 0.60 0.62 0.64 0.66 0.68 0.70 0.72 0.74 0.76 0.78 0.80 0.82 -0.5 0.0 0.5 1.0 -0.4 -0.3 -0.2 -0.1 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 Spectral density

Wigner's original paper (1955)...

... was actually not on Gaussian random matrices, but rather, other matrices motivated from nuclear physics.

<img width=700 src="wigner55head.png" />

This paper considers 3 random matrix problems surrounding this type of matrix:

In [3]:
function wigner55(n::Int, N::Int=n, v::Real=1.0)
    T = promote_type(typeof(v*1), typeof(n))
    H = zeros(T, 2n+1, 2n+1)
    for i=1:2n+1
        H[i,i] = i-1-n
    end
    for m = 1:N, j=1:2n+1-m
        H[j,j+m] = v*(2randbool()-1)
    end
    Symmetric(H, :U)
end
wigner55(2, 1, 1)
Out[3]:
5x5 Symmetric{Int64}:
 -2  -1   0   0   0
 -1  -1   1   0   0
  0   1   0  -1   0
  0   0  -1   1  -1
  0   0   0  -1   2

Wigner studied three problems:

  1. The bandwidth $N = 1$, dimension $n \rightarrow \infty$, with nonzero diagonals.
  2. $n = N$ and $v$ are all finite, with zero diagonals.
  3. $N \rightarrow \infty, n \rightarrow \infty, v \rightarrow \infty, v^2/N = $ some constant $q$, with nonzero diagonals.

Problem 1. $N=1$, $n\rightarrow \infty$, nonzero diagonals

Eigenvectors (Eq. 5b)

...

In [4]:
function demo1(n::Int=20)
    @manipulate for v=0.0:0.25:5.0, k=1:2n+1
        E = eigfact!(wigner55(n, 1, v))
        plot(Coord.Cartesian(ymin=-1, ymax=1),
            Guide.ylabel("Eigenvector"), Guide.xlabel("Position"),
            layer(x=-n:n, y=E[:vectors][:,k],
                Geom.point, Theme(default_color=color("black"))),
            layer(x=-n:n, y=map(m->abs(besselj(m-k+(n+1), 2v)), -n:n),
                Geom.line, Theme(default_color=color("lightgrey"))),
            layer(x=-n:n, y=map(m->-abs(besselj(m-k+(n+1), 2v)), -n:n),
                Geom.line, Theme(default_color=color("lightgrey"))),
        )
    end
end
demo1()
Out[4]:
Position -70 -60 -50 -40 -30 -20 -10 0 10 20 30 40 50 60 70 -60 -58 -56 -54 -52 -50 -48 -46 -44 -42 -40 -38 -36 -34 -32 -30 -28 -26 -24 -22 -20 -18 -16 -14 -12 -10 -8 -6 -4 -2 0 2 4 6 8 10 12 14 16 18 20 22 24 26 28 30 32 34 36 38 40 42 44 46 48 50 52 54 56 58 60 -60 -30 0 30 60 -60 -55 -50 -45 -40 -35 -30 -25 -20 -15 -10 -5 0 5 10 15 20 25 30 35 40 45 50 55 60 -3.5 -3.0 -2.5 -2.0 -1.5 -1.0 -0.5 0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 -3.0 -2.9 -2.8 -2.7 -2.6 -2.5 -2.4 -2.3 -2.2 -2.1 -2.0 -1.9 -1.8 -1.7 -1.6 -1.5 -1.4 -1.3 -1.2 -1.1 -1.0 -0.9 -0.8 -0.7 -0.6 -0.5 -0.4 -0.3 -0.2 -0.1 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8 1.9 2.0 2.1 2.2 2.3 2.4 2.5 2.6 2.7 2.8 2.9 3.0 -4 -2 0 2 4 -3.0 -2.8 -2.6 -2.4 -2.2 -2.0 -1.8 -1.6 -1.4 -1.2 -1.0 -0.8 -0.6 -0.4 -0.2 0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6 1.8 2.0 2.2 2.4 2.6 2.8 3.0 Eigenvector

Problem 2. $N=2n$, $n, v$ finite, 0 diagonals

In [5]:
function wigner55_2(n::Int, N::Int=n, v::Real=1.0)
    T = promote_type(typeof(v*1), typeof(n))
    H = zeros(T, 2n+1, 2n+1) #Diagonals are now 0
    for m = 1:N, j=1:2n+1-m
        H[j,j+m] = v*(2randbool()-1)
    end
    Symmetric(H, :U)
end
wigner55_2(3, 6, 1)
Out[5]:
7x7 Symmetric{Int64}:
  0   1   1   1   1  -1  -1
  1   0  -1   1  -1   1   1
  1  -1   0  -1  -1   1  -1
  1   1  -1   0   1  -1  -1
  1  -1  -1   1   0   1  -1
 -1   1   1  -1   1   0  -1
 -1   1  -1  -1  -1  -1   0
In [6]:
t = 1000 #Number of trials
n = 30 
v = 4 #Magnitude of off-diagonal
evals = zeros(2n+1, t)

@showprogress for i=1:t
    H1 = wigner55_2(n, 2n, v)
    evals[:, i] = eigvals(H1)
end
Progress: 100%|█████████████████████████████████████████| Time: 0:00:04
In [7]:
xg = linspace(-v*√8n, v*√8n, 100)
plot(Coord.Cartesian(aspect_ratio=φ),
    Guide.xlabel("Eigenvalue"), Guide.ylabel("Spectral density"),
    layer(x=xg, y=√abs(8n*v^2 - xg.^2)/(4π*n*v^2), Geom.line,
        Theme(default_color=color("black"))
    ),
    layer(x=evals, Geom.histogram(bincount=round(Int,(t*(2n+1))), density=true))
)
Out[7]:
Eigenvalue -350 -300 -250 -200 -150 -100 -50 0 50 100 150 200 250 300 350 -300 -290 -280 -270 -260 -250 -240 -230 -220 -210 -200 -190 -180 -170 -160 -150 -140 -130 -120 -110 -100 -90 -80 -70 -60 -50 -40 -30 -20 -10 0 10 20 30 40 50 60 70 80 90 100 110 120 130 140 150 160 170 180 190 200 210 220 230 240 250 260 270 280 290 300 -400 -200 0 200 400 -300 -280 -260 -240 -220 -200 -180 -160 -140 -120 -100 -80 -60 -40 -20 0 20 40 60 80 100 120 140 160 180 200 220 240 260 280 300 -0.020 -0.015 -0.010 -0.005 0.000 0.005 0.010 0.015 0.020 0.025 0.030 0.035 -0.0150 -0.0145 -0.0140 -0.0135 -0.0130 -0.0125 -0.0120 -0.0115 -0.0110 -0.0105 -0.0100 -0.0095 -0.0090 -0.0085 -0.0080 -0.0075 -0.0070 -0.0065 -0.0060 -0.0055 -0.0050 -0.0045 -0.0040 -0.0035 -0.0030 -0.0025 -0.0020 -0.0015 -0.0010 -0.0005 0.0000 0.0005 0.0010 0.0015 0.0020 0.0025 0.0030 0.0035 0.0040 0.0045 0.0050 0.0055 0.0060 0.0065 0.0070 0.0075 0.0080 0.0085 0.0090 0.0095 0.0100 0.0105 0.0110 0.0115 0.0120 0.0125 0.0130 0.0135 0.0140 0.0145 0.0150 0.0155 0.0160 0.0165 0.0170 0.0175 0.0180 0.0185 0.0190 0.0195 0.0200 0.0205 0.0210 0.0215 0.0220 0.0225 0.0230 0.0235 0.0240 0.0245 0.0250 0.0255 0.0260 0.0265 0.0270 0.0275 0.0280 0.0285 0.0290 0.0295 0.0300 -0.02 0.00 0.02 0.04 -0.015 -0.014 -0.013 -0.012 -0.011 -0.010 -0.009 -0.008 -0.007 -0.006 -0.005 -0.004 -0.003 -0.002 -0.001 0.000 0.001 0.002 0.003 0.004 0.005 0.006 0.007 0.008 0.009 0.010 0.011 0.012 0.013 0.014 0.015 0.016 0.017 0.018 0.019 0.020 0.021 0.022 0.023 0.024 0.025 0.026 0.027 0.028 0.029 0.030 Spectral density

Problem 3. $N \rightarrow \infty$, $v = \sqrt{qN}$, nonzero diagonals

In [8]:
q = 100
n = 50
N = 2n
t = 10 #Number of trials
evals = zeros(2n+1, t)

@showprogress for i=1:t
    H = wigner55(n, N, (q*N))
    evals[:, i] = eigvals(H)
end
Progress: 100%|█████████████████████████████████████████| Time: 0:00:00
In [9]:
v=(q*N)
xg = linspace(-v*√4N, v*√4N, 100)
plot(Coord.Cartesian(aspect_ratio=φ),
    Guide.xlabel("Eigenvalue"), Guide.ylabel("Spectral density"),
    layer(x=xg, y=√abs(4N*v^2 - xg.^2)/(2π*N*v^2), Geom.line,
        Theme(default_color=color("black"))
    ),
    layer(x=evals, Geom.histogram(bincount=round(Int,(t*(2n+1))), density=true))
)
Out[9]:
Eigenvalue -7×10³ -6×10³ -5×10³ -4×10³ -3×10³ -2×10³ -1×10³ 0 1×10³ 2×10³ 3×10³ 4×10³ 5×10³ 6×10³ 7×10³ -6.0×10³ -5.8×10³ -5.6×10³ -5.4×10³ -5.2×10³ -5.0×10³ -4.8×10³ -4.6×10³ -4.4×10³ -4.2×10³ -4.0×10³ -3.8×10³ -3.6×10³ -3.4×10³ -3.2×10³ -3.0×10³ -2.8×10³ -2.6×10³ -2.4×10³ -2.2×10³ -2.0×10³ -1.8×10³ -1.6×10³ -1.4×10³ -1.2×10³ -1.0×10³ -8.0×10² -6.0×10² -4.0×10² -2.0×10² 0 2.0×10² 4.0×10² 6.0×10² 8.0×10² 1.0×10³ 1.2×10³ 1.4×10³ 1.6×10³ 1.8×10³ 2.0×10³ 2.2×10³ 2.4×10³ 2.6×10³ 2.8×10³ 3.0×10³ 3.2×10³ 3.4×10³ 3.6×10³ 3.8×10³ 4.0×10³ 4.2×10³ 4.4×10³ 4.6×10³ 4.8×10³ 5.0×10³ 5.2×10³ 5.4×10³ 5.6×10³ 5.8×10³ 6.0×10³ -6×10³ -3×10³ 0 3×10³ 6×10³ -6.0×10³ -5.5×10³ -5.0×10³ -4.5×10³ -4.0×10³ -3.5×10³ -3.0×10³ -2.5×10³ -2.0×10³ -1.5×10³ -1.0×10³ -5.0×10² 0 5.0×10² 1.0×10³ 1.5×10³ 2.0×10³ 2.5×10³ 3.0×10³ 3.5×10³ 4.0×10³ 4.5×10³ 5.0×10³ 5.5×10³ 6.0×10³ -0.0005 -0.0004 -0.0003 -0.0002 -0.0001 0.0000 0.0001 0.0002 0.0003 0.0004 0.0005 0.0006 0.0007 0.0008 0.0009 -0.00040 -0.00038 -0.00036 -0.00034 -0.00032 -0.00030 -0.00028 -0.00026 -0.00024 -0.00022 -0.00020 -0.00018 -0.00016 -0.00014 -0.00012 -0.00010 -0.00008 -0.00006 -0.00004 -0.00002 0.00000 0.00002 0.00004 0.00006 0.00008 0.00010 0.00012 0.00014 0.00016 0.00018 0.00020 0.00022 0.00024 0.00026 0.00028 0.00030 0.00032 0.00034 0.00036 0.00038 0.00040 0.00042 0.00044 0.00046 0.00048 0.00050 0.00052 0.00054 0.00056 0.00058 0.00060 0.00062 0.00064 0.00066 0.00068 0.00070 0.00072 0.00074 0.00076 0.00078 0.00080 -0.0005 0.0000 0.0005 0.0010 -0.00040 -0.00035 -0.00030 -0.00025 -0.00020 -0.00015 -0.00010 -0.00005 0.00000 0.00005 0.00010 0.00015 0.00020 0.00025 0.00030 0.00035 0.00040 0.00045 0.00050 0.00055 0.00060 0.00065 0.00070 0.00075 0.00080 Spectral density
In [10]:
@manipulate for logq = -5:0.1:5
    
    q = 10.0^logq
    n = 50
    N = 2n
    t = 10 #Number of trials
    evals = zeros(2n+1, t)

    for i=1:t
        H = wigner55(n, N, (q*N))
        evals[:, i] = eigvals(H)
    end    

    v=(q*N)
    xg = linspace(-v*√4N, v*√4N, 100)
    plot(Coord.Cartesian(aspect_ratio=φ),
        Guide.xlabel("Eigenvalue"), Guide.ylabel("Spectral density"),
        layer(x=xg, y=√abs(4N*v^2 - xg.^2)/(2π*N*v^2), Geom.line,
    Theme(default_color=color("grey"))
        ),
    layer(x=[-n; -n:n; n], y=[0; fill(1/(2n+1), 2n+1); 0], Geom.line,
    Theme(default_color=color("grey"))
        ),
        layer(x=evals, Geom.histogram(bincount=round(Int,(t*(2n+1))), density=true))
    )
end
Out[10]:
Eigenvalue -1000 -900 -800 -700 -600 -500 -400 -300 -200 -100 0 100 200 300 400 500 600 700 800 900 1000 -900 -880 -860 -840 -820 -800 -780 -760 -740 -720 -700 -680 -660 -640 -620 -600 -580 -560 -540 -520 -500 -480 -460 -440 -420 -400 -380 -360 -340 -320 -300 -280 -260 -240 -220 -200 -180 -160 -140 -120 -100 -80 -60 -40 -20 0 20 40 60 80 100 120 140 160 180 200 220 240 260 280 300 320 340 360 380 400 420 440 460 480 500 520 540 560 580 600 620 640 660 680 700 720 740 760 780 800 820 840 860 880 900 -1000 -500 0 500 1000 -900 -850 -800 -750 -700 -650 -600 -550 -500 -450 -400 -350 -300 -250 -200 -150 -100 -50 0 50 100 150 200 250 300 350 400 450 500 550 600 650 700 750 800 850 900 -0.015 -0.010 -0.005 0.000 0.005 0.010 0.015 0.020 0.025 -0.0100 -0.0095 -0.0090 -0.0085 -0.0080 -0.0075 -0.0070 -0.0065 -0.0060 -0.0055 -0.0050 -0.0045 -0.0040 -0.0035 -0.0030 -0.0025 -0.0020 -0.0015 -0.0010 -0.0005 0.0000 0.0005 0.0010 0.0015 0.0020 0.0025 0.0030 0.0035 0.0040 0.0045 0.0050 0.0055 0.0060 0.0065 0.0070 0.0075 0.0080 0.0085 0.0090 0.0095 0.0100 0.0105 0.0110 0.0115 0.0120 0.0125 0.0130 0.0135 0.0140 0.0145 0.0150 0.0155 0.0160 0.0165 0.0170 0.0175 0.0180 0.0185 0.0190 0.0195 0.0200 -0.02 -0.01 0.00 0.01 0.02 -0.011 -0.010 -0.009 -0.008 -0.007 -0.006 -0.005 -0.004 -0.003 -0.002 -0.001 0.000 0.001 0.002 0.003 0.004 0.005 0.006 0.007 0.008 0.009 0.010 0.011 0.012 0.013 0.014 0.015 0.016 0.017 0.018 0.019 0.020 Spectral density

Some free probability experiments

In [11]:
function samplefree(n, N, t, q)
    v=(q*N)

    evals1 = zeros(2n+1, t)
    evals2 = zeros(2n+1, t)
    
    #Sample eigenvalues
    @showprogress for i=1:t
        𝐤 = Diagonal(1.0*[-n:n;])
        𝐯 = wigner55_2(n,2n,v) |> full
        Q = full(qrfact!(randn(2n+1, 2n+1))[:Q])
        evals1 = eigvals(𝐤 + 𝐯)
        evals2 = eigvals(Q*full(𝐤)*Q' + 𝐯)
    end
    
    #Compute normalized histograms
    xv, yv1 = hist(evals1, round(Int, (n*√t)))
    yv1 /= sum(yv1)*diff(xv)[1]
    xv, yv2 = hist(evals2, xv)
    yv2 /= sum(yv2)*diff(xv)[1]

    midpoints(xv), yv1, yv2
end
Out[11]:
samplefree (generic function with 1 method)
In [12]:
function plotdensities(xv, yv1, yv2)
    cm = map(x->AlphaColorValue(x, 0.5), palette("Set1", 3)[1:2])
    plot(Guide.xlabel("Eigenvalue"), Guide.ylabel("Density"),
        layer(x=xv, y=yv1, Geom.bar, Theme(default_color=cm[1])),
        layer(x=xv, y=yv2, Geom.bar, Theme(default_color=cm[2])),
        Guide.manual_color_key("", ["exact", "free"], cm), 
    )
end
Out[12]:
plotdensities (generic function with 1 method)
In [13]:
@manipulate for logq = -5:0.5:5
    xv, yv1, yv2 = samplefree(30, 2n, 100, 10.0^logq)
    plotdensities(xv, yv1, yv2)
end
Progress: 100%|█████████████████████████████████████████| Time: 0:00:01
palette not defined
while loading In[13], in expression starting on line 1

 in plotdensities at In[12]:2
 in anonymous at no file:3
 in lift at /Users/jiahao/.julia/v0.3/Reactive/src/Reactive.jl:298