import numpy as np
import matplotlib.pyplot as plt
import torch
import torch.nn as nn
import torch.optim as optim
%matplotlib inline
# Retina
%config InlineBackend.figure_format = 'retina'
PCA
ML
= 500
N = 2
D
1)
torch.manual_seed(= torch.distributions.MultivariateNormal(torch.zeros(D), torch.tensor([[2.0, 0.5], [0.5, 0.2]])).sample((N,))
X
= X - X.mean(dim=0)
X 0], X[:,1])
plt.scatter(X[:,'equal', adjustable='box') plt.gca().set_aspect(
= torch.matmul(X.t(), X) / (X.size(0)-1)
covariance_matrix covariance_matrix.shape
torch.Size([2, 2])
covariance_matrix
tensor([[1.9740, 0.5022],
[0.5022, 0.2127]])
= torch.linalg.eig(covariance_matrix) eigenvalues, eigenvectors
= eigenvalues.real
eigenvalues eigenvalues
tensor([2.1071, 0.0796])
= eigenvectors.real
eigenvectors eigenvectors
tensor([[ 0.9666, -0.2562],
[ 0.2562, 0.9666]])
# Plot the eigenvectors along with length of eigenvalues
0], X[:,1])
plt.scatter(X[:,'equal', adjustable='box')
plt.gca().set_aspect(0, 0, eigenvectors[0, 0], eigenvectors[1, 0], scale=eigenvalues[0].sqrt().item(), color='red', label='Eigenvector 1')
plt.quiver(0, 0, eigenvectors[0, 1], eigenvectors[1, 1], scale=eigenvalues[1].sqrt().item(), color='green', label='Eigenvector 2')
plt.quiver( plt.legend()
= 1
k = torch.matmul(X, eigenvectors[:, :k])
Z Z.shape
torch.Size([500, 1])
'o') plt.plot(Z, torch.zeros_like(Z),
# Step 4: Sort eigenvectors based on eigenvalues
= torch.argsort(eigenvalues[:, 0], descending=True)
sorted_indices = eigenvectors[:, sorted_indices] sorted_eigenvectors
--------------------------------------------------------------------------- IndexError Traceback (most recent call last) Cell In[93], line 2 1 # Step 4: Sort eigenvectors based on eigenvalues ----> 2 sorted_indices = torch.argsort(eigenvalues[:, 0], descending=True) 3 sorted_eigenvectors = eigenvectors[:, sorted_indices] IndexError: too many indices for tensor of dimension 1
sorted
=0), X.var(dim=0) X.mean(dim
(tensor([-5.7220e-09, 0.0000e+00]), tensor([1.9740, 0.2127]))
= 1 # Number of principal components
k = nn.Parameter(torch.randn(D, k)) B
B.shape, X.shape
(torch.Size([2, 1]), torch.Size([500, 2]))
= nn.MSELoss()
mse_loss
= optim.Adam([B], lr=0.001)
optimizer = 2000 # Number of optimization epochs
num_epochs
for epoch in range(num_epochs):
optimizer.zero_grad()
# Project data onto the current projection matrix
= torch.matmul(X, B)
z = torch.matmul(z, B.t())
x_reconstructed
# Compute reconstruction loss
= mse_loss(x_reconstructed, X)
loss
# Perform backpropagation and optimization step
loss.backward()
optimizer.step()
if (epoch + 1) % 100 == 0:
print(f'Epoch [{epoch + 1}/{num_epochs}], Loss: {loss.item()}')
Epoch [100/2000], Loss: 0.7209210395812988
Epoch [200/2000], Loss: 0.6048457026481628
Epoch [300/2000], Loss: 0.5128658413887024
Epoch [400/2000], Loss: 0.42206627130508423
Epoch [500/2000], Loss: 0.3365266025066376
Epoch [600/2000], Loss: 0.26167812943458557
Epoch [700/2000], Loss: 0.2000642865896225
Epoch [800/2000], Loss: 0.1518339067697525
Epoch [900/2000], Loss: 0.11567910760641098
Epoch [1000/2000], Loss: 0.08960824459791183
Epoch [1100/2000], Loss: 0.07147730886936188
Epoch [1200/2000], Loss: 0.059300925582647324
Epoch [1300/2000], Loss: 0.0514017716050148
Epoch [1400/2000], Loss: 0.0464538112282753
Epoch [1500/2000], Loss: 0.04346401244401932
Epoch [1600/2000], Loss: 0.04172360897064209
Epoch [1700/2000], Loss: 0.0407492071390152
Epoch [1800/2000], Loss: 0.040225498378276825
Epoch [1900/2000], Loss: 0.03995582088828087
Epoch [2000/2000], Loss: 0.039823081344366074
# Project data onto the learned projection matrix
= torch.matmul(X, B) Z
B
Parameter containing:
tensor([[-0.9645],
[-0.2663]], requires_grad=True)
def train_linear_autoencoder(X, K=2, lr=0.001, num_epochs=2000):
= X.size(1)
D = nn.Parameter(torch.randn(D, K))
B
= nn.MSELoss()
mse_loss = optim.Adam([B], lr=lr)
optimizer = num_epochs
num_epochs
for epoch in range(num_epochs):
optimizer.zero_grad()
= torch.matmul(X, B)
z = torch.sigmoid(torch.matmul(z, B.t()))
x_reconstructed
= mse_loss(x_reconstructed, X)
loss
loss.backward()
optimizer.step()
if (epoch + 1) % 100 == 0:
print(f'Epoch [{epoch + 1}/{num_epochs}], Loss: {loss.item()}')
return B
# Get MNIST data
from torchvision import datasets
from torchvision.transforms import ToTensor
= datasets.MNIST(root='~/.datasets', train=True, transform=ToTensor(), download=True) mnist_train
= mnist_train.data.float()
X = X.view(X.size(0), -1)
X = X / 255.0
X
X.shape
torch.Size([60000, 784])
# Work with first 1000 samples
= X[:3000] X
min(), X.max() X.
(tensor(0.), tensor(1.))
# Plot the first 10 images
= plt.subplots(1, 10, figsize=(10, 1))
fig, ax for i in range(10):
28, 28).numpy(), cmap='gray')
ax[i].imshow(X[i].view('off') ax[i].axis(
= train_linear_autoencoder(X, K=5, lr=0.001, num_epochs=2000) B
Epoch [100/2000], Loss: 0.4135352075099945
Epoch [200/2000], Loss: 0.3644017279148102
Epoch [300/2000], Loss: 0.3089104890823364
Epoch [400/2000], Loss: 0.26013466715812683
Epoch [500/2000], Loss: 0.22550225257873535
Epoch [600/2000], Loss: 0.20584215223789215
Epoch [700/2000], Loss: 0.19726672768592834
Epoch [800/2000], Loss: 0.1943318247795105
Epoch [900/2000], Loss: 0.1914585679769516
Epoch [1000/2000], Loss: 0.18916504085063934
Epoch [1100/2000], Loss: 0.18805715441703796
Epoch [1200/2000], Loss: 0.18697765469551086
Epoch [1300/2000], Loss: 0.18595127761363983
Epoch [1400/2000], Loss: 0.18386924266815186
Epoch [1500/2000], Loss: 0.1826687455177307
Epoch [1600/2000], Loss: 0.18168936669826508
Epoch [1700/2000], Loss: 0.18077248334884644
Epoch [1800/2000], Loss: 0.17985987663269043
Epoch [1900/2000], Loss: 0.17898057401180267
Epoch [2000/2000], Loss: 0.17811761796474457
B.shape
torch.Size([784, 5])
# Project data onto the 2d space
with torch.no_grad():
= torch.matmul(X, B)
Z = torch.sigmoid(torch.matmul(Z, B.t()))
X_reconstructed
# Plot the first 10 images and their reconstructions
= plt.subplots(2, 10, figsize=(10, 2))
fig, ax for i in range(10):
0, i].imshow(X[i].view(28, 28).numpy(), cmap='gray')
ax[0, i].axis('off')
ax[1, i].imshow(X_reconstructed[i].view(28, 28).numpy(), cmap='gray')
ax[1, i].axis('off') ax[
from sklearn.decomposition import PCA
= PCA(n_components=32)
pca
= pca.fit_transform(X.numpy())
Z_pca
= pca.inverse_transform(Z_pca)
X_reconstructed_pca
# Plot the first 10 images and their reconstructions
= plt.subplots(2, 10, figsize=(10, 2))
fig, ax for i in range(10):
0, i].imshow(X[i].view(28, 28).numpy(), cmap='gray')
ax[0, i].axis('off')
ax[1, i].imshow(X_reconstructed_pca[i].reshape(28, 28), cmap='gray')
ax[1, i].axis('off') ax[
import pandas as pd
= pd.DataFrame({"X": X[0], "X_reconstructed": X_reconstructed[0], "X_reconstructed_pca": X_reconstructed_pca[0]}) df
df.plot()
# Plot the original data and the projected data
=(6, 6))
plt.figure(figsize0], X[:,1], label='Original data')
plt.scatter(X[:,0], Z[:,0], label='Projected data') plt.scatter(Z[:,
--------------------------------------------------------------------------- RuntimeError Traceback (most recent call last) Cell In[25], line 4 2 plt.figure(figsize=(6, 6)) 3 plt.scatter(X[:,0], X[:,1], label='Original data') ----> 4 plt.scatter(Z[:,0], Z[:,0], label='Projected data') File ~/miniconda3/lib/python3.9/site-packages/matplotlib/pyplot.py:3684, in scatter(x, y, s, c, marker, cmap, norm, vmin, vmax, alpha, linewidths, edgecolors, plotnonfinite, data, **kwargs) 3665 @_copy_docstring_and_deprecators(Axes.scatter) 3666 def scatter( 3667 x: float | ArrayLike, (...) 3682 **kwargs, 3683 ) -> PathCollection: -> 3684 __ret = gca().scatter( 3685 x, 3686 y, 3687 s=s, 3688 c=c, 3689 marker=marker, 3690 cmap=cmap, 3691 norm=norm, 3692 vmin=vmin, 3693 vmax=vmax, 3694 alpha=alpha, 3695 linewidths=linewidths, 3696 edgecolors=edgecolors, 3697 plotnonfinite=plotnonfinite, 3698 **({"data": data} if data is not None else {}), 3699 **kwargs, 3700 ) 3701 sci(__ret) 3702 return __ret File ~/miniconda3/lib/python3.9/site-packages/matplotlib/__init__.py:1465, in _preprocess_data.<locals>.inner(ax, data, *args, **kwargs) 1462 @functools.wraps(func) 1463 def inner(ax, *args, data=None, **kwargs): 1464 if data is None: -> 1465 return func(ax, *map(sanitize_sequence, args), **kwargs) 1467 bound = new_sig.bind(ax, *args, **kwargs) 1468 auto_label = (bound.arguments.get(label_namer) 1469 or bound.kwargs.get(label_namer)) File ~/miniconda3/lib/python3.9/site-packages/matplotlib/axes/_axes.py:4649, in Axes.scatter(self, x, y, s, c, marker, cmap, norm, vmin, vmax, alpha, linewidths, edgecolors, plotnonfinite, **kwargs) 4646 x, y = self._process_unit_info([("x", x), ("y", y)], kwargs) 4647 # np.ma.ravel yields an ndarray, not a masked array, 4648 # unless its argument is a masked array. -> 4649 x = np.ma.ravel(x) 4650 y = np.ma.ravel(y) 4651 if x.size != y.size: File ~/miniconda3/lib/python3.9/site-packages/numpy/ma/core.py:6773, in _frommethod.__call__(self, a, *args, **params) 6770 args = list(args) 6771 a, args[0] = args[0], a -> 6773 marr = asanyarray(a) 6774 method_name = self.__name__ 6775 method = getattr(type(marr), method_name, None) File ~/miniconda3/lib/python3.9/site-packages/numpy/ma/core.py:8005, in asanyarray(a, dtype) 8003 if isinstance(a, MaskedArray) and (dtype is None or dtype == a.dtype): 8004 return a -> 8005 return masked_array(a, dtype=dtype, copy=False, keep_mask=True, subok=True) File ~/miniconda3/lib/python3.9/site-packages/numpy/ma/core.py:2826, in MaskedArray.__new__(cls, data, mask, dtype, copy, subok, ndmin, fill_value, keep_mask, hard_mask, shrink, order) 2817 """ 2818 Create a new masked array from scratch. 2819 (...) 2823 2824 """ 2825 # Process data. -> 2826 _data = np.array(data, dtype=dtype, copy=copy, 2827 order=order, subok=True, ndmin=ndmin) 2828 _baseclass = getattr(data, '_baseclass', type(_data)) 2829 # Check that we're not erasing the mask. File ~/miniconda3/lib/python3.9/site-packages/torch/_tensor.py:1030, in Tensor.__array__(self, dtype) 1028 return handle_torch_function(Tensor.__array__, (self,), self, dtype=dtype) 1029 if dtype is None: -> 1030 return self.numpy() 1031 else: 1032 return self.numpy().astype(dtype, copy=False) RuntimeError: Can't call numpy() on Tensor that requires grad. Use tensor.detach().numpy() instead.