In [ ]:
using LinearAlgebra
using Statistics
using SparseArrays
using MLDatasets
using PyPlot

Test of Johnson-Lindenstrauss transform¶

In [2]:
d = 100000;
n = 100;

Xs = randn(d,n); # the dataset
In [3]:
epsilon = 0.1;
D = Int64(ceil(8 * log(n) / epsilon^2))
Out[3]:
3685
In [4]:
A = randn(D, d) / sqrt(D);
In [5]:
# let's loook at this for a single example
norm(Xs[:,1] - Xs[:,2])
Out[5]:
447.3594790493648
In [6]:
norm(A*Xs[:,1] - A*Xs[:,2])
Out[6]:
453.60422837944606
In [7]:
norm(A*Xs[:,1] - A*Xs[:,2]) / norm(Xs[:,1] - Xs[:,2])
Out[7]:
1.0139591304589126
In [8]:
AXs = A*Xs;
In [9]:
# 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[9]:
(0.9519405550668073, 1.046375370822268)
In [10]:
# how much did we decrease the dimension?
D/d
Out[10]:
0.03685

Principal Component Analysis¶

In [11]:
n = 60000;
d = 28*28;
In [12]:
train_x, train_y = MNIST.traindata();
In [13]:
train_x = Float64.(reshape(train_x, (d, n)));
In [14]:
train_x_mean = mean(train_x; dims=2);
In [15]:
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 [16]:
Sigma = (train_x_minus_mean * train_x_minus_mean')/n;
In [17]:
# find the eigendecomposition of Sigma
ESigma = eigen(Sigma; sortby = (x -> -x))
Out[17]:
Eigen{Float64,Float64,Array{Float64,2},Array{Float64,1}}
values:
784-element Array{Float64,1}:
  5.116787728342092
  3.7413284788648067
  3.2526542396844844
  2.8415733372665466
  2.5670749588127744
  2.273625491620609
  1.7251262295281955
  1.5205348976536206
  1.4562809808876922
  1.2427293764173304
  1.1120709732036005
  1.0666227659356209
  0.9046657547344794
  ⋮
 -4.341596436079342e-19
 -7.816133441131003e-19
 -8.328205317261318e-19
 -9.06846445474923e-19
 -2.7674683969166128e-18
 -7.155496110387914e-18
 -1.057007120922721e-17
 -2.366763579350205e-17
 -4.458784081324597e-17
 -6.16504989725162e-17
 -6.378346028260054e-17
 -9.992685906775396e-17
vectors:
784×784 Array{Float64,2}:
  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      -1.82916e-5   -1.69162e-5
  ⋮                                     ⋱                
 -8.769e-5    0.00143078   0.00113599       1.16734e-13   8.06924e-14
 -4.5962e-5   0.000920176  0.000863292      9.46519e-14  -2.08525e-13
 -1.27477e-5  0.000500157  0.000487834     -4.86915e-13   4.10592e-13
  1.72528e-5  0.000219159  0.000211216  …   1.07941e-12  -3.80997e-13
  5.51266e-6  0.000109011  9.13385e-5      -1.36204e-12   9.24027e-14
 -2.59364e-7  3.54844e-5   3.45235e-5       2.02711e-12  -1.66652e-13
 -4.95273e-7  2.12795e-5   1.89427e-5      -5.01332e-13   1.7726e-13
  1.37842e-7  2.22338e-6   2.25569e-6       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
In [18]:
# largest eigenvalue
ESigma.values[1]
Out[18]:
5.116787728342092
In [19]:
# corresponding eigenvector
ESigma.vectors[:, 1]
Out[19]:
784-element Array{Float64,1}:
  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.1302299848127289e-6
  ⋮
 -8.768996662290797e-5
 -4.596197162656801e-5
 -1.2747685973350732e-5
  1.7252770690812103e-5
  5.512661031686905e-6
 -2.593635886755097e-7
 -4.952733260038533e-7
  1.3784202612423542e-7
  0.0
  0.0
  0.0
  0.0
In [20]:
imshow(reshape(ESigma.vectors[:, 1], (28,28))');
In [21]:
imshow(reshape(ESigma.vectors[:, 2], (28,28))');
In [22]:
imshow(reshape(ESigma.vectors[:, 3], (28,28))');
In [23]:
imshow(reshape(ESigma.vectors[:, 4], (28,28))');
In [24]:
imshow(reshape(ESigma.vec
tors[:, 5], (28,28))');

syntax: missing comma or ) in argument list

Stacktrace:
 [1] top-level scope at In[24]:2
 [2] include_string(::Function, ::Module, ::String, ::String) at ./loading.jl:1091
In [25]:
plot(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 [26]:
D = 5;
In [27]:
A = ESigma.vectors[:, 1:D]';
In [28]:
# original image
i = 1337;
imshow(reshape(train_x[:,i], (28,28))');
In [29]:
# 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 [ ]:
In [30]:
Matrix(X)

UndefVarError: X not defined

Stacktrace:
 [1] top-level scope at In[30]:1
 [2] include_string(::Function, ::Module, ::String, ::String) at ./loading.jl:1091

Sparsity¶

In [31]:
n = 4096;
X = sprand(n,n,0.1);
Y = sprand(n,n,0.1);
In [32]:
# visualize part of it
imshow(X[1:64,1:64]);
colorbar();
In [33]:
# get the size of X in bytes
Base.summarysize(X)
Out[33]:
26896696
In [34]:
# type of X
typeof(X)
Out[34]:
SparseMatrixCSC{Float64,Int64}
In [35]:
Xdense = Matrix(X);
Ydense = Matrix(Y);
typeof(Xdense)
Out[35]:
Array{Float64,2}
In [36]:
Base.summarysize(Xdense)
Out[36]:
134217768
In [37]:
# what's the ratio?
Base.summarysize(X) / Base.summarysize(Xdense)
Out[37]:
0.20039594161631416
In [38]:
# what fraction of X's entries are nonzero?
nnz(X) / length(X)
Out[38]:
0.10007530450820923

Why are these numbers not equal?

In [39]:
@time X * Y;

  1.312289 seconds (110.34 k allocations: 287.229 MiB, 0.31% gc time)
In [40]:
dense_mmpy_time = @elapsed Xdense * Ydense
Out[40]:
1.174083337

What happened here???

In [41]:
# make X and Y less dense
Xsparser = sprand(n,n,0.01);
Ysparser = sprand(n,n,0.01);
@time Xsparser * Ysparser;

  0.121884 seconds (8 allocations: 94.578 MiB)
In [42]:
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 [43]:
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 [46]:
spv = X[:,1]
Out[46]:
4096-element SparseVector{Float64,Int64} with 427 stored entries:
  [16  ]  =  0.494497
  [18  ]  =  0.826087
  [34  ]  =  0.687849
  [59  ]  =  0.318373
  [63  ]  =  0.543537
  [73  ]  =  0.460325
  [75  ]  =  0.0852525
  [90  ]  =  0.141923
  [108 ]  =  0.0701148
  [116 ]  =  0.994332
          ⋮
  [4029]  =  0.570289
  [4035]  =  0.54624
  [4040]  =  0.202171
  [4042]  =  0.403688
  [4062]  =  0.420507
  [4066]  =  0.538902
  [4068]  =  0.265724
  [4072]  =  0.191174
  [4080]  =  0.67434
  [4084]  =  0.384076
  [4089]  =  0.188717
In [47]:
spv.n
Out[47]:
4096
In [48]:
spv.nzind
Out[48]:
427-element Array{Int64,1}:
   16
   18
   34
   59
   63
   73
   75
   90
  108
  116
  125
  127
  134
    ⋮
 4021
 4029
 4035
 4040
 4042
 4062
 4066
 4068
 4072
 4080
 4084
 4089
In [49]:
spv.nzval
Out[49]:
427-element Array{Float64,1}:
 0.4944967485040994
 0.8260867791524058
 0.6878492561841358
 0.31837265760532185
 0.5435373410214652
 0.4603252124070474
 0.08525252863099708
 0.14192308123381947
 0.07011482669433211
 0.994332364738097
 0.07054927165346192
 0.3335951909204913
 0.9303865975954675
 ⋮
 0.2442214611309792
 0.5702885838128535
 0.5462397434084996
 0.20217057750630207
 0.40368809673279826
 0.420507447700027
 0.5389015609303036
 0.2657240323430008
 0.19117367805893304
 0.674340219210197
 0.3840757164610593
 0.18871747643515846
In [ ]: