04_sparse-matrices.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img src=\"dsci512_header.png\" width=\"600\">\n",
    "\n",
    "# Lecture 4 Sparse matrices"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Topics\n",
    "- Dask and parallelized computation on CPUs \n",
    "- Pytorch and GPU computation\n",
    "- Data types\n",
    "- Sparse matrices\n",
    "- Bag-of-word matrices"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Learning objectives\n",
    "\n",
    "- Identify situations where matrix parallelization is necessary\n",
    "- Understand the date types\n",
    "- Understand sparse matrix representations\n",
    "- Understand bag-of-word matrices"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Block matrix operations\n",
    "\n",
    "A block matrix is a matrix that is divided into smaller rectangular or square submatrices, which are called blocks.\n",
    "\n",
    "When doing matrix multiplication, you can split a large matrix into smaller blocks and perform operations on these blocks independently. You can get the exact same result as doing multiplication with the whole matrix. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "A^2 computed using block operations:\n",
      " [[ 90 100 110 120]\n",
      " [202 228 254 280]\n",
      " [314 356 398 440]\n",
      " [426 484 542 600]]\n",
      "\n",
      "A^2 computed using standard matrix multiplication:\n",
      " [[ 90 100 110 120]\n",
      " [202 228 254 280]\n",
      " [314 356 398 440]\n",
      " [426 484 542 600]]\n",
      "\n",
      "Are the results equivalent? True\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "\n",
    "# Create a large matrix A\n",
    "A = np.array([[1, 2, 3, 4],\n",
    "              [5, 6, 7, 8],\n",
    "              [9, 10, 11, 12],\n",
    "              [13, 14, 15, 16]])\n",
    "\n",
    "# Partition A into four blocks\n",
    "A11 = A[:2, :2]\n",
    "A12 = A[:2, 2:]\n",
    "A21 = A[2:, :2]\n",
    "A22 = A[2:, 2:]\n",
    "\n",
    "# Compute A^2 using block matrix operations\n",
    "A_squared_block = np.block([[np.dot(A11, A11) + np.dot(A12, A21), np.dot(A11, A12) + np.dot(A12, A22)],\n",
    "                            [np.dot(A21, A11) + np.dot(A22, A21), np.dot(A21, A12) + np.dot(A22, A22)]])\n",
    "\n",
    "# Compute A^2 using standard matrix multiplication\n",
    "A_squared = np.dot(A, A)\n",
    "\n",
    "# Check if the two results are equivalent\n",
    "equivalent = np.array_equal(A_squared, A_squared_block)\n",
    "\n",
    "print(\"A^2 computed using block operations:\\n\", A_squared_block)\n",
    "print(\"\\nA^2 computed using standard matrix multiplication:\\n\", A_squared)\n",
    "print(\"\\nAre the results equivalent?\", equivalent)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can utilize this nice property to speed up matrix computation. We can do operations on individual blocls means we can parallelize the operations, that is, distribute the task to individual processes. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[Dask](https://www.dask.org/) is a flexible parallel computing library for Python that is designed to scale computations across multiple cores or even distributed systems.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To use dask, we first need to start the client. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We need to Start the Client. This is just for setting up dask. You need to do it once at the start of your program and ignore it afterwards. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "    <div style=\"width: 24px; height: 24px; background-color: #e1e1e1; border: 3px solid #9D9D9D; border-radius: 5px; position: absolute;\"> </div>\n",
       "    <div style=\"margin-left: 48px;\">\n",
       "        <h3 style=\"margin-bottom: 0px;\">Client</h3>\n",
       "        <p style=\"color: #9D9D9D; margin-bottom: 0px;\">Client-96e4d429-8b86-11ef-a4bb-cb798e5c5873</p>\n",
       "        <table style=\"width: 100%; text-align: left;\">\n",
       "\n",
       "        <tr>\n",
       "        \n",
       "            <td style=\"text-align: left;\"><strong>Connection method:</strong> Cluster object</td>\n",
       "            <td style=\"text-align: left;\"><strong>Cluster type:</strong> distributed.LocalCluster</td>\n",
       "        \n",
       "        </tr>\n",
       "\n",
       "        \n",
       "            <tr>\n",
       "                <td style=\"text-align: left;\">\n",
       "                    <strong>Dashboard: </strong> <a href=\"http://127.0.0.1:8787/status\" target=\"_blank\">http://127.0.0.1:8787/status</a>\n",
       "                </td>\n",
       "                <td style=\"text-align: left;\"></td>\n",
       "            </tr>\n",
       "        \n",
       "\n",
       "        </table>\n",
       "\n",
       "        \n",
       "\n",
       "        \n",
       "            <details>\n",
       "            <summary style=\"margin-bottom: 20px;\"><h3 style=\"display: inline;\">Cluster Info</h3></summary>\n",
       "            <div class=\"jp-RenderedHTMLCommon jp-RenderedHTML jp-mod-trusted jp-OutputArea-output\">\n",
       "    <div style=\"width: 24px; height: 24px; background-color: #e1e1e1; border: 3px solid #9D9D9D; border-radius: 5px; position: absolute;\">\n",
       "    </div>\n",
       "    <div style=\"margin-left: 48px;\">\n",
       "        <h3 style=\"margin-bottom: 0px; margin-top: 0px;\">LocalCluster</h3>\n",
       "        <p style=\"color: #9D9D9D; margin-bottom: 0px;\">a3ed19ab</p>\n",
       "        <table style=\"width: 100%; text-align: left;\">\n",
       "            <tr>\n",
       "                <td style=\"text-align: left;\">\n",
       "                    <strong>Dashboard:</strong> <a href=\"http://127.0.0.1:8787/status\" target=\"_blank\">http://127.0.0.1:8787/status</a>\n",
       "                </td>\n",
       "                <td style=\"text-align: left;\">\n",
       "                    <strong>Workers:</strong> 2\n",
       "                </td>\n",
       "            </tr>\n",
       "            <tr>\n",
       "                <td style=\"text-align: left;\">\n",
       "                    <strong>Total threads:</strong> 6\n",
       "                </td>\n",
       "                <td style=\"text-align: left;\">\n",
       "                    <strong>Total memory:</strong> 15.56 GiB\n",
       "                </td>\n",
       "            </tr>\n",
       "            \n",
       "            <tr>\n",
       "    <td style=\"text-align: left;\"><strong>Status:</strong> running</td>\n",
       "    <td style=\"text-align: left;\"><strong>Using processes:</strong> True</td>\n",
       "</tr>\n",
       "\n",
       "            \n",
       "        </table>\n",
       "\n",
       "        <details>\n",
       "            <summary style=\"margin-bottom: 20px;\">\n",
       "                <h3 style=\"display: inline;\">Scheduler Info</h3>\n",
       "            </summary>\n",
       "\n",
       "            <div style=\"\">\n",
       "    <div>\n",
       "        <div style=\"width: 24px; height: 24px; background-color: #FFF7E5; border: 3px solid #FF6132; border-radius: 5px; position: absolute;\"> </div>\n",
       "        <div style=\"margin-left: 48px;\">\n",
       "            <h3 style=\"margin-bottom: 0px;\">Scheduler</h3>\n",
       "            <p style=\"color: #9D9D9D; margin-bottom: 0px;\">Scheduler-edcdf493-87e8-4cc8-9690-6971f1b56c25</p>\n",
       "            <table style=\"width: 100%; text-align: left;\">\n",
       "                <tr>\n",
       "                    <td style=\"text-align: left;\">\n",
       "                        <strong>Comm:</strong> tcp://127.0.0.1:40135\n",
       "                    </td>\n",
       "                    <td style=\"text-align: left;\">\n",
       "                        <strong>Workers:</strong> 2\n",
       "                    </td>\n",
       "                </tr>\n",
       "                <tr>\n",
       "                    <td style=\"text-align: left;\">\n",
       "                        <strong>Dashboard:</strong> <a href=\"http://127.0.0.1:8787/status\" target=\"_blank\">http://127.0.0.1:8787/status</a>\n",
       "                    </td>\n",
       "                    <td style=\"text-align: left;\">\n",
       "                        <strong>Total threads:</strong> 6\n",
       "                    </td>\n",
       "                </tr>\n",
       "                <tr>\n",
       "                    <td style=\"text-align: left;\">\n",
       "                        <strong>Started:</strong> Just now\n",
       "                    </td>\n",
       "                    <td style=\"text-align: left;\">\n",
       "                        <strong>Total memory:</strong> 15.56 GiB\n",
       "                    </td>\n",
       "                </tr>\n",
       "            </table>\n",
       "        </div>\n",
       "    </div>\n",
       "\n",
       "    <details style=\"margin-left: 48px;\">\n",
       "        <summary style=\"margin-bottom: 20px;\">\n",
       "            <h3 style=\"display: inline;\">Workers</h3>\n",
       "        </summary>\n",
       "\n",
       "        \n",
       "        <div style=\"margin-bottom: 20px;\">\n",
       "            <div style=\"width: 24px; height: 24px; background-color: #DBF5FF; border: 3px solid #4CC9FF; border-radius: 5px; position: absolute;\"> </div>\n",
       "            <div style=\"margin-left: 48px;\">\n",
       "            <details>\n",
       "                <summary>\n",
       "                    <h4 style=\"margin-bottom: 0px; display: inline;\">Worker: 0</h4>\n",
       "                </summary>\n",
       "                <table style=\"width: 100%; text-align: left;\">\n",
       "                    <tr>\n",
       "                        <td style=\"text-align: left;\">\n",
       "                            <strong>Comm: </strong> tcp://127.0.0.1:37213\n",
       "                        </td>\n",
       "                        <td style=\"text-align: left;\">\n",
       "                            <strong>Total threads: </strong> 3\n",
       "                        </td>\n",
       "                    </tr>\n",
       "                    <tr>\n",
       "                        <td style=\"text-align: left;\">\n",
       "                            <strong>Dashboard: </strong> <a href=\"http://127.0.0.1:41589/status\" target=\"_blank\">http://127.0.0.1:41589/status</a>\n",
       "                        </td>\n",
       "                        <td style=\"text-align: left;\">\n",
       "                            <strong>Memory: </strong> 7.78 GiB\n",
       "                        </td>\n",
       "                    </tr>\n",
       "                    <tr>\n",
       "                        <td style=\"text-align: left;\">\n",
       "                            <strong>Nanny: </strong> tcp://127.0.0.1:42535\n",
       "                        </td>\n",
       "                        <td style=\"text-align: left;\"></td>\n",
       "                    </tr>\n",
       "                    <tr>\n",
       "                        <td colspan=\"2\" style=\"text-align: left;\">\n",
       "                            <strong>Local directory: </strong> /home/lukeum/Documents/DSCI_512_alg-data-struct_instructors/dask_tmp/dask-scratch-space/worker-54pi5lst\n",
       "                        </td>\n",
       "                    </tr>\n",
       "\n",
       "                    \n",
       "                    <tr>\n",
       "                        <td style=\"text-align: left;\">\n",
       "                            <strong>GPU: </strong>NVIDIA GeForce RTX 2080 Ti\n",
       "                        </td>\n",
       "                        <td style=\"text-align: left;\">\n",
       "                            <strong>GPU memory: </strong> 11.00 GiB\n",
       "                        </td>\n",
       "                    </tr>\n",
       "                    \n",
       "\n",
       "                    \n",
       "\n",
       "                </table>\n",
       "            </details>\n",
       "            </div>\n",
       "        </div>\n",
       "        \n",
       "        <div style=\"margin-bottom: 20px;\">\n",
       "            <div style=\"width: 24px; height: 24px; background-color: #DBF5FF; border: 3px solid #4CC9FF; border-radius: 5px; position: absolute;\"> </div>\n",
       "            <div style=\"margin-left: 48px;\">\n",
       "            <details>\n",
       "                <summary>\n",
       "                    <h4 style=\"margin-bottom: 0px; display: inline;\">Worker: 1</h4>\n",
       "                </summary>\n",
       "                <table style=\"width: 100%; text-align: left;\">\n",
       "                    <tr>\n",
       "                        <td style=\"text-align: left;\">\n",
       "                            <strong>Comm: </strong> tcp://127.0.0.1:42963\n",
       "                        </td>\n",
       "                        <td style=\"text-align: left;\">\n",
       "                            <strong>Total threads: </strong> 3\n",
       "                        </td>\n",
       "                    </tr>\n",
       "                    <tr>\n",
       "                        <td style=\"text-align: left;\">\n",
       "                            <strong>Dashboard: </strong> <a href=\"http://127.0.0.1:40475/status\" target=\"_blank\">http://127.0.0.1:40475/status</a>\n",
       "                        </td>\n",
       "                        <td style=\"text-align: left;\">\n",
       "                            <strong>Memory: </strong> 7.78 GiB\n",
       "                        </td>\n",
       "                    </tr>\n",
       "                    <tr>\n",
       "                        <td style=\"text-align: left;\">\n",
       "                            <strong>Nanny: </strong> tcp://127.0.0.1:44997\n",
       "                        </td>\n",
       "                        <td style=\"text-align: left;\"></td>\n",
       "                    </tr>\n",
       "                    <tr>\n",
       "                        <td colspan=\"2\" style=\"text-align: left;\">\n",
       "                            <strong>Local directory: </strong> /home/lukeum/Documents/DSCI_512_alg-data-struct_instructors/dask_tmp/dask-scratch-space/worker-nf9wilat\n",
       "                        </td>\n",
       "                    </tr>\n",
       "\n",
       "                    \n",
       "                    <tr>\n",
       "                        <td style=\"text-align: left;\">\n",
       "                            <strong>GPU: </strong>NVIDIA GeForce RTX 2080 Ti\n",
       "                        </td>\n",
       "                        <td style=\"text-align: left;\">\n",
       "                            <strong>GPU memory: </strong> 11.00 GiB\n",
       "                        </td>\n",
       "                    </tr>\n",
       "                    \n",
       "\n",
       "                    \n",
       "\n",
       "                </table>\n",
       "            </details>\n",
       "            </div>\n",
       "        </div>\n",
       "        \n",
       "\n",
       "    </details>\n",
       "</div>\n",
       "\n",
       "        </details>\n",
       "    </div>\n",
       "</div>\n",
       "            </details>\n",
       "        \n",
       "\n",
       "    </div>\n",
       "</div>"
      ],
      "text/plain": [
       "<Client: 'tcp://127.0.0.1:40135' processes=2 threads=6, memory=15.56 GiB>"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import dask\n",
    "from dask.distributed import Client\n",
    "# set a temporary folder to store some intermediate results generated by dask.\n",
    "dask.config.set(temporary_directory='../dask_tmp')\n",
    "client = Client(n_workers=2) \n",
    "# set how many processes you want. this depends on how many cores your computer has. generally speaking, if the scale of your problem is large enough, more cores lead to better performance (assuming that your computation can be parallelized).\n",
    "client # get client information"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here is an example of numpy vs. dask. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/lukeum/anaconda3/envs/py10/lib/python3.10/site-packages/distributed/client.py:3163: UserWarning: Sending large graph of size 762.94 MiB.\n",
      "This may cause some slowdown.\n",
      "Consider scattering data ahead of time and using futures.\n",
      "  warnings.warn(\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "NumPy time: 5.186022996902466\n",
      "Dask time: 3.5584545135498047\n",
      "True\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "import dask.array as da\n",
    "import time\n",
    "\n",
    "# Create a large NumPy array\n",
    "numpy_array = np.random.rand(10**8)\n",
    "\n",
    "# Create a large Dask array\n",
    "# here we split the numpy array into chunks of 10**6\n",
    "dask_array = da.from_array(numpy_array, chunks=10**6).persist()\n",
    "\n",
    "# Perform a computationally intensive operation using NumPy\n",
    "def complex_operation(arr):\n",
    "    return np.sin(arr) + np.cos(arr) * np.tan(arr)\n",
    "\n",
    "start_time = time.time()\n",
    "result_numpy = complex_operation(numpy_array)\n",
    "print(\"NumPy time:\", time.time() - start_time)\n",
    "\n",
    "# Perform the same operation using Dask with parallelism\n",
    "# Dask will automatically manage the parallization on individual chunks for you\n",
    "# its syntax is almost the same as numpy\n",
    "def complex_operation_dask(arr):\n",
    "    return da.sin(arr) + da.cos(arr) * da.tan(arr)\n",
    "\n",
    "start_time = time.time()\n",
    "result_dask = complex_operation_dask(dask_array)\n",
    "result_dask = result_dask.compute()  # Compute the result in parallel. If you don't call compute(), dask will not run and nothing will be returned. This is called lazy evaluation. Remember `range(n)` in Python? \n",
    "print(\"Dask time:\", time.time() - start_time)\n",
    "\n",
    "# Compare the results\n",
    "print(np.allclose(result_numpy, result_dask))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's choose a different chunk size. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/lukeum/anaconda3/envs/py10/lib/python3.10/site-packages/distributed/client.py:3163: UserWarning: Sending large graph of size 762.96 MiB.\n",
      "This may cause some slowdown.\n",
      "Consider scattering data ahead of time and using futures.\n",
      "  warnings.warn(\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Dask time: 37.812100410461426\n",
      "False\n"
     ]
    }
   ],
   "source": [
    "# Create a large NumPy array\n",
    "numpy_array = np.random.rand(10**8)\n",
    "\n",
    "# Create a large Dask array\n",
    "dask_array = da.from_array(numpy_array, chunks=10**5).persist()\n",
    "\n",
    "# Perform the same operation using Dask with parallelism\n",
    "# Dask will automatically manage the parallization on individual chunks for you\n",
    "# its syntax is almost the same as numpy\n",
    "def complex_operation_dask(arr):\n",
    "    return da.sin(arr) + da.cos(arr) * da.tan(arr)\n",
    "\n",
    "start_time = time.time()\n",
    "result_dask = complex_operation_dask(dask_array)\n",
    "result_dask = result_dask.compute()  # Compute the result in parallel. If you don't call compute(), dask will not run and nothing will be returned. This is called lazy evaluation. Remember `range(n)` in Python? \n",
    "print(\"Dask time:\", time.time() - start_time)\n",
    "\n",
    "# Compare the results\n",
    "print(np.allclose(result_numpy, result_dask))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It's now taking a much longer time! \n",
    "Why? \n",
    "\n",
    "There is communication overhead associated with parallelized computation. If the chunksize is small, the processing time might be overwhelmed by the time for communication between different processes, bringing no speed-up at all!\n",
    "\n",
    "Generally speaking, the chunks should be more than the number of cores you allocated (why? if you have 4 cores and split the matrix into 2 chunks, what will happen?).\n",
    "The chunks should not be too small, such that the time needed to process each chunk is much more than the communication overhead. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Paralleling matrices across CPUs can speed up the computation. We can get even much more speed-up by using GPUs!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Massively parallelized matrix operations"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "PyTorch is a machine learning framework based on the Torch library, used for applications such as computer vision and natural language processing, originally developed by Meta AI and now part of the Linux Foundation umbrella. It is free and open-source software released under the modified BSD license."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Numpy arrays are mainly used in typical machine learning algorithms (such as k-means or Decision Tree in scikit-learn) whereas pytorch tensors are mainly used in deep learning which requires heavy matrix computation."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Example:\n",
    "\n",
    "Cosine similarity measures the cosine of the angle between two non-zero vectors $a$ and $b$. It quantifies how similar the two vectors are, regardless of their magnitude (length). Because we normalize them to be length 1.\n",
    "\n",
    "$$\n",
    "sim(a,b) = \\frac{a}{||a||}\\cdot\\frac{b}{||b||} = \\frac{ab}{||a||\\cdot||b||}\n",
    "$$\n",
    "\n",
    "The value of cosine similarity ranges from -1 to 1:\n",
    " - 1 means the vectors are identical (pointing in the same direction).\n",
    " - 0 means the vectors are orthogonal (perpendicular), indicating no similarity.\n",
    " - -1 means the vectors are diametrically opposed (pointing in opposite directions).\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "NumPy Runtime: 0.0028 seconds\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "import torch\n",
    "import time\n",
    "\n",
    "# Define the matrix sizes\n",
    "size = 500000\n",
    "num_iterations = 10\n",
    "\n",
    "# NumPy Matrix Multiplication\n",
    "a = np.random.rand(size)\n",
    "b = np.random.rand(size)\n",
    "\n",
    "start_time = time.time()\n",
    "for _ in range(num_iterations):\n",
    "    np_result = np.dot(a, b)/(np.linalg.norm(a)*np.linalg.norm(b))\n",
    "np_numpy_time = (time.time() - start_time) \n",
    "print(f\"NumPy Runtime: {np_numpy_time:.4f} seconds\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "PyTorch Runtime: 0.0874 seconds\n"
     ]
    }
   ],
   "source": [
    "\n",
    "# PyTorch Matrix Multiplication\n",
    "a = torch.rand(size)\n",
    "b = torch.rand(size)\n",
    "\n",
    "start_time = time.time()\n",
    "for _ in range(num_iterations):\n",
    "    torch_result = torch.dot(a, b)/(torch.norm(a)*torch.norm(b))\n",
    "torch_pytorch_time = (time.time() - start_time)\n",
    "\n",
    "\n",
    "print(f\"PyTorch Runtime: {torch_pytorch_time:.4f} seconds\")\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Pytorch isn't that fast on CPU. However, it is optimized for GPU computation. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Tips: if you don't have a GPU on your computer, you can open this notebook on [Google Colab](https://colab.research.google.com/) and select T4 as your runtime. Then you can run you jobs on GPU. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "PyTorch on GPU Runtime: 0.0014 seconds\n"
     ]
    }
   ],
   "source": [
    "device = 'cuda'\n",
    "# PyTorch Matrix Multiplication\n",
    "a = torch.rand(size).to(device)\n",
    "b = torch.rand(size).to(device)\n",
    "\n",
    "start_time = time.time()\n",
    "for _ in range(num_iterations):\n",
    "    torch_result = torch.dot(a, b)/(torch.norm(a)*torch.norm(b))\n",
    "torch_pytorch_time = (time.time() - start_time) \n",
    "\n",
    "print(f\"PyTorch on GPU Runtime: {torch_pytorch_time:.4f} seconds\")\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What is a GPU? \n",
    "\n",
    "A GPU, or Graphics Processing Unit, is a special type of computer chip originally designed to handle complex tasks related to graphics and images. However, GPU's ability to handle complex matrix operations makes it possible to train neural networks like [AlexNet](https://en.wikipedia.org/wiki/AlexNet). Hence begins the neural network revolution that you and I are currently experiencing. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "PyTorch's significant speed advantage when using GPUs compared to CPUs is due to several key factors:\n",
    "\n",
    "- **Parallelism**: Modern GPUs are designed with thousands of cores, which allow them to perform many operations in parallel. PyTorch is capable of harnessing this parallelism, making it highly efficient for tasks like deep learning where many matrix multiplications occur simultaneously.\n",
    "\n",
    "- **Optimized Libraries**: PyTorch leverages highly optimized GPU libraries, such as NVIDIA cuBLAS and cuDNN, to accelerate mathematical operations. These libraries are finely tuned for specific GPU architectures and provide substantial speed improvements over CPU-based implementations.\n",
    "\n",
    "- **Data Transfer**: When using PyTorch on a GPU, data transfer between the CPU and GPU is typically minimized. This is achieved by loading data directly onto the GPU and keeping it there as much as possible, avoiding costly data transfers that can occur when using CPUs.\n",
    "\n",
    "\n",
    "It's important to note that not all tasks benefit equally from GPU acceleration. Simple, small-scale operations may not see a substantial speed improvement on GPUs compared to CPUs. However, for computationally intensive tasks such as deep learning, scientific simulations, and large-scale linear algebra operations, GPUs can provide orders of magnitude in performance improvements."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Pytorch has similar operations in Numpy, but with a slightly different syntax. For example, an array in numpy is called tensor in pytorch. Below are some simple examples. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create a tensor\n",
    "x = torch.Tensor(2, 3, 4)\n",
    "print(x)\n",
    "\n",
    "# Create a tensor from a (nested) list\n",
    "x = torch.Tensor([[1, 2], [3, 4]])\n",
    "print(x)\n",
    "\n",
    "# Create a tensor with random values between 0 and 1 with the shape [2, 3, 4]\n",
    "x = torch.rand(2, 3, 4)\n",
    "print(x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "shape = x.shape\n",
    "print(\"Shape:\", x.shape)\n",
    "\n",
    "size = x.size()\n",
    "print(\"Size:\", size)\n",
    "\n",
    "dim1, dim2, dim3 = x.size()\n",
    "print(\"Size:\", dim1, dim2, dim3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = torch.arange(12).view(3, 4)\n",
    "print(\"X\", x)\n",
    "\n",
    "print(x[:, 1])   # Second column\n",
    "print(x[0])      # First row\n",
    "print(x[:2, -1]) # First two rows, last column\n",
    "print(x[1:3, :]) # Middle two rows"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = torch.arange(6)\n",
    "x = x.view(2, 3) # view is reshape in numpy\n",
    "print(\"X\", x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "W = torch.arange(9).view(3, 3) # We can also stack multiple operations in a single line. \n",
    "print(\"W\", W)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "h = torch.matmul(x, W) # Verify the result by calculating it by hand too!\n",
    "print(\"h\", h)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Speeding up matrix multiplication with half-precision\n",
    "\n",
    "Note. A GPU is needed to run this experiment.  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In computer engineering, decimal numbers like 1.0151 or 566132.8 are traditionally represented as floating point numbers. Since we can have infinitely precise numbers (think π), but limited space in which to store them, we have to make a compromise between precision (the number of decimals we can include in a number before we have to start rounding it) and size (how many bits we use to store the number).\n",
    "\n",
    "The technical standard for floating point numbers, IEEE 754 (for a deep dive I recommend the PyCon 2019 talk Floats are Friends: making the most of IEEE754.00000000000000002), sets the following standards:\n",
    "\n",
    " - `fp64`, aka double-precision or “double”, max rounding error of ~2^-52\n",
    "\n",
    " - `fp32`, aka single-precision or “single”, max rounding error of ~2^-23\n",
    "\n",
    " - `fp16`, aka half-precision or “half”, max rounding error of ~2^-10\n",
    "\n",
    "Python uses fp64 for the float type. PyTorch, which is much more memory-sensitive, uses fp32 as its default dtype instead.\n",
    "\n",
    "If we halve the precision (fp32 → fp16), we halve the time and space at the cost of precision (only if this is supported by specific hardware, i.e. GPU).\n",
    "\n",
    "It's rare that we need to do highly accurate matrix computations at scale. In fact, often we're doing some kind of machine learning, and less accurate approaches can prevent overfitting.\n",
    "\n",
    "If we accept some decrease in accuracy, then we can often increase speed by orders of magnitude (and/or decrease memory use)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's look at an example with cosine similarity. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "float64\n",
      "NumPy Runtime: 2.5368 seconds\n"
     ]
    }
   ],
   "source": [
    "# Define the matrix sizes\n",
    "size = 500000\n",
    "num_iterations = 10000\n",
    "\n",
    "# NumPy Matrix Multiplication\n",
    "a = np.random.rand(size)\n",
    "b = np.random.rand(size)\n",
    "print(a.dtype)\n",
    "\n",
    "start_time = time.time()\n",
    "for _ in range(num_iterations):\n",
    "    np_result = np.dot(a, b)/(np.linalg.norm(a)*np.linalg.norm(b)) # this line is vectorized\n",
    "np_numpy_time = (time.time() - start_time) \n",
    "print(f\"NumPy Runtime: {np_numpy_time:.4f} seconds\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "float32\n",
      "NumPy Runtime: 1.5467 seconds\n"
     ]
    }
   ],
   "source": [
    "# NumPy Matrix Multiplication\n",
    "a = np.random.rand(size).astype(np.float32)\n",
    "b = np.random.rand(size).astype(np.float32)\n",
    "\n",
    "print(a.dtype)\n",
    "\n",
    "start_time = time.time()\n",
    "for _ in range(num_iterations):\n",
    "    np_result = np.dot(a, b)/(np.linalg.norm(a)*np.linalg.norm(b))\n",
    "np_numpy_time = (time.time() - start_time) \n",
    "print(f\"NumPy Runtime: {np_numpy_time:.4f} seconds\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# NumPy Matrix Multiplication\n",
    "a = np.random.rand(size).astype(np.float16)\n",
    "b = np.random.rand(size).astype(np.float16)\n",
    "\n",
    "start_time = time.time()\n",
    "for _ in range(num_iterations):\n",
    "    np_result = np.dot(a, b)/(np.linalg.norm(a)*np.linalg.norm(b))\n",
    "np_numpy_time = (time.time() - start_time) \n",
    "print(f\"NumPy Runtime: {np_numpy_time:.4f} seconds\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This example above probably gives you a numerical error. Because low precision computation can easily lead to underflow or overflow, triggering errors. \n",
    "\n",
    " - Overflow: When numbers get too big for the computer to handle, they wrap around or behave unexpectedly.\n",
    " - Underflow: When numbers get too small, they are rounded down to zero, losing important information."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's look at how  `pytorch` does it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "torch.float32\n",
      "PyTorch Runtime: 0.6180 seconds\n"
     ]
    }
   ],
   "source": [
    "# Set the size of the vectors to be created\n",
    "size = 500000  \n",
    "\n",
    "# Set the number of iterations \n",
    "num_iterations = 10000  \n",
    "\n",
    "# Specify the device for computation (GPU). You need this to enable GPU computation.\n",
    "device = 'cuda' \n",
    "\n",
    "# Create two random vectors (arrays) using PyTorch\n",
    "a = torch.rand(size).to(device)  \n",
    "b = torch.rand(size).to(device)  \n",
    "\n",
    "\n",
    "print(a.dtype)  \n",
    "\n",
    "# Start measuring the time for the dot product computation\n",
    "start_time = time.time()  \n",
    "\n",
    "\n",
    "for _ in range(num_iterations):  \n",
    "    # Calculate the cosine similarity between vectors 'a' and 'b'\n",
    "    torch_result = torch.dot(a, b) / (torch.norm(a) * torch.norm(b))  \n",
    "\n",
    "\n",
    "torch_pytorch_time = (time.time() - start_time)\n",
    "\n",
    "\n",
    "print(f\"PyTorch Runtime: {torch_pytorch_time:.4f} seconds\")  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "torch.float16\n",
      "PyTorch Runtime: 0.6301 seconds\n"
     ]
    }
   ],
   "source": [
    "\n",
    "# PyTorch Matrix Multiplication\n",
    "a = torch.rand(size).half().to(device)\n",
    "b = torch.rand(size).half().to(device)\n",
    "print(a.dtype)\n",
    "\n",
    "start_time = time.time()\n",
    "for _ in range(num_iterations):\n",
    "    torch_result = torch.dot(a, b)/(torch.norm(a)*torch.norm(b))\n",
    "torch_pytorch_time = (time.time() - start_time)\n",
    "\n",
    "\n",
    "print(f\"PyTorch Runtime: {torch_pytorch_time:.4f} seconds\")\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Using `float16`, aslo known as half precision, can reduce the computation time/memory by half! You might not see the speed up in this example. We need a problem of bigger scale to see the actual improvement. If you increase the size or the iterations by 100 times, you should see the speed up.\n",
    "Of course, a price we must pay is the reduced precision, which makes numerical underflow and overflow more frequent."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Intro to sparse matrices"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "import scipy.sparse\n",
    "import networkx as nx"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Sparse matrices are a class of matrices with most of its elements being 0. \n",
    "Let's initialize a sparse matrix.  \n",
    "(Here we are initializing the graph adjcency matrix, which we will cover in week 3 on graph data structure. Here you just need to know that it is a matrix with a lot of zeros in it.  )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "G = nx.ladder_graph(5)\n",
    "# Return the Ladder graph of length n.\n",
    "# This is two rows of n nodes, with each pair connected by a single edge.\n",
    "am_ladder = nx.adjacency_matrix(G)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "scipy.sparse._csr.csr_matrix"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "type(am_ladder) # Compressed Sparse Row"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- Sparse matrices are a conceptual data structure like a list, dictionary, set, etc.\n",
    "- [`scipy.sparse`](https://docs.scipy.org/doc/scipy/reference/sparse.html) matrices are the standard Python implementation of this conceptual data structure, like `list`, `dict`, `set`, etc.\n",
    "- Going to that link, we can see there are many types of scipy sparse matrix.\n",
    "  - This one is a [`csr_matrix`](https://docs.scipy.org/doc/scipy/reference/generated/scipy.sparse.csr_matrix.html#scipy.sparse.csr_matrix)\n",
    "  - More later on these types.\n",
    "- You can convert them to numpy arrays with `toarray()`, but this is often a bad idea."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[0 1 0 0 0 1 0 0 0 0]\n",
      " [1 0 1 0 0 0 1 0 0 0]\n",
      " [0 1 0 1 0 0 0 1 0 0]\n",
      " [0 0 1 0 1 0 0 0 1 0]\n",
      " [0 0 0 1 0 0 0 0 0 1]\n",
      " [1 0 0 0 0 0 1 0 0 0]\n",
      " [0 1 0 0 0 1 0 1 0 0]\n",
      " [0 0 1 0 0 0 1 0 1 0]\n",
      " [0 0 0 1 0 0 0 1 0 1]\n",
      " [0 0 0 0 1 0 0 0 1 0]]\n"
     ]
    }
   ],
   "source": [
    "print(am_ladder.toarray())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[0 1 0 0 0 1 0 0 0 0]\n",
      " [1 0 1 0 0 0 1 0 0 0]\n",
      " [0 1 0 1 0 0 0 1 0 0]\n",
      " [0 0 1 0 1 0 0 0 1 0]\n",
      " [0 0 0 1 0 0 0 0 0 1]\n",
      " [1 0 0 0 0 0 1 0 0 0]\n",
      " [0 1 0 0 0 1 0 1 0 0]\n",
      " [0 0 1 0 0 0 1 0 1 0]\n",
      " [0 0 0 1 0 0 0 1 0 1]\n",
      " [0 0 0 0 1 0 0 0 1 0]]\n",
      "---\n",
      "[[0 1 0 0 0 1 0 0 0 0]\n",
      " [1 0 1 0 0 0 1 0 0 0]\n",
      " [0 1 0 1 0 0 0 1 0 0]\n",
      " [0 0 1 0 1 0 0 0 1 0]\n",
      " [0 0 0 1 0 0 0 0 0 1]\n",
      " [1 0 0 0 0 0 1 0 0 0]\n",
      " [0 1 0 0 0 1 0 1 0 0]\n",
      " [0 0 1 0 0 0 1 0 1 0]\n",
      " [0 0 0 1 0 0 0 1 0 1]\n",
      " [0 0 0 0 1 0 0 0 1 0]]\n"
     ]
    }
   ],
   "source": [
    "A = am_ladder.toarray() # undirect graph = symmetric matrix;\n",
    "\n",
    "print(A)\n",
    "print(\"---\")\n",
    "print(A.T)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "scipy.sparse._csr.csr_array"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# make a bigger sparse matrix\n",
    "G = nx.fast_gnp_random_graph(100_000,1e-4) \n",
    "am = nx.adjacency_matrix(G)\n",
    "type(am)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(100000, 100000)"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "am.shape # 10 billion"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "999120"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "am.nnz # non-zero elements of the matrix "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Stored in full form, the matrix would take up:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "3.2 ms ± 1.23 ms per loop (mean ± std. dev. of 7 runs, 100 loops each)\n"
     ]
    }
   ],
   "source": [
    "%timeit np.sum(am, axis = 0) # row-wise sum"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "8\n",
      "[1 7]\n",
      "[2 6]\n",
      "[1 7]\n",
      "[2 6]\n"
     ]
    }
   ],
   "source": [
    "# https://numpy.org/doc/stable/reference/generated/numpy.sum.html\n",
    "print(np.sum([[0, 1], [2, 5]], axis=None))  # default\n",
    "\n",
    "print(np.sum([[0, 1], [2, 5]], axis=1))     # row-wise sum\n",
    "print(np.sum([[0, 1], [2, 5]], axis=0))     # column-wise sum\n",
    "\n",
    "print(np.sum([[0, 1], [2, 5]], axis=-1))     # last axis ->  row\n",
    "print(np.sum([[0, 1], [2, 5]], axis=-2))     # second last axis ->  column"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The full matrix would take up 80 GB\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "\n",
    "full_size = int(np.prod(am.shape))*8/(1e9)\n",
    "print(\"The full matrix would take up %d GB\" % full_size)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Size of the array:  10000000000\n",
      "Memory size of one array element in bytes:  8\n",
      "Memory size of numpy array in bytes: 80000000000 Bytes\n",
      "Memory size of numpy array in gigabytes: 80.0 GB\n"
     ]
    }
   ],
   "source": [
    "x = am.toarray()\n",
    "\n",
    "print(\"Size of the array: \",\n",
    "      x.size)\n",
    " \n",
    "print(\"Memory size of one array element in bytes: \",\n",
    "      x.itemsize)\n",
    "\n",
    "print(\"Memory size of numpy array in bytes:\",\n",
    "      x.size * x.itemsize, \"Bytes\")\n",
    "\n",
    "print(\"Memory size of numpy array in gigabytes:\",\n",
    "      (x.size * x.itemsize) / (1e9), \"GB\")                        # 1 byte = 1e-9 gigabyte"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "That's a lot! How big is the sparse matrix?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The sparse matrix takes up 7 NB\n"
     ]
    }
   ],
   "source": [
    "sparse_size = am.data.nbytes/1e6\n",
    "print(\"The sparse matrix takes up %d NB\" % sparse_size)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So, the fraction of space saved is:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The sparse matrix is 10009x smaller\n"
     ]
    }
   ],
   "source": [
    "frac_nz = am.nnz / np.prod(am.shape)\n",
    "print(\"The sparse matrix is %dx smaller\" % (1/frac_nz))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- Right, so we definitely don't want to store the full matrix."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- Regular numpy functions work on sparse matrices, although they might be fast/slow depending on various factors. \n",
    "- You definitely do not want to iterate through the rows - make sure you use built-in numpy functions. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Sparse datasets\n",
    "\n",
    "Sparse matrices come up _a lot_ in practice. For example:\n",
    "\n",
    "- Word counts: we might represent a document by which words are in it, but only a small fraction of all words would appear in a given document.\n",
    "- Ratings: we might represent an Amazon item by the user ratings, but only a small fraction of all users have rated a given item.\n",
    "- Physical processes: in a 2019 Capstone project, MDS students examined images from a particle physics dataset, in which most of the sensors got zero signal.\n",
    "- etc."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Building sparse matrices"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To create an empty sparse martrix, just provide the shape as the argument to the constructor"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(10, 10)"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "from scipy.sparse import csr_matrix,csc_matrix,lil_matrix\n",
    "# csr: Compressed Sparse Row array\n",
    "# csc: Compressed Sparse Column array\n",
    "# lil: Row-based list of lists sparse array\n",
    "\n",
    "shape = (10,10)\n",
    "#my code here\n",
    "x = csr_matrix(shape)\n",
    "#my code here\n",
    "x.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    " However, this is a bad way to build a sparse matrix of any significant size. In fact, scipy will warn you that it's a very bad idea to try to populate a CSR sparse matrix directly"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/jungyeul/Library/Python/3.9/lib/python/site-packages/scipy/sparse/_index.py:103: SparseEfficiencyWarning: Changing the sparsity structure of a csr_matrix is expensive. lil_matrix is more efficient.\n",
      "  self._set_intXint(row, col, x.flat[0])\n"
     ]
    }
   ],
   "source": [
    "x[5,5] = 1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A better way? Build the matrix in another sparse format, and then convert. Let build a diagonal matrix using the lil (linked list, i.e. the treasure boxes) format, which doesn't support efficient row or column operations like the csr and csc matrices, but does support efficient assignment."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  (0, 0)\t1.0\n",
      "  (1, 1)\t1.0\n",
      "  (2, 2)\t1.0\n",
      "  (3, 3)\t1.0\n",
      "  (4, 4)\t1.0\n",
      "  (5, 5)\t1.0\n",
      "  (6, 6)\t1.0\n",
      "  (7, 7)\t1.0\n",
      "  (8, 8)\t1.0\n",
      "  (9, 9)\t1.0\n"
     ]
    }
   ],
   "source": [
    "#my code here\n",
    "shape = (10,10)\n",
    "\n",
    "x = lil_matrix(shape)\n",
    "for i in range(shape[0]):\n",
    "    x[i,i] = 1\n",
    "    \n",
    "x = csr_matrix(x)\n",
    "print(x)\n",
    "#my code here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The Scikit learn package allows you to create csr sparse matrices from Python dictionaries (using [DictVectorizer](https://scikit-learn.org/stable/modules/generated/sklearn.feature_extraction.DictVectorizer.html)) and even directly from raw text strings (using [CountVectorizer](https://scikit-learn.org/stable/modules/generated/sklearn.feature_extraction.text.CountVectorizer.html);).\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "['bar' 'baz' 'foo']\n",
      "[[2. 0. 1.]\n",
      " [0. 1. 3.]]\n",
      "<class 'numpy.ndarray'>\n",
      "----\n"
     ]
    }
   ],
   "source": [
    "from sklearn.feature_extraction import DictVectorizer\n",
    "\n",
    "dv = DictVectorizer(sparse=False)\n",
    "D1 = [{'foo': 1, 'bar': 2}, {'foo': 3, 'baz': 1}]\n",
    "X1 = dv.fit_transform(D1)\n",
    "print(dv.get_feature_names_out())\n",
    "print(X1)\n",
    "print(type(X1))\n",
    "# print(dv.vocabulary_)\n",
    "print(\"----\")\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "['bar' 'baz' 'foo']\n",
      "  (0, 2)\t1\n",
      "  (0, 0)\t1\n",
      "  (0, 1)\t1\n",
      "  (1, 2)\t1\n",
      "  (1, 0)\t1\n",
      "  (2, 2)\t1\n",
      "  (3, 0)\t1\n",
      "  (4, 0)\t1\n",
      "<class 'scipy.sparse._csr.csr_matrix'>\n",
      "[[1 1 1]\n",
      " [1 0 1]\n",
      " [0 0 1]\n",
      " [1 0 0]\n",
      " [1 0 0]]\n"
     ]
    }
   ],
   "source": [
    "from sklearn.feature_extraction.text import CountVectorizer\n",
    "\n",
    "cv = CountVectorizer()\n",
    "D2 = ['foo bar baz', 'bar foo', 'foo', 'bar', 'bar']\n",
    "X2 = cv.fit_transform(D2)\n",
    "print(cv.get_feature_names_out())\n",
    "print(X2)\n",
    "print(type(X2))\n",
    "print(X2.toarray())\n",
    "# print(cv.vocabulary_)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Optional example"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "[nltk_data] Downloading package brown to /home/lukeum/nltk_data...\n",
      "[nltk_data]   Unzipping corpora/brown.zip.\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "True"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# #pip install nltk\n",
    "import nltk\n",
    "nltk.download(\"brown\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "from nltk.corpus import brown"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Then, we apply the vectorizer to convert this to a csr matrix. First, initialize the vectorizer, then use its `fit_transform` method"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`CountVectorizer()`\n",
    "\n",
    "Convert a collection of text documents to a matrix of token counts.\n",
    "\n",
    "This implementation produces a sparse representation of the counts using scipy.sparse.csr_matrix."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "284\n"
     ]
    }
   ],
   "source": [
    "from sklearn.feature_extraction.text import CountVectorizer\n",
    "\n",
    "count_vec = CountVectorizer()\n",
    "count_vec_output = count_vec.fit_transform(brown.words('ca01')) \n",
    "\n",
    "# print the identified unique words along with their indices\n",
    "print(count_vec.vocabulary_['fulton'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "<class 'scipy.sparse._csr.csr_matrix'>\n",
      "[[0 0 0 ... 0 0 0]\n",
      " [0 0 0 ... 0 0 0]\n",
      " [0 0 0 ... 0 0 0]\n",
      " ...\n",
      " [0 0 0 ... 0 0 0]\n",
      " [0 0 0 ... 0 0 0]\n",
      " [0 0 0 ... 0 0 0]]\n"
     ]
    }
   ],
   "source": [
    "#print(count_vec_output)\n",
    "print(type(count_vec_output))\n",
    "print(count_vec_output.toarray())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Pitfalls "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Matrix vs. Array "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- Note that when you do an operation which on sparse matrix that removes one of the dimensions (as we did above), you actually don't have a sparse matrix anymore. \n",
    "- Surprisingly, you also don't have a numpy array! \n",
    "- What you have a numpy [matrix](https://docs.scipy.org/doc/numpy/reference/generated/numpy.matrix.html), which is distinct from an array. \n",
    "    - See [here](https://numpy.org/devdocs/user/numpy-for-matlab-users.html#array-or-matrix-which-should-i-use) for a discussion of the differences between arrays and matrices. \n",
    "- Rule of thumb: (non-sparse) matrices will cause you headaches and you should stay away from them. \n",
    "- You can access an array version of the matrix (without copying it) using `np.asarray`, and you should probably do that right away"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "numpy.matrix"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "type(wpt_matrix)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[2242.],\n",
       "       [2277.],\n",
       "       [2275.],\n",
       "       [2217.],\n",
       "       [2244.],\n",
       "       [2263.],\n",
       "       [2270.],\n",
       "       [2187.],\n",
       "       [2234.],\n",
       "       [2282.]])"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "wpt_array = np.asarray(wpt_matrix)\n",
    "wpt_array[:10]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "An example of where a numpy matrix will cause you trouble: operations like `flatten` (which is suppose to reduce a dataset to 1 dimension) won't work properly with a matrix, which is *always* 2-d"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "matrix([[2242., 2277., 2275., 2217., 2244., 2263., 2270., 2187., 2234.,\n",
       "         2282., 2259., 2338., 2241., 2329., 2314., 2374., 2293., 2314.,\n",
       "         2260., 2317., 2245., 2309., 2327., 2252., 2327., 2278., 2288.,\n",
       "         2220., 2349., 2274., 2311., 2323., 2284., 2235., 2217., 2229.,\n",
       "         2324., 2221., 2452., 2340., 2313., 2338., 2291., 2277., 2200.,\n",
       "         2234., 2244., 2230., 2241., 2231., 2239., 2395., 2312., 2298.,\n",
       "         2293., 2296., 2267., 2264., 2311., 2260., 2288., 2194., 2287.,\n",
       "         2299., 2278., 2294., 2357., 2346., 2368., 2316., 2262., 2415.,\n",
       "         2320., 2356., 2342., 2398., 2433., 2421., 2351., 2345., 2461.,\n",
       "         2480., 2342., 2282., 2574., 2503., 2311., 2370., 2213., 2334.,\n",
       "         2332., 2318., 2410., 2317., 2359., 2332., 2314., 2238., 2304.,\n",
       "         2307., 2236., 2276., 2315., 2451., 2343., 2258., 2279., 2257.,\n",
       "         2295., 2271., 2222., 2370., 2226., 2476., 2402., 2333., 2288.,\n",
       "         2261., 2346., 2333., 2223., 2244., 2302., 2260., 2251., 2239.,\n",
       "         2402., 2307., 2253., 2227., 2243., 2323., 2291., 2246., 2294.,\n",
       "         2247., 2339., 2233., 2300., 2272., 2232., 2369., 2264., 2256.,\n",
       "         2258., 2250., 2424., 2205., 2356., 2306., 2332., 2310., 2361.,\n",
       "         2316., 2285., 2236., 2331., 2236., 2328., 2321., 2278., 2212.,\n",
       "         2431., 2299., 2384., 2313., 2258., 2288., 2267., 2302., 2259.,\n",
       "         2220., 2322., 2219., 2300., 2439., 2362., 2222., 2405., 2351.,\n",
       "         2250., 2261., 2355., 2288., 2245., 2261., 2267., 2275., 2222.,\n",
       "         2289., 2311., 2403., 2359., 2394., 2397., 2308., 2322., 2278.,\n",
       "         2301., 2237., 2332., 2311., 2320., 2514., 2226., 2471., 2289.,\n",
       "         2259., 2286., 2245., 2291., 2316., 2211., 2255., 2195., 2359.,\n",
       "         2211., 2240., 2328., 2299., 2314., 2335., 2351., 2285., 2386.,\n",
       "         2275., 2317., 2436., 2318., 2276., 2221., 2293., 2238., 2357.,\n",
       "         2214., 2433., 2204., 2415., 2312., 2278., 2277., 2307., 2328.,\n",
       "         2266., 2282., 2349., 2286., 2312., 2251., 2267., 2285., 2209.,\n",
       "         2347., 2235., 2302., 2389., 2348., 2407., 2262., 2295., 2247.,\n",
       "         2334., 2408., 2288., 2416., 2360., 2207., 2316., 2228., 2411.,\n",
       "         2436., 2318., 2451., 2313., 2254., 2412., 2408., 2422., 2313.,\n",
       "         2346., 2539., 2250., 2194., 2406., 2296., 2559., 2282., 2471.,\n",
       "         2305., 2200., 2248., 2317., 2225., 2214., 2189., 2182., 2205.,\n",
       "         2334., 2276., 2148., 2370., 2206., 2302., 2262., 2277., 2328.,\n",
       "         2294., 2298., 2317., 2258., 2291., 2290., 2333., 2211., 2229.,\n",
       "         2185., 2337., 2223., 2292., 2299., 2257., 2243., 2230., 2267.,\n",
       "         2335., 2170., 2241., 2208., 2401., 2262., 2186., 2174., 2339.,\n",
       "         2171., 2250., 2304., 2250., 2238., 2214., 2266., 2260., 2301.,\n",
       "         2189., 2319., 2257., 2255., 2327., 2255., 2222., 2264., 2239.,\n",
       "         2267., 2270., 2276., 2331., 2206., 2299., 2291., 2386., 2313.,\n",
       "         2210., 2302., 2281., 2187., 2296., 2372., 2486., 2306., 2242.,\n",
       "         2341., 2525., 2371., 2246., 2254., 2440., 2325., 2296., 2310.,\n",
       "         2345., 2527., 2443., 2375., 2281., 2396., 2399., 2425., 2364.,\n",
       "         2340., 2330., 2319., 2380., 2218., 2345., 2378., 2376., 2404.,\n",
       "         2319., 2419., 2230., 2404., 2319., 2517., 2264., 2374., 2370.,\n",
       "         2428., 2285., 2371., 2288., 2331., 2489., 2393., 2391., 2341.,\n",
       "         2437., 2405., 2485., 2358., 2419., 2337., 2324., 2330., 2455.,\n",
       "         2414., 2358., 2340., 2446., 2486., 2412., 2364., 2461., 2369.,\n",
       "         2378., 2448., 2383., 2380., 2349., 2434., 2279., 2403., 2289.,\n",
       "         2483., 2471., 2308., 2424., 2395., 2381., 2318., 2405., 2529.,\n",
       "         2275., 2457., 2356., 2341., 2313., 2349., 2523., 2418., 2341.,\n",
       "         2449., 2281., 2560., 2332., 2320., 2422., 2546., 2346., 2538.,\n",
       "         2385., 2317., 2187., 2386., 2346., 2463., 2338., 2501., 2524.,\n",
       "         2525., 2288., 2469., 2449., 2428., 2520., 2441., 2546., 2524.,\n",
       "         2310., 2356., 2459., 2339., 2417., 2331., 2392., 2587., 2256.,\n",
       "         2464., 2506., 2456., 2371., 2332.]])"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "wpt_matrix.flatten()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "But this will work fine with the array version"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([2242., 2277., 2275., 2217., 2244., 2263., 2270., 2187., 2234.,\n",
       "       2282., 2259., 2338., 2241., 2329., 2314., 2374., 2293., 2314.,\n",
       "       2260., 2317., 2245., 2309., 2327., 2252., 2327., 2278., 2288.,\n",
       "       2220., 2349., 2274., 2311., 2323., 2284., 2235., 2217., 2229.,\n",
       "       2324., 2221., 2452., 2340., 2313., 2338., 2291., 2277., 2200.,\n",
       "       2234., 2244., 2230., 2241., 2231., 2239., 2395., 2312., 2298.,\n",
       "       2293., 2296., 2267., 2264., 2311., 2260., 2288., 2194., 2287.,\n",
       "       2299., 2278., 2294., 2357., 2346., 2368., 2316., 2262., 2415.,\n",
       "       2320., 2356., 2342., 2398., 2433., 2421., 2351., 2345., 2461.,\n",
       "       2480., 2342., 2282., 2574., 2503., 2311., 2370., 2213., 2334.,\n",
       "       2332., 2318., 2410., 2317., 2359., 2332., 2314., 2238., 2304.,\n",
       "       2307., 2236., 2276., 2315., 2451., 2343., 2258., 2279., 2257.,\n",
       "       2295., 2271., 2222., 2370., 2226., 2476., 2402., 2333., 2288.,\n",
       "       2261., 2346., 2333., 2223., 2244., 2302., 2260., 2251., 2239.,\n",
       "       2402., 2307., 2253., 2227., 2243., 2323., 2291., 2246., 2294.,\n",
       "       2247., 2339., 2233., 2300., 2272., 2232., 2369., 2264., 2256.,\n",
       "       2258., 2250., 2424., 2205., 2356., 2306., 2332., 2310., 2361.,\n",
       "       2316., 2285., 2236., 2331., 2236., 2328., 2321., 2278., 2212.,\n",
       "       2431., 2299., 2384., 2313., 2258., 2288., 2267., 2302., 2259.,\n",
       "       2220., 2322., 2219., 2300., 2439., 2362., 2222., 2405., 2351.,\n",
       "       2250., 2261., 2355., 2288., 2245., 2261., 2267., 2275., 2222.,\n",
       "       2289., 2311., 2403., 2359., 2394., 2397., 2308., 2322., 2278.,\n",
       "       2301., 2237., 2332., 2311., 2320., 2514., 2226., 2471., 2289.,\n",
       "       2259., 2286., 2245., 2291., 2316., 2211., 2255., 2195., 2359.,\n",
       "       2211., 2240., 2328., 2299., 2314., 2335., 2351., 2285., 2386.,\n",
       "       2275., 2317., 2436., 2318., 2276., 2221., 2293., 2238., 2357.,\n",
       "       2214., 2433., 2204., 2415., 2312., 2278., 2277., 2307., 2328.,\n",
       "       2266., 2282., 2349., 2286., 2312., 2251., 2267., 2285., 2209.,\n",
       "       2347., 2235., 2302., 2389., 2348., 2407., 2262., 2295., 2247.,\n",
       "       2334., 2408., 2288., 2416., 2360., 2207., 2316., 2228., 2411.,\n",
       "       2436., 2318., 2451., 2313., 2254., 2412., 2408., 2422., 2313.,\n",
       "       2346., 2539., 2250., 2194., 2406., 2296., 2559., 2282., 2471.,\n",
       "       2305., 2200., 2248., 2317., 2225., 2214., 2189., 2182., 2205.,\n",
       "       2334., 2276., 2148., 2370., 2206., 2302., 2262., 2277., 2328.,\n",
       "       2294., 2298., 2317., 2258., 2291., 2290., 2333., 2211., 2229.,\n",
       "       2185., 2337., 2223., 2292., 2299., 2257., 2243., 2230., 2267.,\n",
       "       2335., 2170., 2241., 2208., 2401., 2262., 2186., 2174., 2339.,\n",
       "       2171., 2250., 2304., 2250., 2238., 2214., 2266., 2260., 2301.,\n",
       "       2189., 2319., 2257., 2255., 2327., 2255., 2222., 2264., 2239.,\n",
       "       2267., 2270., 2276., 2331., 2206., 2299., 2291., 2386., 2313.,\n",
       "       2210., 2302., 2281., 2187., 2296., 2372., 2486., 2306., 2242.,\n",
       "       2341., 2525., 2371., 2246., 2254., 2440., 2325., 2296., 2310.,\n",
       "       2345., 2527., 2443., 2375., 2281., 2396., 2399., 2425., 2364.,\n",
       "       2340., 2330., 2319., 2380., 2218., 2345., 2378., 2376., 2404.,\n",
       "       2319., 2419., 2230., 2404., 2319., 2517., 2264., 2374., 2370.,\n",
       "       2428., 2285., 2371., 2288., 2331., 2489., 2393., 2391., 2341.,\n",
       "       2437., 2405., 2485., 2358., 2419., 2337., 2324., 2330., 2455.,\n",
       "       2414., 2358., 2340., 2446., 2486., 2412., 2364., 2461., 2369.,\n",
       "       2378., 2448., 2383., 2380., 2349., 2434., 2279., 2403., 2289.,\n",
       "       2483., 2471., 2308., 2424., 2395., 2381., 2318., 2405., 2529.,\n",
       "       2275., 2457., 2356., 2341., 2313., 2349., 2523., 2418., 2341.,\n",
       "       2449., 2281., 2560., 2332., 2320., 2422., 2546., 2346., 2538.,\n",
       "       2385., 2317., 2187., 2386., 2346., 2463., 2338., 2501., 2524.,\n",
       "       2525., 2288., 2469., 2449., 2428., 2520., 2441., 2546., 2524.,\n",
       "       2310., 2356., 2459., 2339., 2417., 2331., 2392., 2587., 2256.,\n",
       "       2464., 2506., 2456., 2371., 2332.])"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "wpt_array.flatten()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Indexing"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Another issue with sparse matrices is indexing. With a regular numpy array, `x[i,j]` and `x[i][j]` are equivalent:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = np.random.rand(10,10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.5267092734326778"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "x[1,2]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.5267092734326778"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "x[1][2]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is because `x[1]` returns the first row, and then the `[2]` indexes into that row:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0.92771172, 0.76413423, 0.52670927, 0.24739091, 0.32841451,\n",
       "       0.75422245, 0.1076381 , 0.28458965, 0.94934966, 0.90127702])"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "x[1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.5267092734326778"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "row_1 = x[1]\n",
    "row_1[2]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "However, with `scipy.sparse` matrices, things are a bit different:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x_sparse = csr_matrix(x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.5267092734326778"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "x_sparse[1,2]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "ename": "IndexError",
     "evalue": "row index (2) out of range",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m\n",
      "\u001b[0;31mIndexError\u001b[0m                                Traceback (most recent call last)\n",
      "Cell \u001b[0;32mIn [237], line 1\u001b[0m\n",
      "\u001b[0;32m----> 1\u001b[0m x_sparse[\u001b[39m1\u001b[39;49m][\u001b[39m2\u001b[39;49m]\n",
      "\n",
      "File \u001b[0;32m~/Library/Python/3.9/lib/python/site-packages/scipy/sparse/_index.py:47\u001b[0m, in \u001b[0;36mIndexMixin.__getitem__\u001b[0;34m(self, key)\u001b[0m\n",
      "\u001b[1;32m     46\u001b[0m \u001b[39mdef\u001b[39;00m \u001b[39m__getitem__\u001b[39m(\u001b[39mself\u001b[39m, key):\n",
      "\u001b[0;32m---> 47\u001b[0m     row, col \u001b[39m=\u001b[39m \u001b[39mself\u001b[39;49m\u001b[39m.\u001b[39;49m_validate_indices(key)\n",
      "\u001b[1;32m     49\u001b[0m     \u001b[39m# Dispatch to specialized methods.\u001b[39;00m\n",
      "\u001b[1;32m     50\u001b[0m     \u001b[39mif\u001b[39;00m \u001b[39misinstance\u001b[39m(row, INT_TYPES):\n",
      "\n",
      "File \u001b[0;32m~/Library/Python/3.9/lib/python/site-packages/scipy/sparse/_index.py:155\u001b[0m, in \u001b[0;36mIndexMixin._validate_indices\u001b[0;34m(self, key)\u001b[0m\n",
      "\u001b[1;32m    153\u001b[0m row \u001b[39m=\u001b[39m \u001b[39mint\u001b[39m(row)\n",
      "\u001b[1;32m    154\u001b[0m \u001b[39mif\u001b[39;00m row \u001b[39m<\u001b[39m \u001b[39m-\u001b[39mM \u001b[39mor\u001b[39;00m row \u001b[39m>\u001b[39m\u001b[39m=\u001b[39m M:\n",
      "\u001b[0;32m--> 155\u001b[0m     \u001b[39mraise\u001b[39;00m \u001b[39mIndexError\u001b[39;00m(\u001b[39m'\u001b[39m\u001b[39mrow index (\u001b[39m\u001b[39m%d\u001b[39;00m\u001b[39m) out of range\u001b[39m\u001b[39m'\u001b[39m \u001b[39m%\u001b[39m row)\n",
      "\u001b[1;32m    156\u001b[0m \u001b[39mif\u001b[39;00m row \u001b[39m<\u001b[39m \u001b[39m0\u001b[39m:\n",
      "\u001b[1;32m    157\u001b[0m     row \u001b[39m+\u001b[39m\u001b[39m=\u001b[39m M\n",
      "\n",
      "\u001b[0;31mIndexError\u001b[0m: row index (2) out of range"
     ]
    }
   ],
   "source": [
    "x_sparse[1][2]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Why?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "row_1_sparse = x_sparse[1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(1, 10)"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "row_1_sparse.shape\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- The sparse matrix returns a different shape.\n",
    "- This is because sparse matrices must always be 2d (same problem as non-sparse matrices)\n",
    "- In general, use the `x[1,2]` notation when possible because chaining the `[]` can be problematic in several places (e.g., also pandas).\n",
    "- However, **this is only for numpy**, not, say, a list of lists:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "lst = [[1,2,3],[4,5,6],[7,9]]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "lst[0][1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "ename": "TypeError",
     "evalue": "list indices must be integers or slices, not tuple",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m\n",
      "\u001b[0;31mTypeError\u001b[0m                                 Traceback (most recent call last)\n",
      "Cell \u001b[0;32mIn [150], line 1\u001b[0m\n",
      "\u001b[0;32m----> 1\u001b[0m lst[\u001b[39m0\u001b[39;49m,\u001b[39m1\u001b[39;49m]\n",
      "\n",
      "\u001b[0;31mTypeError\u001b[0m: list indices must be integers or slices, not tuple"
     ]
    }
   ],
   "source": [
    "lst[0,1]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Rows versus columns operations"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Generally, you should avoid looping directly through sparse matrices.\n",
    "\n",
    "But if you do need to, you need to be very sensitive to the difference between Row (CSR) and Column (CSC) sparse matricies. \n",
    "\n",
    "Note. CSR (Compressed Sparse Row) and CSC (Compressed Sparse Column) are two common formats for storing sparse matrices. Here’s a breakdown of the differences:\n",
    "\n",
    "| Feature               | CSR                          | CSC                          |\n",
    "|----------------------|------------------------------|------------------------------|\n",
    "| Storage Arrays       | Values, Column Indices, Row Pointers | Values, Row Indices, Column Pointers |\n",
    "| Optimized For        | Row operations               | Column operations             |\n",
    "| Access Pattern       | Fast row slicing             | Fast column slicing           |\n",
    "| Use Cases            | Matrix-vector multiplication, row-wise operations | Column-wise computations, matrix factorization |\n",
    "\n",
    "You should not use CSR matrices for operations which require looping over columns. Observe:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "num_rows, num_cols = am.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "scipy.sparse._csr.csr_matrix"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "type(am)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2.25 ms ± 39.4 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)\n"
     ]
    }
   ],
   "source": [
    "%%timeit \n",
    "(np.sum(am, axis = 0))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "3.22 ms ± 14.6 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)\n"
     ]
    }
   ],
   "source": [
    "%%timeit \n",
    "(np.sum(am, axis = 1))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "3.92 s ± 0 ns per loop (mean ± std. dev. of 1 run, 1 loop each)\n"
     ]
    }
   ],
   "source": [
    "%%timeit -r1 -n1 \n",
    "\n",
    "# iterate through the rows - this is much slower than numpy functions\n",
    "for i in range(num_rows):\n",
    "    am[i,:]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "6min 17s ± 0 ns per loop (mean ± std. dev. of 1 run, 1 loop each)\n"
     ]
    }
   ],
   "source": [
    "# %%timeit -r1 -n1\n",
    "\n",
    "# # iterate through the columns - this is much MUCH slower\n",
    "# for i in range(num_cols):\n",
    "#     am[:,i]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Above: \n",
    "\n",
    "- Looping is slower, as per usual.\n",
    "- But at least looping through the rows of a `csr_matrix` isn't that bad.\n",
    "- However, looping through the columns of a `csr_matrix` is a disaster - it took several min on my laptop!!\n",
    "  - Because it is stored _row by row_. \n",
    "  - To grab a single column, it needs to loop through each row and look for items in that column.\n",
    "- If you want to loop through columns you should first convert the matrix type."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "am_csc = am.tocsc()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "3.98 s ± 0 ns per loop (mean ± std. dev. of 1 run, 1 loop each)\n"
     ]
    }
   ],
   "source": [
    "%%timeit -r1 -n1\n",
    "\n",
    "for i in range(num_cols):\n",
    "    am_csc[:,i]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Optional: Einstein Summation (einsum)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The Einstein summation convention is the ultimate generalization of products such as matrix multiplication to multiple dimensions. It offers a compact and elegant way of specifying almost any product of scalars/vectors/matrices/tensors. Despite its generality, it can reduce the number of errors made by computer scientists and reduce the time they spend reasoning about linear algebra. It does so by being simultaneously clearer, more explicit, more self-documenting, more declarative in style and less cognitively burdensome to use. Its advantages over such things as matrix multiplication are that it liberates its user from having to think about:\n",
    "\n",
    "- The correct order in which to supply the argument tensors\n",
    "- The correct transpositions to apply to the argument tensors\n",
    "- nsuring that the correct tensor dimensions are lined up with one another\n",
    "- The correct transposition to apply to the resulting tensor\n",
    "\n",
    "\n",
    "The only things it requires are knowledge of:\n",
    "\n",
    "- Along which dimensions to compute (inner/element-wise/outer) products.\n",
    "- The desired output shape.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[EINSUM IS ALL YOU NEED - EINSTEIN SUMMATION IN DEEP LEARNING](https://rockt.github.io/2018/04/30/einsum)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.11.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}