In [1]:
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]:
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]:
In [6]:
norm(A*Xs[:,1] - A*Xs[:,2])
Out[6]:
In [7]:
norm(A*Xs[:,1] - A*Xs[:,2]) / norm(Xs[:,1] - Xs[:,2])
Out[7]:
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]:
In [10]:
# how much did we decrease the dimension?
D/d
Out[10]:
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]:
In [18]:
# largest eigenvalue
ESigma.values[1]
Out[18]:
In [19]:
# corresponding eigenvector
ESigma.vectors[:, 1]
Out[19]:
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))');
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)
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]:
In [34]:
# type of X
typeof(X)
Out[34]:
In [35]:
Xdense = Matrix(X);
Ydense = Matrix(Y);
typeof(Xdense)
Out[35]:
In [36]:
Base.summarysize(Xdense)
Out[36]:
In [37]:
# what's the ratio?
Base.summarysize(X) / Base.summarysize(Xdense)
Out[37]:
In [38]:
# what fraction of X's entries are nonzero?
nnz(X) / length(X)
Out[38]:
Why are these numbers not equal?
In [39]:
@time X * Y;
In [40]:
dense_mmpy_time = @elapsed Xdense * Ydense
Out[40]:
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;
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]:
In [47]:
spv.n
Out[47]:
In [48]:
spv.nzind
Out[48]:
In [49]:
spv.nzval
Out[49]:
In [ ]: