In [4]:
using LinearAlgebra
using Statistics
using SparseArrays
using MLDatasets
using PyPlot
Test of Johnson-Lindenstrauss transform¶
In [6]:
D = 100000;
n = 100;
Xs = randn(D,n); # the dataset
In [12]:
epsilon = 0.1;
d = Int64(ceil(1 * log(n) / epsilon^2))
Out[12]:
461
In [13]:
A = randn(d, D) / sqrt(d);
In [14]:
# let's loook at this for a single example
norm(Xs[:,1] - Xs[:,2])
Out[14]:
448.5205995177306
In [15]:
norm(A*Xs[:,1] - A*Xs[:,2])
Out[15]:
437.62667373218517
In [16]:
norm(A*Xs[:,1] - A*Xs[:,2]) / norm(Xs[:,1] - Xs[:,2])
Out[16]:
0.9757114259696008
In [17]:
AXs = A*Xs;
In [18]:
# over the whole dataset
extrema((i == j) ? 1.0 : norm(AXs[:,i] - AXs[:,j]) / norm(Xs[:,i] - Xs[:,j]) for i = 1:n, j = 1:n)
Out[18]:
(0.9001255539501946, 1.1102947985872953)
In [20]:
# how much did we decrease the dimension?
D/d
Out[20]:
216.91973969631238
Principal Component Analysis¶
In [22]:
n = 60000;
d = 28*28;
In [23]:
train_x, train_y = MNIST.traindata();
┌ Warning: MNIST.traindata() is deprecated, use `MNIST(split=:train)[:]` instead. └ @ MLDatasets ~/.julia/packages/MLDatasets/268Pz/src/datasets/vision/mnist.jl:187
This program has requested access to the data dependency MNIST.
which is not currently installed. It can be installed automatically, and you will not see this message again.
Dataset: THE MNIST DATABASE of handwritten digits
Authors: Yann LeCun, Corinna Cortes, Christopher J.C. Burges
Website: http://yann.lecun.com/exdb/mnist/
[LeCun et al., 1998a]
Y. LeCun, L. Bottou, Y. Bengio, and P. Haffner.
"Gradient-based learning applied to document recognition."
Proceedings of the IEEE, 86(11):2278-2324, November 1998
The files are available for download at the offical
website linked above. Note that using the data
responsibly and respecting copyright remains your
responsibility. The authors of MNIST aren't really
explicit about any terms of use, so please read the
website to make sure you want to download the
dataset.
Do you want to download the dataset from ["https://ossci-datasets.s3.amazonaws.com/mnist/train-images-idx3-ubyte.gz", "https://ossci-datasets.s3.amazonaws.com/mnist/train-labels-idx1-ubyte.gz", "https://ossci-datasets.s3.amazonaws.com/mnist/t10k-images-idx3-ubyte.gz", "https://ossci-datasets.s3.amazonaws.com/mnist/t10k-labels-idx1-ubyte.gz"] to "/Users/cmd353/.julia/datadeps/MNIST"?
[y/n]
In [24]:
train_x = Float64.(reshape(train_x, (d, n)));
In [25]:
train_x_mean = mean(train_x; dims=2);
In [26]:
train_x_minus_mean = train_x .- train_x_mean;
$$\Sigma = \frac{1}{n} \sum_{i=1}^n \left( x_i - \frac{1}{n} \sum_{j=1}^n x_j \right) \left( x_i - \frac{1}{n} \sum_{j=1}^n x_j \right)^T$$
In [27]:
Sigma = (train_x_minus_mean * train_x_minus_mean')/n;
In [28]:
# find the eigendecomposition of Sigma
ESigma = eigen(Sigma; sortby = (x -> -x))
Out[28]:
Eigen{Float64, Float64, Matrix{Float64}, Vector{Float64}}
values:
784-element Vector{Float64}:
5.116787728342094
3.7413284788648093
3.2526542396844857
2.8415733372665475
2.567074958812777
2.2736254916206096
1.7251262295281966
1.5205348976536206
1.4562809808876926
1.2427293764173302
1.1120709732035996
1.0666227659356202
0.9046657547344794
⋮
-2.622930252904868e-18
-2.7778853055429405e-18
-3.646087122795034e-18
-4.041423018984843e-18
-8.850995824136279e-18
-1.5896951408324886e-17
-1.7253002922741064e-17
-2.3746777919648282e-17
-3.262713073070834e-17
-4.8281967692188785e-17
-5.943455237030603e-17
-2.0801005196067098e-16
vectors:
784×784 Matrix{Float64}:
0.0 0.0 0.0 … 0.0 0.0
0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 … 0.0 0.0
0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 … 0.0 0.0
0.0 0.0 0.0 0.0 0.0
-1.13023e-6 1.79706e-6 7.96314e-7 0.000174968 7.72231e-6
⋮ ⋱
-8.769e-5 0.00143078 0.00113599 -7.5788e-14 -1.29722e-13
-4.5962e-5 0.000920176 0.000863292 6.40207e-14 3.66268e-13
-1.27477e-5 0.000500157 0.000487834 -1.98988e-14 -3.58485e-13
1.72528e-5 0.000219159 0.000211216 … -1.05846e-13 -2.26809e-14
5.51266e-6 0.000109011 9.13385e-5 1.16982e-13 2.45336e-13
-2.59364e-7 3.54844e-5 3.45235e-5 1.84993e-13 1.90297e-12
-4.95273e-7 2.12795e-5 1.89427e-5 -8.8161e-14 -1.49707e-11
1.37842e-7 2.22338e-6 2.25569e-6 0.0 5.97173e-11
0.0 0.0 0.0 … 0.0 0.0
0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0
In [29]:
# largest eigenvalue
ESigma.values[1]
Out[29]:
5.116787728342094
In [30]:
# corresponding eigenvector
ESigma.vectors[:, 1]
Out[30]:
784-element Vector{Float64}:
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
-1.1302299848127304e-6
⋮
-8.768996662291144e-5
-4.5961971626349e-5
-1.2747685974018059e-5
1.7252770690534547e-5
5.512661031492616e-6
-2.593635885644874e-7
-4.952733259032394e-7
1.378420261242352e-7
0.0
0.0
0.0
0.0
In [96]:
imshow(reshape(ESigma.vectors[:, 1], (28,28))');
In [97]:
imshow(reshape(ESigma.vectors[:, 2], (28,28))');
In [98]:
imshow(reshape(ESigma.vectors[:, 3], (28,28))');
In [31]:
imshow(reshape(ESigma.vectors[:, 4], (28,28))');
In [32]:
imshow(reshape(ESigma.vectors[:, 5], (28,28))');
In [115]:
semilogx(cumsum(ESigma.values))
ylabel("eigenvalue")
xlabel("index of eigenvalue")
title("Eigenvalues of MNIST Covariance Matrix");
PCA can represent objects in low dimension without losing information.
In [39]:
D = 100;
In [40]:
A = ESigma.vectors[:, 1:D]';
In [41]:
# original image
i = 1337;
imshow(reshape(train_x[:,i], (28,28))');
In [42]:
# dimension-reduced image
x_dr = A * train_x_minus_mean[:,i]
x_recovered = A' * x_dr + train_x_mean;
# original image
imshow(reshape(x_recovered, (28,28))');
Still enough to classify the image!
In [ ]:
Sparsity¶
In [43]:
n = 4096;
X = sprand(n,n,0.1);
Y = sprand(n,n,0.1);
In [47]:
X
Out[47]:
4096×4096 SparseMatrixCSC{Float64, Int64} with 1677115 stored entries:
⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿
⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿
⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿
⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿
⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿
⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿
⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿
⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿
⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿
⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿
⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿
⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿
⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿
⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿
⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿
⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿
⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿
⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿
⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿
⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿
In [44]:
# visualize part of it
imshow(X[1:64,1:64]);
colorbar();
In [45]:
# get the size of X in bytes
Base.summarysize(X)
Out[45]:
26866776
In [46]:
26866776/(n*n*8)
Out[46]:
0.20017307996749878
In [48]:
# type of X
typeof(X)
Out[48]:
SparseMatrixCSC{Float64, Int64}
In [49]:
Xdense = Matrix(X);
Ydense = Matrix(Y);
typeof(Xdense)
Out[49]:
Matrix{Float64} (alias for Array{Float64, 2})
In [50]:
Base.summarysize(Xdense)
Out[50]:
134217768
In [51]:
# what's the ratio?
Base.summarysize(X) / Base.summarysize(Xdense)
Out[51]:
0.20017302031128995
In [52]:
# what fraction of X's entries are nonzero?
nnz(X) / length(X)
Out[52]:
0.09996384382247925
Why are these numbers not equal?
In [88]:
@time X * Y;
0.827856 seconds (18 allocations: 281.658 MiB)
In [59]:
dense_mmpy_time = @elapsed Xdense * Ydense
Out[59]:
0.490777375
What happened here???
In [60]:
# make X and Y less dense
Xsparser = sprand(n,n,0.01);
Ysparser = sprand(n,n,0.01);
In [62]:
Xdense = Matrix(Xsparser);
Ydense = Matrix(Ysparser);
In [83]:
@time Xsparser * Ysparser';
0.090760 seconds (18 allocations: 97.028 MiB)
In [67]:
dense_mmpy_time = @elapsed Xdense * Ydense
Out[67]:
0.505052542
In [137]:
densities = [0.001, 0.003, 0.01, 0.03, 0.1, 0.3];
sparse_mm_times = Float64[];
for d in densities
Xsparse = sprand(n,n,d);
Ysparse = sprand(n,n,d);
push!(sparse_mm_times, @elapsed(Xsparse * Ysparse));
end
In [138]:
loglog(densities, sparse_mm_times; label="sparse mmpy");
loglog(densities, densities * 0 .+ dense_mmpy_time; label="dense mmpy");
xlabel("density");
ylabel("matrix multiply runtime");
legend();
In [66]:
spv = X[:,1]
Out[66]:
4096-element SparseVector{Float64, Int64} with 377 stored entries:
[9 ] = 0.476618
[14 ] = 0.731828
[24 ] = 0.993639
[25 ] = 0.29867
[26 ] = 0.6277
[49 ] = 0.99446
[91 ] = 0.446904
[121 ] = 0.463486
[126 ] = 0.532175
[128 ] = 0.857214
⋮
[3991] = 0.403587
[3992] = 0.484139
[4015] = 0.40613
[4017] = 0.406976
[4021] = 0.0818852
[4025] = 0.18207
[4042] = 0.184404
[4043] = 0.768553
[4048] = 0.450478
[4074] = 0.85291
[4076] = 0.68382
In [56]:
spv.n
Out[56]:
4096
In [57]:
spv.nzind
Out[57]:
429-element Vector{Int64}:
9
17
18
24
28
54
60
67
70
79
83
89
93
⋮
4010
4011
4032
4034
4036
4040
4057
4061
4062
4065
4066
4094
In [58]:
spv.nzval
Out[58]:
429-element Vector{Float64}:
0.3665466206937954
0.9898197098105949
0.8472553677537141
0.4200382543710388
0.575582096772398
0.478161953650292
0.5218372708763237
0.5337721068240063
0.6665316753775424
0.2454782188665694
0.6461326583138556
0.8885971852707479
0.3711005317479982
⋮
0.5222351482640311
0.15576228978743545
0.0989029063487934
0.0017576590602756959
0.5716640778660788
0.2710602459932675
0.8826343694387182
0.9838845254047324
0.8704305124740629
0.9049698008038184
0.6041680082714094
0.24830600332937425
In [ ]: