In [42]:
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 [33]:
n = 60000;
d = 28*28;
In [34]:
train_x, train_y = MNIST.traindata();
In [35]:
train_x = Float64.(reshape(train_x, (d, n)));
In [36]:
train_x_mean = mean(train_x; dims=2);
In [37]:
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 [38]:
Sigma = (train_x_minus_mean * train_x_minus_mean')/n;
In [39]:
# find the eigendecomposition of Sigma
ESigma = eigen(Sigma; sortby = (x -> -x))
Out[39]:
In [40]:
# largest eigenvalue
ESigma.values[1]
Out[40]:
In [41]:
# corresponding eigenvector
ESigma.vectors[:, 1]
Out[41]:
In [104]:
imshow(reshape(ESigma.vectors[:, 1], (28,28))');
In [103]:
imshow(reshape(ESigma.vectors[:, 2], (28,28))');
In [105]:
imshow(reshape(ESigma.vectors[:, 3], (28,28))');
In [106]:
imshow(reshape(ESigma.vectors[:, 4], (28,28))');
In [107]:
imshow(reshape(ESigma.vectors[:, 5], (28,28))');
In [98]:
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 [204]:
D = 5;
In [205]:
A = ESigma.vectors[:, 1:D]';
In [206]:
# original image
i = 1337;
imshow(reshape(train_x[:,i], (28,28))');
In [207]:
# 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 [ ]:
Sparsity¶
In [158]:
n = 4096;
X = sprand(n,n,0.1);
Y = sprand(n,n,0.1);
In [159]:
# visualize it
imshow(X);
In [160]:
# get the size of X in bytes
Base.summarysize(X)
Out[160]:
In [161]:
# type of X
typeof(X)
Out[161]:
In [162]:
Xdense = Matrix(X);
Ydense = Matrix(Y);
typeof(Xdense)
Out[162]:
In [163]:
Base.summarysize(Xdense)
Out[163]:
In [164]:
# what's the ratio?
Base.summarysize(X) / Base.summarysize(Xdense)
Out[164]:
In [165]:
# what fraction of X's entries are nonzero?
nnz(X) / length(X)
Out[165]:
Why are these numbers not equal?
In [170]:
@time X * Y;
In [198]:
dense_mmpy_time = @elapsed Xdense * Ydense
Out[198]:
What happened here???
In [196]:
# make X and Y less dense
Xsparser = sprand(n,n,0.01);
Ysparser = sprand(n,n,0.01);
@time Xsparser * Ysparser;
In [186]:
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 [199]:
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 [ ]: