{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"pygments_lexer":"ipython3","nbconvert_exporter":"python","version":"3.6.4","file_extension":".py","codemirror_mode":{"name":"ipython","version":3},"name":"python","mimetype":"text/x-python"}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"<h1 style='text-align: center'>Visualizing Neural Networks using t-SNE and UMAP</h1>\n</br>\nA general issue in machine learning is that we often have to trade off model interpretability for prediction accuracy. This is especially true for neural networks, where we have a very complex shape of the function $f$, that connects the input variables to the output variables. \n\nIn this notebook, I want to demonstrate two techniques that can used to visualize neural networks.\nI will demonstrate the techniques on the well-known MNIST dataset. \n\n**Outline**\n\n1. Dataset Overview\n2. MNIST Classifier\n3. t-SNE Visualization\n4. UMAP Visualization\n5. Performance Comparison","metadata":{}},{"cell_type":"code","source":"import numpy as np\nimport torch\nimport torch.nn as nn\nimport torchvision\nimport pandas as pd\nimport matplotlib.pyplot as plt\nimport seaborn as sns\nimport itertools\nimport pytorch_lightning as pl\nfrom torch.utils.data import Dataset\nfrom torch.utils.data import DataLoader\nfrom pytorch_lightning import Trainer\nfrom pytorch_lightning.callbacks.early_stopping import EarlyStopping\nimport cuml\nimport sklearn\nimport umap\nimport timeit\n\ntorch.manual_seed(42)\nnp.random.seed(42)","metadata":{"execution":{"iopub.status.busy":"2021-09-24T14:03:13.889708Z","iopub.execute_input":"2021-09-24T14:03:13.889999Z","iopub.status.idle":"2021-09-24T14:03:36.520262Z","shell.execute_reply.started":"2021-09-24T14:03:13.889930Z","shell.execute_reply":"2021-09-24T14:03:36.519388Z"},"_kg_hide-input":true,"_kg_hide-output":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# 1. Dataset Overview\n\nLet's first have a quick overview of the data we are dealing with. The MNIST dataset contains handwritten digits 0 to 9.","metadata":{}},{"cell_type":"code","source":"train = pd.read_csv('../input/digit-recognizer/train.csv')\ntest = pd.read_csv('../input/digit-recognizer/test.csv')\n\nfig, ax = plt.subplots(2, 5, figsize=(15,6))\nfor a in ax.ravel():\n    a.set(xticks=[], yticks=[])\nr = 0\nfor c, num in zip(itertools.cycle([0,1,2,3,4]), range(10)):\n    img = train.loc[train.label == num].sample()\n    img = img.iloc[:, 1:].to_numpy()\n    img = np.resize(img, (28,28))\n    ax[r, c].imshow(img, cmap='gray')\n    ax[r, c].set_xlabel(f'Actual number: {num}', fontsize=12)\n    if c == 4:\n        r += 1\nfig.suptitle('MNIST Example Digits', fontsize=20, fontweight='bold')\nplt.tight_layout()\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2021-09-24T14:03:36.522005Z","iopub.execute_input":"2021-09-24T14:03:36.522399Z","iopub.status.idle":"2021-09-24T14:03:42.006469Z","shell.execute_reply.started":"2021-09-24T14:03:36.522369Z","shell.execute_reply":"2021-09-24T14:03:42.005752Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# 2. MNIST Classifier\n\nIn this section, I will construct and train a simple CNN to predict the digit of a given image. The MNIST Classifier will consist of 2 convolutional layers and 2 fully connected layers. Since this is not an important part, I will not describe my procedure in greater detail, since this is not the reason for this notebook. You can take a look at the code to get a better understanding of the architecture.","metadata":{}},{"cell_type":"code","source":"# Create a custom Dataset https://pytorch.org/tutorials/beginner/basics/data_tutorial.html\nclass MNISTDataset(Dataset):\n    '''\n    csv_path: is the relative path to the desired csv file\n    transform: transformation to be applied to the images; default None\n    '''\n    def __init__(self, csv_path, transform=None):\n        self.data=pd.read_csv(csv_path)\n        self.transform=transform\n        \n    def __len__(self):\n        return self.data.shape[0]\n    \n    def __getitem__(self, idx):\n        img = self.data.iloc[idx, 1:].to_numpy()\n        label = self.data.iloc[idx, 0]\n        return (img, label)","metadata":{"_kg_hide-input":true,"_kg_hide-output":true,"execution":{"iopub.status.busy":"2021-09-24T14:03:42.007812Z","iopub.execute_input":"2021-09-24T14:03:42.008120Z","iopub.status.idle":"2021-09-24T14:03:42.014987Z","shell.execute_reply.started":"2021-09-24T14:03:42.008085Z","shell.execute_reply":"2021-09-24T14:03:42.013990Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"class MNISTClassifier(pl.LightningModule):\n    def __init__(self, hparams):\n        super(MNISTClassifier, self).__init__()\n        self.hparams = hparams\n        self.input_size = 28*28\n        \n        self.conv = nn.Sequential(\n            nn.Conv2d(in_channels=1, out_channels=6, kernel_size=5), \n            nn.AvgPool2d(2), \n            nn.ReLU(),   \n            nn.Conv2d(in_channels=6, out_channels=16, kernel_size=5), \n            nn.MaxPool2d(2), \n            nn.ReLU(), \n            nn.Flatten()\n            )\n        self.fc1 = nn.Linear(256, 120)\n        self.fc2 = nn.Linear(120, 84)\n        self.fc3 = nn.Linear(84, 10)\n        self.relu = nn.ReLU()\n        \n    def forward(self, x):\n        out = self.conv(x)\n        out = self.fc1(out)\n        out = self.relu(out)\n        out = self.fc2(out)\n        out = self.relu(out)\n        out = self.fc3(out)\n        return out\n        \n    def training_step(self, batch, batch_idx):\n        images, labels = batch\n        images = images.reshape(-1,1,28, 28).float()        \n        \n        outputs = self.forward(images)\n        criterion = nn.CrossEntropyLoss()\n        loss = criterion(outputs, labels)\n        \n        return loss\n        \n    def validation_step(self, batch, batch_idx):\n        images, labels = batch\n        images = images.reshape(-1,1,28, 28).float() \n        \n        outputs = self.forward(images)\n        criterion = nn.CrossEntropyLoss()\n        loss = criterion(outputs, labels)\n        \n        self.log('val_loss', loss, on_step=False, on_epoch=True, prog_bar=True, logger=True)\n        \n#     def test_step(self, batch, batch_idx):\n#         images, labels = batch\n#         images = images.reshape(-1,1,28, 28).float() \n        \n#         prediction = self.forward(images)\n#         accuracy = torch.sum(prediction == labels).item() / (len(prediction) * 1.0)\n#         return accuracy\n    \n#     def test_epoch_end(self, step_accuracy):\n#         acc = 0\n#         for i in step_accuracy:\n#             acc += i.mean() / len(self.dataset['test'])\n#         print(f'Test accuracy: {acc:.3f}')\n    \n    def configure_optimizers(self):\n        return torch.optim.Adam(self.parameters(), lr=self.hparams['lr'])\n    \n    def prepare_data(self):\n        train_path = '../input/digit-recognizer/train.csv'\n        test_path = '../input/digit-recognizer/test.csv'\n        \n        mean = train.iloc[:, 1:].mean().mean()\n        std = train.iloc[:, 1:].std().std()\n        \n        transform = torchvision.transforms.Compose([\n            torchvision.transforms.ToTensor(),\n            torchvision.transforms.Normalize(mean, std)\n        ])\n        \n        dataset = MNISTDataset(train_path, transform=transform)\n        train_data, val_data = torch.utils.data.random_split(dataset, [33600, 8400]) #80/20 Split\n        test_data = MNISTDataset(test_path, transform=transform)\n        \n        self.dataset = {}\n        self.dataset['train'], self.dataset['val'], self.dataset['test'] = train_data, val_data, test_data\n        \n    def train_dataloader(self):\n        return DataLoader(self.dataset['train'], shuffle=True, batch_size=self.hparams['batch_size'], num_workers=4)\n\n    def val_dataloader(self):\n        return DataLoader(self.dataset['val'], shuffle=False, batch_size=self.hparams['batch_size'], num_workers=4)\n        \n    def test_dataloader(self):\n        return DataLoader(self.dataset['test'], shuffle=False, batch_size=self.hparams['batch_size'], num_workers=4)","metadata":{"_kg_hide-input":true,"_kg_hide-output":true,"execution":{"iopub.status.busy":"2021-09-24T14:03:42.016772Z","iopub.execute_input":"2021-09-24T14:03:42.017188Z","iopub.status.idle":"2021-09-24T14:03:42.037430Z","shell.execute_reply.started":"2021-09-24T14:03:42.017149Z","shell.execute_reply":"2021-09-24T14:03:42.036530Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"hparams = {\"lr\": 1e-3,\n           \"batch_size\": 32}\n\nmodel = MNISTClassifier(hparams)\n\ntrainer = Trainer(\n                  progress_bar_refresh_rate=50,\n                  callbacks=[EarlyStopping(monitor='val_loss', patience=3)],\n                  gpus= 1 if torch.cuda.is_available() else 0\n                 )\n\ntrainer.fit(model)","metadata":{"_kg_hide-input":true,"_kg_hide-output":true,"execution":{"iopub.status.busy":"2021-09-24T14:05:43.927616Z","iopub.execute_input":"2021-09-24T14:05:43.929711Z","iopub.status.idle":"2021-09-24T14:07:22.630751Z","shell.execute_reply.started":"2021-09-24T14:05:43.929665Z","shell.execute_reply":"2021-09-24T14:07:22.625573Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# 3. t-SNE Visualization\n\nNow that we have our model fully trained, we can start with the t-SNE visualization. t-SNE stands for t-distributed stochastic neighborhood embedding\nand tries to map high-dimensional data into lower dimensions while preserving both the local and global neighborhood structure. Details of the exact algorithm can be found in the paper [1].\n\nThe algorithm has many parameters that can be tuned, but the most important one is called $\\textit{perplexity}$. It is used to balance the attention between local or global structure. It can be seen as a guess of how many close neighbors one point has in the high dimensional space.\nA very good study of t-SNE, the hyperparameters and pitfalls can be found on [distill.pub](https://distill.pub/2016/misread-tsne/).\n\nBut now to the useful part. t-SNE gives us a tool to visualize high dimensional data. This high-dimensional data can for example be found after the CNN layers in our neural network. Both fully connected layers are very high in dimension (120-and 84-dimensional respectively). So when visualizing this layer, we can see how well the network has already been able to separate the digits in our dataset. To obtain the embeddings, we need to propagate the images through the network, but not until the final layer, but before, at the fully connected layer of interest. The following function does exactly that, capture the output after the second fully connected layer.","metadata":{}},{"cell_type":"code","source":"def embedding_fc2(data):\n    model.eval()\n    img = torch.Tensor(data.reshape(-1,1,28, 28)).float() #transform the image\n    out = model.conv(img) #pass through CNN\n    out = model.relu(model.fc1(out)) #pass through FC 1\n    out = model.relu(model.fc2(out)) #pass through FC 2\n    return out","metadata":{"execution":{"iopub.status.busy":"2021-09-24T14:07:22.633105Z","iopub.execute_input":"2021-09-24T14:07:22.633759Z","iopub.status.idle":"2021-09-24T14:07:22.644700Z","shell.execute_reply.started":"2021-09-24T14:07:22.633718Z","shell.execute_reply":"2021-09-24T14:07:22.643850Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"The obtained embedding has to be fit and transformed using the t-SNE algorithm. [Sklearn](https://scikit-learn.org/stable/modules/generated/sklearn.manifold.TSNE.html) has a function that performs the algorithm on the CPU. However, it takes very long to converge as the number of data points becomes larger. \nThat is why we are going to use an alternative that runs on the GPU - RAPIDS [cuML](https://docs.rapids.ai/api/cuml/stable/api.html?highlight=tsne#cuml.TSNE) library. I am going to compare the performance in section 6.\n\nWe are going to use the cuML library in the following. The three hyperparameters that I have set are: \n- n_components: dimensionality of reduction\n- perplexity: focus of attention on local or global structure (typically between 5 and 50)\n- n_neighbors: the amout of nearest neighbors considered when doing the reduction (3\\*perplexity is recommended in the documentation)\n\nTo get a better sense of the perplexity parameter, I am going the visualize the embeddings for 3 different perplexity levels: 5, 30 and 50.","metadata":{"_kg_hide-output":true}},{"cell_type":"code","source":"fig, ax = plt.subplots(1, 3, figsize=(20,7), constrained_layout=True)\n\nimg = train.iloc[:, 1:].to_numpy() #all images\nlabel = train.iloc[:, 0] #all labels\n\nembedding = embedding_fc2(img)\n\nfor c, per in zip(itertools.count(), [5, 30, 50]):\n    tsne = cuml.manifold.TSNE(n_components=2,\n                perplexity=per,\n                n_neighbors=per*3)\n    tsne = tsne.fit_transform(embedding.detach().numpy())\n    scatter = ax[c].scatter(tsne[:, 0], tsne[:, 1], c=label, cmap='tab10', s=0.3)\n    ax[c].set_title(f'Perplexity: {per}', fontsize=16)    \n\nfig.suptitle('t-SNE Dimensionality reduction', fontweight='bold', fontsize=25)\ncbar = fig.colorbar(scatter, boundaries=np.arange(11)-0.5, location='right')\ncbar.set_ticks(np.arange(10))\ncbar.set_ticklabels(np.arange(10))\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2021-09-24T14:28:56.340390Z","iopub.execute_input":"2021-09-24T14:28:56.340758Z","iopub.status.idle":"2021-09-24T14:29:09.293487Z","shell.execute_reply.started":"2021-09-24T14:28:56.340718Z","shell.execute_reply":"2021-09-24T14:29:09.292431Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"We can see that across all the three different perplexity levels, the data is well structured and separated. That indicates that our network can successfully separate between the different digits. I want to mention that those three t-SNE results represent the same underlying data, so even though the right-most scatterplot seems to be better, in reality that has just to do with the way the dimensionality reduction has been constructed.","metadata":{}},{"cell_type":"markdown","source":"# 4. UMAP Visualization\n\nAn alternative to t-SNE is the UMAP dimensionality reduction algorithm. Uniform Manifold Approximation and Projection aims to preserve the manifold of the high dimensional data when mapping into the lower-dimensional representation. For details, please take a look at the paper ([2]). While t-SNE is still to most popular algorithm, UMAP turns out to be a very good alternative. What suprised me about it is that it can run quickly even on CPU, whereas t-SNE may take up to an several minutes. I am going to compare performance in the next section.\n\nNow to the most important parameters of UMAP:\n- n_components: dimensionality of reduction\n- n_neighbors: tradeoff between local and global structure (lower values only take into consideration the local structure)\n- min_dist: the minimum allowed distance in the lower structure; can prevent or encourage forming of clusters\n- metric: what distance measure to use\n\nAs we can see, both algorithms share some hyperparameters, which is no suprise as they both are based on neighborhood structure.\nYou might want to also take a look at the [docs](https://umap-learn.readthedocs.io/en/latest/), they offer some great tutorials.","metadata":{}},{"cell_type":"code","source":"umap_hparams = {'n_neighbors':5,\n                'min_dist':0.1,\n                'n_components':2,\n                'metric':'euclidean'}\n\nimg = train.iloc[:, 1:].to_numpy() #all images\nlabel = train.iloc[:, 0] #all labels\nembedding = embedding_fc2(img)\n\nfig, ax = plt.subplots(figsize=(11,8), constrained_layout=False)\nax.set(xticks=[], yticks=[])\n\numap_embedding = umap.UMAP(n_neighbors=umap_hparams['n_neighbors'], min_dist=umap_hparams['min_dist'], n_components=umap_hparams['n_components'], metric=umap_hparams['metric']).fit_transform(embedding.detach().numpy())\nscatter = ax.scatter(x = umap_embedding[:,0], y = umap_embedding[:,1], s=0.3, c=label, cmap='tab10')\n\ncbar = plt.colorbar(scatter, boundaries=np.arange(11)-0.5)\ncbar.set_ticks(np.arange(10))\ncbar.set_ticklabels(np.arange(10))\n\nplt.title('UMAP Dimensionality reduction', fontsize=25, fontweight='bold')\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2021-09-24T14:28:25.724810Z","iopub.execute_input":"2021-09-24T14:28:25.725177Z","iopub.status.idle":"2021-09-24T14:28:56.338659Z","shell.execute_reply.started":"2021-09-24T14:28:25.725145Z","shell.execute_reply":"2021-09-24T14:28:56.337921Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# 5. Performance Comparison\n\nNow I am going to compare the three libraries Sklearn, cuML and umap based on performance. Obviously it is hard to compare t-SNE to UMAP as they are build on different mathematical foundation, but I still find is interesting to find out which algorithms performs better (especially since both are a great tool for visualizing high-dimensional data). I am going to time the algorithms on a randomly created 64-dimensional dataset. The time is tracked across different dataset sizes, ranging from 100 to 51200. I am going to use default settings for all algorithms/versions.","metadata":{}},{"cell_type":"code","source":"sizes = [100, 200, 400, 800, 1600, 3200, 6400, 12800, 25600, 51200]\ndimension = 64\n\ntsne_cpu_ = []\ntsne_gpu_ = []\numap_ = []\n\nfor size in sizes:\n    data = np.random.rand(size, dimension)\n\n    a = timeit.default_timer()\n    sklearn.manifold.TSNE().fit_transform(data)\n    b = timeit.default_timer()\n    tsne_cpu_.append(b-a)\n\nfor size in sizes:\n    data = np.random.rand(size, dimension)\n    \n    a = timeit.default_timer()\n    cuml.manifold.TSNE().fit_transform(data)\n    b = timeit.default_timer()\n    tsne_gpu_.append(b-a)\n    \nfor size in sizes:\n    data = np.random.rand(size, dimension)\n    \n    a = timeit.default_timer()\n    umap.UMAP().fit_transform(data)\n    b = timeit.default_timer()\n    umap_.append(b-a)","metadata":{"execution":{"iopub.status.busy":"2021-09-24T14:09:13.815605Z","iopub.execute_input":"2021-09-24T14:09:13.815988Z","iopub.status.idle":"2021-09-24T14:24:22.565323Z","shell.execute_reply.started":"2021-09-24T14:09:13.815954Z","shell.execute_reply":"2021-09-24T14:24:22.564096Z"},"_kg_hide-input":true,"_kg_hide-output":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import numpy as np\nimport matplotlib.pyplot as plt\nloc = np.arange(len(sizes))\nwidth = 0.2\n\nfig, ax = plt.subplots(figsize=(12,8))\n\nr1 = ax.bar(loc-width, tsne_gpu_, width)\nr2 = ax.bar(loc, umap_, width)\nr3 = ax.bar(loc+width, tsne_cpu_, width)\nax.legend((r1[0], r2[0], r3[0]), ('t-SNE on GPU (cuML)', 'UMAP', 't-SNE on CPU (Sklearn)'), loc='upper left', fontsize=12)\n\nax.set(xticks=loc, xticklabels=sizes)\nax.yaxis.set_major_formatter('{x:.0f}s')\nax.set_ylabel('Time in Seconds', fontsize=12)\nax.set_xlabel('Data Size', fontsize=12)\n\nplt.title(f'Performance Comparison on {dimension}-dimensional data', fontweight='bold', fontsize=20)\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2021-09-24T14:24:29.533030Z","iopub.execute_input":"2021-09-24T14:24:29.533363Z","iopub.status.idle":"2021-09-24T14:24:29.768123Z","shell.execute_reply.started":"2021-09-24T14:24:29.533333Z","shell.execute_reply":"2021-09-24T14:24:29.767132Z"},"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"It is not surprising that the t-SNE algorithm running on GPU outperforms both CPU variations. However we can clearly see that UMAP scales better to large data in terms of time than Sklearn's t-SNE. Finally Sklearn's t-SNE has some weird variation in the graph which I cannot explain. I have timed it several times and there always seemed to be a peek in between the two extremes.\n\nOverall, both algorithms offer some great insight into the complex data we are dealing with, and how well the model has been able to learn to segment different images.","metadata":{}},{"cell_type":"markdown","source":"# 7. Resources\n\n[1] L. van der Maaten, G. Hinton; Visualizing Data using t-SNE; https://www.cs.toronto.edu/~hinton/absps/tsne.pdf\n\n[2] L. McInnes, J. Healy, and J. Melville; UMAP: Uniform Manifold Approximation and Projection for Dimension Reduction; http://arxiv.org/abs/1802.03426","metadata":{}},{"cell_type":"markdown","source":"<h2 style='text-align: center'>That you for reading this notebook to the end!<br>Feel free to upvote and leave a comment.</h2><h4 style='text-align: center'>Also please let me know what can be improved...<h4>","metadata":{}}]}