In [4]:

    

using PyPlot

push!(LOAD_PATH, "../src")

using HDStat

cca(U, V) = sum(svd(U' * V)[2]) / size(U, 2)


    

    
Out[4]:





cca (generic function with 1 method)

In [52]:

    

npt = 40
n, p, k = 1000, 800, 2

model = MPModel(n, p, 1)
    
sthresh = sv_infloor(model)
sv_out = zeros(npt)
sv_in = linspace(0, sthresh * 4, npt)
left, right = zeros(npt), zeros(npt)
for (ksv, sv) in enumerate(sv_in)
    U, V = qr(randn(n, k))[1], qr(randn(p, k))[1]
    X = rand(model) + sv * U * V'
    sv_out[ksv] = sqrt(eigs(X * X'; nev=k)[1][k])
    Up, _, Vp = svd(X)
    left[ksv], right[ksv] = cca(Up[:, 1:k], U), cca(Vp[:, 1:k], V)
end
    
sv_out_theory = map(x -> sv_xfer(model, x), sv_in)
tmp = map(x -> svec_overlap(model, x), sv_in)
left_theory, right_theory = map(x -> x[1], tmp), map(x -> x[2], tmp)


    

    
Out[52]:





(Real[0,0,0,0,0,0,0,0,0,0  …  0.9404266215010971,0.9442490859036637,0.9477129659773743,0.9508620643299124,0.9537336403548294,0.9563595593168388,0.958767210540193,0.9609802467426369,0.963019183599772,0.9649018891875982],Real[0,0,0,0,0,0,0,0,0,0  …  0.9505195451967994,0.9537950192792397,0.9567522268643246,0.9594314812598415,0.961866858951566,0.9640873306533126,0.9661176586036544,0.9679791140497407,0.9696900551733535,0.9712663958009776])

In [55]:

    

f = figure(figsize=(4, 3))
plot(sv_in, sv_out, ".")
plot(sv_in, sv_out_theory, "k-", linewidth=2)
axvline(sthresh, color="r", linestyle="--", linewidth=2)
axhline(sv_xfer(model, sthresh), color="r", linestyle="--", linewidth=2)
plot(sv_in, sv_in, "k-")
yticks([]); xticks([minimum(sv_in), int(maximum(sv_in))])
tight_layout()
savefig("../figures/sv_in_out.eps")

g = figure(figsize=(4, 3))
plot(sv_in, left, ".")
plot(sv_in, left_theory, "k-", linewidth=2)
axvline(sthresh, color="r", linestyle="--", linewidth=2)
yticks([]); xticks([minimum(sv_in), int(maximum(sv_in))])
tight_layout()
savefig("../figures/svec_in_out.eps")


    

In [7]:

    

f = figure(figsize=(4, 3))
N = 1000
k = 0.1
m = 0.3
K, M = int(k * N), int(m * N)
A = qr(randn(N, K))[1][randperm(N)[1:M], :]
es = eigs(A * A'; nev=K)[1]
# es = eig(A * A')[1]

lp, lm = (sqrt(k * (1 - m)) + sqrt(m * (1 - k)))^2, (sqrt(k * (1 - m)) - sqrt(m * (1 - k)))^2
ls = linspace(sqrt(lm), sqrt(lp), 101).^2
rho(l) = sqrt((l - lm) * (lp - l)) / 2 / pi / m / l / (1 - l)
rhos = map(rho, ls)
            
plt.hist(sqrt(abs(es)), normed=true, histtype="step", bins=20, linewidth=2, color="k")
plot(sqrt(ls), rhos * 2 .* sqrt(ls) * m / k, linewidth=2, color="r")

axvline(sqrt(lp), linestyle="--", linewidth=2, color="r")
axvline(sqrt(lm), linestyle="--", linewidth=2, color="r")
xticks([0, 1])
# ylim([0, maximum(sqrt(rhos)) * 1.2])
yticks([])
# title("Eigenvalue Spectrum of AA'")
            
tight_layout()
savefig("../figures/m.geq.k.eps")