11-LinearAlgebra-Part2.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "e689afa4",
   "metadata": {},
   "source": [
    "# Linear Algebra (Part 2): Matrices and Eigenstuff\n",
    "\n",
    "\n",
    "### After completing this notebook, you'll be able to:\n",
    "* Construct and multiply matrices in Python (and by hand)\n",
    "* Create and manipulate special cases of matrices (unit matrix, diagonal matrix)\n",
    "* Explain matrices as a linear transformation and relate matrix properties to properties of that linear transformation\n",
    "* Define what eigenvalues/eigenvectors are and determine them using Python\n",
    "\n",
    "<hr>\n",
    "\n",
    "## Setup\n",
    "Below, we'll import a custom module with some **helper functions** to easily visualize vectors."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a54e6752",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Import modules, including a custom one!\n",
    "from matrices import *\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e639f00b",
   "metadata": {},
   "source": [
    "## Matrices\n",
    "\n",
    "A **matrix** is a rectangular array of numbers. The numbers in the matrix are **entries**. We can refer to each of the entries by their row and column.\n",
    "```\n",
    "A = [-1 -6  2\n",
    "      0  8  1]\n",
    "```\n",
    "You will sometimes see matrices referred to in mathematical notation. For example, a matrix $A$ could be denoted by $A = (a_{ij})$ where $a_{ij}$  is the entry in the $i^{th}$ row and $j^{th}$ column of matrix $A$. So for the matrix above, $A_{23}$ == 1.\n",
    "\n",
    "The **size** of the matrix is the number of rows multipled by the number of columns. $A$ is a 2x3 matrix, with a total size of 6. If a matrix has the same number of rows and columns, it is a **square** matrix. If it only has a dimension of one in one direction (e.g. 3x1), it is a **column matrix** (or **column vector**). \n",
    "\n",
    "The **transpose** of a matrix switches its rows and columns. The transpose (^T, or $^T$) of matrix $A$ above would be:\n",
    "```\n",
    "A^T = [-1 0\n",
    "       -6 8\n",
    "        2 1]\n",
    "```\n",
    "\n",
    "### Buiding matrices in NumPy\n",
    "We can build matrices in Python using numpy, using the following notation (notice there are parentheses to denote the array function, with brackets inside to indicate a list, with brackets inside *those* brackets for each row):\n",
    "\n",
    "```\n",
    "my_matrix = np.array([[row_1],[row_2],...[row_n]])\n",
    "```\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "dd3b1ae8",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Build a 3x3 matrix\n",
    "my_matrix = np.array([[1,2,3],[4,5,6],[7,8,9]])\n",
    "my_matrix"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9827147e",
   "metadata": {},
   "source": [
    "Other useful matrix functions:\n",
    "* `np.random.randint()` builds a random matrix\n",
    "* `np.eye()` builds an **identity matrix**\n",
    "* `np.zeros()` builds a matrix of zeros\n",
    "* `np.ones()` builds a matrix of ones\n",
    "* `np.diag()` builds a **diagonal matrix**\n",
    "\n",
    "**Note**: Each of these differ in the inputs they take to instruct their size and shape -- always refer to the documentation!"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9007b5ad",
   "metadata": {},
   "source": [
    "We can also do matrix multiplication. In Python, we can use the `@` operator for matrix/vector multiplications. We can also use the NumPy [`np.dot`](https://numpy.org/doc/stable/reference/generated/numpy.dot.html#numpy.dot), [`np.matmul`](https://numpy.org/doc/stable/reference/generated/numpy.matmul.html). See the documentation for their differences."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "687476e4",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create a random matrix\n",
    "random_matrix = ...\n",
    "\n",
    "# Use matmul for matrix multiplication\n",
    "np.matmul(my_matrix,random_matrix)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0d6379e9",
   "metadata": {},
   "source": [
    "### Solving linear equations (abstract example):\n",
    "\n",
    "Let's say we have the following system of linear equations that we'd like to solve:\n",
    "\n",
    "$4x + 3y + 2z = 25$\n",
    "\n",
    "$-2x + 2y + 3z = -10$\n",
    "\n",
    "$3x - 5y + 2z = -4$\n",
    "\n",
    "We can use the inverse method, implemented using `np.linalg.inv()`. The **inverse** of a matrix (e.g., of $M$) is one that can be multiplied by it to produce an identity matrix. This is how we \"cancel\" a matrix in order to solve problems in the form $ Ax = b $."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "aa08fdca",
   "metadata": {},
   "outputs": [],
   "source": [
    "M = np.array([[4,3,2],[-2,2,3],[3,-5,2]]) # 3x3 co-efficients matrix\n",
    "K = np.array([[25],[-10],[-4]])           # column matrix\n",
    "V = np.array([['x'],['y'],['z']])         # unknown matrix"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0426ba75",
   "metadata": {},
   "source": [
    "$ M \\cdot V = K $\n",
    "\n",
    "so,\n",
    "\n",
    "$ V = M^{-1} \\cdot K$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "14563ea9",
   "metadata": {},
   "outputs": [],
   "source": [
    "V = np.linalg.inv(M).dot(K)\n",
    "V"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d7ee41b6",
   "metadata": {},
   "source": [
    "### Coding Exercise 1.1: Understanding neural transformations using linear equations\n",
    "\n",
    "We will look at a group of 2 LGN neurons which get input from 2 retinal neurons: we will call the population of LGN neurons population $p$. Below, we have the system of linear equations that dictates the neuron models for each population. $r_1$ and $r_2$ correspond to the retinal neural activities (of neuron 1 and 2). $g_{p_1}$ and  $g_{p_2}$ correspond to the responses of the LGN neurons 1 and 2 in population $p$.\n",
    "\n",
    "\\begin{align}\n",
    "r_1 + 3r_2 &= g_{p_1} \\\\\n",
    "2r_1 + r_2 &= g_{p_2}\n",
    "\\end{align}\n",
    "\n",
    "<br>\n",
    "\n",
    "<div class=\"alert alert-success\"><b>Tasks:</b>\n",
    "    \n",
    "1. Cast each equation (i.e., $g_{p_1}$ and $g_{p_2}$) as a matrix-vector multiplication:\n",
    "\n",
    "\\begin{equation}\n",
    "\\mathbf{g}_p = \\mathbf{P}\\mathbf{r}\n",
    "\\end{equation}\n",
    "\n",
    "where $P$ is the weight matrix to population $p$.\n",
    "\n",
    "2. Let's say we only recorded from the LGN cells (and know the weight matrix) and are trying to figure out how the retinal cells responded. Solve the matrix equation for the given LGN activities:\n",
    "\n",
    "\\begin{equation}\n",
    "\\mathbf{g}_p =\n",
    "\\begin{bmatrix}\n",
    "16 \\\\\n",
    "7\n",
    "\\end{bmatrix}\n",
    "\\end{equation}\n",
    "</div>\n",
    "<br>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2b972665",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create P (using np array)\n",
    "P = ...\n",
    "\n",
    "# Create g_p (using np array)\n",
    "g_p = ...\n",
    "\n",
    "# Solve for r (using np.linalg.inv)\n",
    "r = ...\n",
    "\n",
    "# Print r\n",
    "print(r)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d00e7449",
   "metadata": {},
   "source": [
    "You can recover how the retinal neurons respond given the weight matrix and LGN responses! You have solved the system of equations using matrices. We can't always do this though: let's say we have a different group of 2 LGN neurons -  population q - with the following weight matrix from the retinal neurons.\n",
    "\n",
    "\\begin{equation}Q =\n",
    "\\begin{bmatrix}\n",
    "4 & 1 \\\\\n",
    "8 & 2\n",
    "\\end{bmatrix}\n",
    "\\end{equation}\n",
    "\n",
    "As you can see if you run the next code cell, we get an error if we try to invert this matrix to solve the equation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c3de46fd",
   "metadata": {},
   "outputs": [],
   "source": [
    "g_q = np.array([16, 7])\n",
    "Q = np.array([[4, 1], [8, 2]])\n",
    "\n",
    "print(np.linalg.inv(Q) @ g_q)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "77990fce",
   "metadata": {},
   "source": [
    "## Matrices as linear transformations\n",
    "\n",
    "For now, let's start to think about all of this as **[linear transformations of matrices](https://www.youtube.com/watch?v=N6UUV9tVIr8)**. \n",
    "\n",
    "Matrices can be thought of as enacting linear transformations. When multiplied with a vector, they transform it into another vector. In fact, they are transforming a grid of space in a linear manner: the origin stays in place and grid lines remain straight, parallel, and evenly spaced."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e7389af4",
   "metadata": {},
   "source": [
    "### Coding Exercise 1.2: Creating matrices for transformations\n",
    "\n",
    "<div class=\"alert alert-success\"><b>Tasks:</b>\n",
    "\n",
    "1. Come up with a matrix $A$ for which the corresponding linear transformation is reflection through the $y$ axis (flipping across the $y$ axis). For example, $\\mathbf{x} = \\begin{bmatrix}\n",
    "2 \\\\\n",
    "6  \\\\\n",
    "\\end{bmatrix}$ should become $\\mathbf{b} = \\begin{bmatrix}\n",
    "-2 \\\\\n",
    "6  \\\\\n",
    "\\end{bmatrix}$ when multiplied with $A$.\n",
    "2. Come up with a matrix $B$ for which the corresponding linear transformation is projecting onto the $x$ axis. For example, $\\bar{x} = \\begin{bmatrix}\n",
    "2 \\\\\n",
    "3  \\\\\n",
    "\\end{bmatrix}$ should become $\\bar{b} = \\begin{bmatrix}\n",
    "2 \\\\\n",
    "0  \\\\\n",
    "\\end{bmatrix}$ when multiplied with $B$.\n",
    "\n",
    "</div>\n",
    "\n",
    "**Remember to think about where your basis vectors should end up! Then your matrix consists of the transformed basis vectors. Drawing out what you want to happen can help**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "88e6fe4a",
   "metadata": {},
   "outputs": [],
   "source": [
    "A = ...\n",
    "\n",
    "# Uncomment to visualize transformation\n",
    "# plot_linear_transformation(A)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3c299a6f",
   "metadata": {},
   "outputs": [],
   "source": [
    "B = ...\n",
    "\n",
    "# Uncomment to visualize transformation\n",
    "# plot_linear_transformation(A)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0074acc3",
   "metadata": {},
   "source": [
    "## Eigenvalues & Eigenvectors"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ddc7347c",
   "metadata": {},
   "source": [
    "[**This video**](https://www.youtube.com/watch?v=l-c7ptT7znM) covers eigenvalues and eigenvectors.\n",
    "\n",
    "Eigenvectors $\\mathbf{v}$ of a matrix $\\mathbf{W}$ are vectors that, when multipled by the matrix, equal a scalar multiple of themselves. That scalar multiple is the corresponding eigenvalue $\\lambda$.\n",
    "\n",
    "\\begin{equation}\n",
    "\\mathbf{W}\\mathbf{v} = \\lambda\\mathbf{v}\n",
    "\\end{equation}\n",
    "\n",
    "If we have one eigenvector for a matrix, we technically have an infinite amount: every vector along the span of that eigenvector is also an eigenvector. So, we often use the unit vector in that direction to summarize all the eigenvectors along that line.\n",
    "\n",
    "We can find the eigenvalues and eigenvectors of a matrix in numpy using `np.linalg.eig`."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b468b3a8",
   "metadata": {},
   "source": [
    "### Identifying transformations from eigenvectors\n",
    "\n",
    "Earlier, we learned how to think about linear transformations in terms of where the standard basis vectors end up. We can also think about them in terms of eigenvectors.\n",
    "\n",
    "Just by looking at eigenvectors before and after a transformation, **can you describe what the transformation is in words (e.g.contraction, expansion, horizontal vs vertical, projection onto an axis, reflection, and rotation)**? Try for each of the two plots below.\n",
    "\n",
    "Note that I show an eigenvector for every eigenvalue. The $x/y$ limits do not change in before vs after (so eigenvectors are showed scaled by the eigenvalues)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ff477884",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Example #1\n",
    "W = np.array([[3, 0], [0, 1]])\n",
    "plot_eig_vec_transform(W)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "dc021f60",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Example #2\n",
    "W = np.array([[0, 1], [1, 0]])\n",
    "plot_eig_vec_transform(W)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fbbba205",
   "metadata": {},
   "source": [
    "As we saw above, looking at how just the eigenvectors change after a transformation can be very informative about what that transformation was."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7ed6458d",
   "metadata": {},
   "source": [
    "<hr>\n",
    "\n",
    "## About this notebook\n",
    "Most of the content here is directly adapted from [Neuromatch Academy Materials](https://compneuro.neuromatch.io/tutorials/W0D3_LinearAlgebra/student/W0D3_Tutorial1.html), shared under a Creative Commons Attribution 4.0 International License."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}