Eigenvalues of Triangle Normalized Laplacian in Julia Code

triangle-law.jl

##
using LinearAlgebra
using SparseArrays
using StatsBase
using Plots
using MatrixNetworks 
# for Julia 1.0, you need to run 
# using Pkg; Pkg.add("MatrixNetworks#2f5c59c")
# to get the latest 1.0 compatible version. 
pyplot()
##
ergraph(n,d) = Float64.(sparse(Symmetric(triu!(sprand(Bool, n, n, d/n),1))))
function nlaplacian(A)
  id = 1.0./sqrt.(vec(sum(A,dims=2)))
  I,J,V = findnz(A)
  return sparse(I,J,V.*(id[I].*id[J]),size(A)...)
end
##
avgd = [10,20,30]
n = 200
nt = 150
lams = Vector{Vector{Float64}}()
for d in avgd
  lamsd = zeros(0)
  for t=1:nt
    A = ergraph(n,d)
    T = (A*A).*A
    L = nlaplacian(largest_component(T)[1])
    push!(lamsd, eigvals!(Matrix(L))[1:end-1]...)
  end
  push!(lams, lamsd)
end
##
plot(size=(400,250))
for (i,d) in enumerate(avgd)
  lamsd = lams[i]
  h = fit(Histogram, (abs.(1.0.-lamsd)); nbins = 50)
  h = normalize(h)
  plot!(h.edges, h.weights, label="$d")
end
mprho = 0.038
mpa = (1-sqrt(mprho))^2
mpb = (1+sqrt(mprho))^2
@show 2.0-mpa, 2.0-mpb
xs = collect(range(0,stop=2,length=100))
plot!(xs, map(lam -> begin
  lam = 2.0-lam
  if mpa <= lam <= mpb
    return sqrt.((mpb - lam)*(lam-mpa))./(2π*lam*mprho)
  else
    return 0.0
  end
end, xs))
savefig("triangle-law-$n-1.pdf")
gui()


## Now look at the regular eigenvalue laws
avgd = [10,20,30]
n = 200
nt = 150
lams = Vector{Vector{Float64}}()
for d in avgd
  lamsd = zeros(0)
  for t=1:nt
    A = ergraph(n,d)
    T = A
    L = nlaplacian(largest_component(T)[1])
    push!(lamsd, eigvals!(Matrix(L))[1:end-1]...)
  end
  push!(lams, lamsd)
end
##
plot(size=(400,250))
for (i,d) in enumerate(avgd)
  lamsd = lams[i]
  h = fit(Histogram, (abs.(1.0.-lamsd)); nbins = 50)
  h = normalize(h)
  plot!(h.edges, h.weights, label="$d")
end
C = 0.19
xs = collect(range(0,stop=2,length=100))
mpa = 1-sqrt(C)
mpb = 1+sqrt(C)
plot!(xs, map(lam -> begin
  if mpa <= lam <= mpb
    lam = (1.0-lam)/sqrt(C)
    return 2*sqrt.((1-lam.^2))/(sqrt(C)*π)
  else
    return 0.0
  end
end, xs))
gui()
savefig("eigenvalue-law-$n-1.pdf")

##
plot(size=(400,250))
for (i,d) in enumerate(avgd)
  lamsd = lams[i]
  h = fit(Histogram, (abs.(1.0.-lamsd)); nbins = 50)
  h = normalize(h)
  plot!(h.edges, h.weights, label="$d")
end
C = 0.11
xs = collect(range(0,stop=2,length=100))
mpa = 1-sqrt(C)
mpb = 1+sqrt(C)
plot!(xs, map(lam -> begin
  if mpa <= lam <= mpb
    lam = (1.0-lam)/sqrt(C)
    return 2*sqrt.((1-lam.^2))/(sqrt(C)*π)
  else
    return 0.0
  end
end, xs))
gui()
savefig("eigenvalue-law-$n-2.pdf")

zonotope-sampling.jl

##
using Plots
pyplot()
##
function num_vertices(m::Int,n::Int)
  N::Int = 0
  for i=0:n-1
    N = N + binomial(m-1,i)
  end
  N = 2*N
end
function zonotope_vertices_rand_fast(W::Array{Float64,2};nmax=10^8)
  m,n = size(W)
  rand_size = 100

  nverts = num_vertices(m,n)
  ntrials = 0

  Vcur=Array{Float64,2}(undef,0,m)
  while size(Vcur,1) < nverts && ntrials < nmax
    Z = randn(rand_size,n)
    Y = Z*W'
    Y = sign.(Y)
    V = unique(Y,dims=1)
    Vcur = unique([Vcur;V;-V],dims=1)
    ntrials += rand_size
  end
  if ntrials >= nmax && size(Vcur,1) < nverts
    warn(string("The zonotope sampling algorithm did not complete, ",
                "missing $(nverts-size(Vcur,1)) vertices"))
  end
  return Vcur*W
end
##
using Random
#Random.seed!(5)
Random.seed!(10)
G = randn(2,10)
#G = [1.5 1 2.5; -0.5 3 2]
V = zonotope_vertices_rand_fast(copy(G'))'
p = scatter(V[1,:],V[2,:],marker=:circle,
  markersize=12,label="",
  markercolor=nothing,
  size=(300,300))
gui()
##
p = scatter(V[1,:],V[2,:],marker=:circle,
  markersize=12,label="",
  markercolor=nothing,
  size=(300,300))
plot!(framestyle=:grid)
p2 = scatter!([V[1,1]],[V[2,1]],label="",markercolor=:red,markerstrokecolor=nothing,alpha=0.1)
xlims!(-12,12)
ylims!(-10,10)
wd = 0.2/25
anim = @animate for i=1:5000
  #x = G'*sign.(G*randn(3))
  x = G*sign.(G'*randn(2))
  #push!(p[2], [V[1,j]+wd*randn()],[V[2,j]+wd*randn()])
  push!(p2[1][2], x[1]+sqrt(i)*wd*randn(),x[2]+sqrt(i)*wd*randn())
  title!("$i samples")
end every 50
gui()
gif(anim, "sampling-complex.gif")