{ "nbformat": 4, "nbformat_minor": 0, "metadata": { "colab": { "name": "Week2Tutorial2.ipynb", "provenance": [], "collapsed_sections": [], "toc_visible": true, "authorship_tag": "ABX9TyN0VWmUDsDhb1k04J7dHjJp", "include_colab_link": true }, "kernelspec": { "name": "python3", "display_name": "Python 3" } }, "cells": [ { "cell_type": "markdown", "metadata": { "id": "view-in-github", "colab_type": "text" }, "source": [ "\"Open" ] }, { "cell_type": "markdown", "metadata": { "id": "S9l0lyqLvkW7" }, "source": [ "# Tutorial 2\n", "\n", "**Linear Algebra II, Matrices**\n", "\n", "\n", "**[insert your name]**\n", "\n", "**Important reminders**: Before starting, click \"File -> Save a copy in Drive\". Produce a pdf for submission by \"File -> Print\" and then choose \"Save to PDF\".\n", "\n", "To complete this tutorial, you should have watched Videos 2.1 through 2.5.\n", "\n" ] }, { "cell_type": "markdown", "metadata": { "id": "49jtwhuUNHM8" }, "source": [ "**Credits:**\n", "The videos you watched for this week were from 3Blue1Brown. Some elements of this problem set are from or inspired by https://openedx.seas.gwu.edu/courses/course-v1:GW+EngComp4+2019/about. In particular, we are using their `plot_linear_transformation` and `plot_linear_transformations` functions, and the demonstration of the additional transformation of a matrix inverse (end of Exercise 2)" ] }, { "cell_type": "code", "metadata": { "id": "cv9HSBNPyLV9" }, "source": [ "# Imports\n", "import numpy as np\n", "import matplotlib.pyplot as plt\n", "import matplotlib\n", "import scipy.linalg\n", "\n", "# Plotting parameters\n", "matplotlib.rcParams.update({'font.size': 22})\n" ], "execution_count": null, "outputs": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Plotting functions\n" ] }, { "cell_type": "code", "metadata": { "id": "ZIdPVYl9TzmK", "cellView": "form", "tags": [ "hide-input" ] }, "source": [ "# @title Plotting functions\n", "import numpy\n", "from numpy.linalg import inv, eig\n", "from math import ceil\n", "from matplotlib import pyplot, ticker, get_backend, rc\n", "from mpl_toolkits.mplot3d import Axes3D\n", "from itertools import cycle\n", "\n", "def plot_range_3_by_2_matrix(matrix):\n", "\n", " fig = plt.figure()\n", " ax = fig.add_subplot(111, projection='3d')\n", " # Make data\n", "\n", " x = []\n", " y = []\n", " z = []\n", " for a1 in range(-10, 10, 1):\n", " for a2 in range(-10, 10, 1):\n", " vec = a1*matrix[:,0]+a2*matrix[:,1]\n", " x1, y1, z1 = vec\n", " x.append(x1)\n", " y.append(y1)\n", " z.append(z1)\n", "\n", " ax.scatter(np.array(x), np.array(y), np.array(z), color='b')\n", "\n", "def plot_null_2_by_3_matrix(matrix):\n", "\n", " fig = plt.figure()\n", " ax = fig.add_subplot(111, projection='3d')\n", " # Make data\n", "\n", " basis = scipy.linalg.null_space(matrix)\n", " x = []\n", " y = []\n", " z = []\n", " for a1 in range(-10, 10, 1):\n", " for a2 in range(-10, 10, 1):\n", "\n", " vec = a1*basis[:,0]\n", " if basis.shape[1]==2:\n", " vec += a2*basis[:,1]\n", " x1, y1, z1 = vec\n", " x.append(x1)\n", " y.append(y1)\n", " z.append(z1)\n", "\n", " ax.scatter(np.array(x), np.array(y), np.array(z), color='b')\n", "\n", "\n", "_int_backends = ['GTK3Agg', 'GTK3Cairo', 'MacOSX', 'nbAgg',\n", " 'Qt4Agg', 'Qt4Cairo', 'Qt5Agg', 'Qt5Cairo',\n", " 'TkAgg', 'TkCairo', 'WebAgg', 'WX', 'WXAgg', 'WXCairo']\n", "_backend = get_backend() # get current backend name\n", "\n", "# shrink figsize and fontsize when using %matplotlib notebook\n", "if _backend in _int_backends:\n", " fontsize = 4\n", " fig_scale = 0.75\n", "else:\n", " fontsize = 5\n", " fig_scale = 1\n", "\n", "grey = '#808080'\n", "gold = '#cab18c' # x-axis grid\n", "lightblue = '#0096d6' # y-axis grid\n", "green = '#008367' # x-axis basis vector\n", "red = '#E31937' # y-axis basis vector\n", "darkblue = '#004065'\n", "\n", "pink, yellow, orange, purple, brown = '#ef7b9d', '#fbd349', '#ffa500', '#a35cff', '#731d1d'\n", "\n", "quiver_params = {'angles': 'xy',\n", " 'scale_units': 'xy',\n", " 'scale': 1,\n", " 'width': 0.012}\n", "\n", "grid_params = {'linewidth': 0.5,\n", " 'alpha': 0.8}\n", "\n", "def set_rc(func):\n", " def wrapper(*args, **kwargs):\n", " rc('font', family='serif', size=fontsize)\n", " rc('figure', dpi=200)\n", " rc('axes', axisbelow=True, titlesize=5)\n", " rc('lines', linewidth=1)\n", " func(*args, **kwargs)\n", " return wrapper\n", "\n", "@set_rc\n", "def plot_vector(vectors, tails=None):\n", " ''' Draw 2d vectors based on the values of the vectors and the position of their tails.\n", " \n", " Parameters\n", " ----------\n", " vectors : list.\n", " List of 2-element array-like structures, each represents a 2d vector.\n", " \n", " tails : list, optional.\n", " List of 2-element array-like structures, each represents the coordinates of the tail\n", " of the corresponding vector in vectors. If None (default), all tails are set at the\n", " origin (0,0). If len(tails) is 1, all tails are set at the same position. Otherwise,\n", " vectors and tails must have the same length.\n", " \n", " Examples\n", " --------\n", " >>> v = [(1, 3), (3, 3), (4, 6)]\n", " >>> plot_vector(v) # draw 3 vectors with their tails at origin\n", " >>> t = [numpy.array((2, 2))]\n", " >>> plot_vector(v, t) # draw 3 vectors with their tails at (2,2)\n", " >>> t = [[3, 2], [-1, -2], [3, 5]]\n", " >>> plot_vector(v, t) # draw 3 vectors with 3 different tails\n", "\n", " ''' \n", " vectors = numpy.array(vectors)\n", " assert vectors.shape[1] == 2, \"Each vector should have 2 elements.\" \n", " if tails is not None:\n", " tails = numpy.array(tails)\n", " assert tails.shape[1] == 2, \"Each tail should have 2 elements.\"\n", " else:\n", " tails = numpy.zeros_like(vectors)\n", " \n", " # tile vectors or tails array if needed\n", " nvectors = vectors.shape[0]\n", " ntails = tails.shape[0]\n", " if nvectors == 1 and ntails > 1:\n", " vectors = numpy.tile(vectors, (ntails, 1))\n", " elif ntails == 1 and nvectors > 1:\n", " tails = numpy.tile(tails, (nvectors, 1))\n", " else:\n", " assert tails.shape == vectors.shape, \"vectors and tail must have a same shape\"\n", "\n", " # calculate xlimit & ylimit\n", " heads = tails + vectors\n", " limit = numpy.max(numpy.abs(numpy.hstack((tails, heads))))\n", " limit = numpy.ceil(limit * 1.2) # add some margins\n", " \n", " figsize = numpy.array([2,2]) * fig_scale\n", " figure, axis = pyplot.subplots(figsize=figsize)\n", " axis.quiver(tails[:,0], tails[:,1], vectors[:,0], vectors[:,1], color=darkblue, \n", " angles='xy', scale_units='xy', scale=1)\n", " axis.set_xlim([-limit, limit])\n", " axis.set_ylim([-limit, limit])\n", " axis.set_aspect('equal')\n", "\n", " # if xticks and yticks of grid do not match, choose the finer one\n", " xticks = axis.get_xticks()\n", " yticks = axis.get_yticks()\n", " dx = xticks[1] - xticks[0]\n", " dy = yticks[1] - yticks[0]\n", " base = max(int(min(dx, dy)), 1) # grid interval is always an integer\n", " loc = ticker.MultipleLocator(base=base)\n", " axis.xaxis.set_major_locator(loc)\n", " axis.yaxis.set_major_locator(loc)\n", " axis.grid(True, **grid_params)\n", " \n", " # show x-y axis in the center, hide frames\n", " axis.spines['left'].set_position('center')\n", " axis.spines['bottom'].set_position('center')\n", " axis.spines['right'].set_color('none')\n", " axis.spines['top'].set_color('none')\n", "\n", "@set_rc\n", "def plot_transformation_helper(axis, matrix, *vectors, unit_vector=True, unit_circle=False, title=None):\n", " \"\"\" A helper function to plot the linear transformation defined by a 2x2 matrix.\n", " \n", " Parameters\n", " ----------\n", " axis : class matplotlib.axes.Axes.\n", " The axes to plot on.\n", "\n", " matrix : class numpy.ndarray.\n", " The 2x2 matrix to visualize.\n", "\n", " *vectors : class numpy.ndarray.\n", " The vector(s) to plot along with the linear transformation. Each array denotes a vector's\n", " coordinates before the transformation and must have a shape of (2,). Accept any number of vectors. \n", " \n", " unit_vector : bool, optional.\n", " Whether to plot unit vectors of the standard basis, default to True.\n", " \n", " unit_circle: bool, optional.\n", " Whether to plot unit circle, default to False.\n", " \n", " title: str, optional.\n", " Title of the plot.\n", "\n", " \"\"\"\n", " assert matrix.shape == (2,2), \"the input matrix must have a shape of (2,2)\"\n", " grid_range = 20\n", " x = numpy.arange(-grid_range, grid_range+1)\n", " X_, Y_ = numpy.meshgrid(x,x)\n", " I = matrix[:,0]\n", " J = matrix[:,1]\n", " X = I[0]*X_ + J[0]*Y_\n", " Y = I[1]*X_ + J[1]*Y_\n", " origin = numpy.zeros(1)\n", " \n", " # draw grid lines\n", " for i in range(x.size):\n", " axis.plot(X[i,:], Y[i,:], c=gold, **grid_params)\n", " axis.plot(X[:,i], Y[:,i], c=lightblue, **grid_params)\n", " \n", " # draw (transformed) unit vectors\n", " if unit_vector:\n", " axis.quiver(origin, origin, [I[0]], [I[1]], color=green, **quiver_params)\n", " axis.quiver(origin, origin, [J[0]], [J[1]], color=red, **quiver_params)\n", "\n", " # draw optional vectors\n", " color_cycle = cycle([pink, darkblue, orange, purple, brown])\n", " if vectors:\n", " for vector in vectors:\n", " color = next(color_cycle)\n", " vector_ = matrix @ vector.reshape(-1,1)\n", " axis.quiver(origin, origin, [vector_[0]], [vector_[1]], color=color, **quiver_params)\n", "\n", " # draw optional unit circle\n", " if unit_circle:\n", " alpha = numpy.linspace(0, 2*numpy.pi, 41)\n", " circle = numpy.vstack((numpy.cos(alpha), numpy.sin(alpha)))\n", " circle_trans = matrix @ circle\n", " axis.plot(circle_trans[0], circle_trans[1], color=red, lw=0.8)\n", "\n", " # hide frames, set xlimit & ylimit, set title\n", " limit = 4\n", " axis.spines['left'].set_position('center')\n", " axis.spines['bottom'].set_position('center')\n", " axis.spines['left'].set_linewidth(0.3)\n", " axis.spines['bottom'].set_linewidth(0.3)\n", " axis.spines['right'].set_color('none')\n", " axis.spines['top'].set_color('none')\n", " axis.set_xlim([-limit, limit])\n", " axis.set_ylim([-limit, limit])\n", " if title is not None:\n", " axis.set_title(title)\n", "\n", "@set_rc\n", "def plot_linear_transformation(matrix, *vectors, unit_vector=True, unit_circle=False):\n", " \"\"\" Plot the linear transformation defined by a 2x2 matrix using the helper\n", " function plot_transformation_helper(). It will create 2 subplots to visualize some\n", " vectors before and after the transformation.\n", " \n", " Parameters\n", " ----------\n", " matrix : class numpy.ndarray.\n", " The 2x2 matrix to visualize.\n", "\n", " *vectors : class numpy.ndarray.\n", " The vector(s) to plot along with the linear transformation. Each array denotes a vector's\n", " coordinates before the transformation and must have a shape of (2,). Accept any number of vectors.\n", " \n", " unit_vector : bool, optional.\n", " Whether to plot unit vectors of the standard basis, default to True.\n", " \n", " unit_circle: bool, optional.\n", " Whether to plot unit circle, default to False.\n", " \n", " \"\"\"\n", " figsize = numpy.array([4,2]) * fig_scale\n", " figure, (axis1, axis2) = pyplot.subplots(1, 2, figsize=figsize)\n", " plot_transformation_helper(axis1, numpy.identity(2), *vectors, unit_vector=unit_vector, unit_circle=unit_circle, title='Before transformation')\n", " plot_transformation_helper(axis2, matrix, *vectors, unit_vector=unit_vector, unit_circle=unit_circle, title='After transformation')\n", "\n", "@set_rc\n", "def plot_linear_transformations(*matrices, unit_vector=True, unit_circle=False):\n", " \"\"\" Plot the linear transformation defined by a sequence of n 2x2 matrices using the helper\n", " function plot_transformation_helper(). It will create n+1 subplots to visualize some\n", " vectors before and after each transformation.\n", "\n", " Parameters\n", " ----------\n", " *matrices : class numpy.ndarray.\n", " The 2x2 matrices to visualize. Accept any number of matrices.\n", " \n", " unit_vector : bool, optional.\n", " Whether to plot unit vectors of the standard basis, default to True.\n", " \n", " unit_circle: bool, optional.\n", " Whether to plot unit circle, default to False.\n", " \n", " \"\"\"\n", " nplots = len(matrices) + 1\n", " nx = 2\n", " ny = ceil(nplots/nx)\n", " figsize = numpy.array([2*nx, 2*ny]) * fig_scale\n", " figure, axes = pyplot.subplots(nx, ny, figsize=figsize)\n", "\n", " for i in range(nplots): # fig_idx \n", " if i == 0:\n", " matrix_trans = numpy.identity(2)\n", " title = 'Before transformation'\n", " else:\n", " matrix_trans = matrices[i-1] @ matrix_trans\n", " if i == 1:\n", " title = 'After {} transformation'.format(i)\n", " else:\n", " title = 'After {} transformations'.format(i)\n", " plot_transformation_helper(axes[i//nx, i%nx], matrix_trans, unit_vector=unit_vector, unit_circle=unit_circle, title=title)\n", " # hide axes of the extra subplot (only when nplots is an odd number)\n", " if nx*ny > nplots:\n", " axes[-1,-1].axis('off')\n", " " ], "execution_count": null, "outputs": [] }, { "cell_type": "markdown", "metadata": { "id": "60wK4bTuUtPL" }, "source": [ "# Exercise 1: Don't be a square\n", "\n", "In linear algebra, square matrices get a lot of the attention. We will often be dealing with m x n non-square matrices though, especially in neuroscience contexts! In this problem set, we'll explore nonsquare matrices a bit more." ] }, { "cell_type": "markdown", "metadata": { "id": "H5yEMjS5VqZT" }, "source": [ "## A) Rank & range if n < m\n", "\n", "i) If we have a 4 x 2 matrix A, what is the maximum possible rank of this matrix? (Hint: think about the columns of A)\n", "\n", "ii) What does this mean for the range of the associated linear transformation? What is the maximum dimensionality subspace is can be?\n", "\n", "iii) If we have an m x n matrix where n < m, can we ever reach every vector in $R^m$?" ] }, { "cell_type": "markdown", "metadata": { "id": "JLXGvr67svet" }, "source": [ "**Your text answer**" ] }, { "cell_type": "markdown", "metadata": { "id": "PS6kY_k1XL_r" }, "source": [ "Play around below with inputting a 3 x 2 matrix and seeing the resulting range/column space in $R^3$. In this plot, we're taking all vectors with elements between -10 and 10 in $R^2$ and plotting the transformed vector in $R^3$.\n", "\n", "First test a random matrix. Then create a matrix that results in a 1D line." ] }, { "cell_type": "code", "metadata": { "id": "RitMNJ-CSy1w" }, "source": [ "matrix = ...\n", "\n", "plot_range_3_by_2_matrix(matrix)" ], "execution_count": null, "outputs": [] }, { "cell_type": "markdown", "metadata": { "id": "FQYf38poWnm9" }, "source": [ "## B) Rank & range if m < n\n", "\n", "i) If we have a 2 x 4 matrix A, what is the maximum possible rank of this matrix?\n", "\n", "ii) What is the minimum possible dimension of the null space for this matrix? (Hint: think about how space squishes)" ] }, { "cell_type": "markdown", "metadata": { "id": "2X7b0LsksxhY" }, "source": [ "**Your text answer**" ] }, { "cell_type": "markdown", "metadata": { "id": "Z18hBYQQXvj1" }, "source": [ "Below, we'll visualize the null space of a 2 x 3 matrix. Try out a random matrix and see what happens! Can you figure out how to make the null space 2D (Hint: take a look at Extra info: The Rank Theorem)?" ] }, { "cell_type": "code", "metadata": { "id": "yCmPQWamX606" }, "source": [ "matrix = ...\n", "plot_null_2_by_3_matrix(matrix)" ], "execution_count": null, "outputs": [] }, { "cell_type": "markdown", "metadata": { "id": "k4cMhXFEXwBc" }, "source": [ "## Extra info: Inverses of nonsquare matrices\n", "\n", "Nonsquare matrices do not have inverses. By definition, for matrix A, an inverse of A exists only if $A^{-1}A = AA^{-1} = I_n$. If you have a nonsquare matrix A, you can never get the dimensions of $A^{-1}$ to work. \n", "\n", "Left and right inverses of nonsquare matrices do exist. If $A$ is m x n and has rank equal to n, $A$ has a left inverse $B$ such that $BA = I_n$, where $B$ is n x m. If $A$ has rank m, then is has a right inverse $C$ such that $AC = I_m$, where C is n x m.\n", "\n", "Thinking through why the rank conditions in this definition are what they are is difficult but will help consolidate your knowledge if you want to try! Think about going from larger to smaller spaces vs smaller to larger spaces" ] }, { "cell_type": "markdown", "metadata": { "id": "fLcN75QeYjL1" }, "source": [ "## Extra info: Determinants of nonsquare matrices\n", "\n", "Nonsquare matrices do not have determinants. Intuitively, determinants denote how much the associated linear transformation changes the area or volume of space covered. Because nonsquare matrices are transforming between dimensions, you cannot really compute the area changes since it changes dimensions (how do you compute the difference between an area and a volume)." ] }, { "cell_type": "markdown", "metadata": { "id": "Ctk7cfcjMJ-a" }, "source": [ "## Extra info: The Rank Theorem\n", "\n", "There is a theorem that says for an m x n matrix $A$:\n", "\n", "$$rank A + dim Nul A = n$$\n", "\n", "Basically the rank of a matrix plus the dimensionality of the null space must add up to the number of columns. We will see why this is true next week once we learn briefly how to solve matrix-vector equations by hand." ] }, { "cell_type": "markdown", "metadata": { "id": "6NXkIn0tU9-Y" }, "source": [ "# (Optional) Exercise 2: Changing basis\n", "\n", "In Video 1.3: Vector Spaces, we talked about changing basis (or changing coordinate systems). We looked at transforming a vector from a basis $\\beta = \\{\\bar{b}_1, \\bar{b}_2\\}$ where $\\bar{b}_1 = \\begin{bmatrix}\n", "1 \\\\\n", "1 \\\\\n", "\\end{bmatrix}$ and $\\bar{b}_2 = \\begin{bmatrix}\n", "-1 \\\\\n", "1 \\\\\n", "\\end{bmatrix}$ to the standard basis. In particular, we had the $\\beta$ coordinates of vector $\\bar{a}$ (coordinates of a vector with respect to basis $\\beta$) as $[\\bar{a}]_{\\beta} = \\begin{bmatrix} 2 \\\\ 1 \\\\ \\end{bmatrix}$. We computed the standard coordinates as:\n", "\n", "$$ \\begin{align} \\bar{a} &= 2 \\bar{b}_1 + \\bar{b}_2 \\\\\n", "&= \\begin{bmatrix}\n", "1 \\\\\n", "3\\\\\n", "\\end{bmatrix}\n", " \\end{align}$$\n" ] }, { "cell_type": "markdown", "metadata": { "id": "HeSqoZZikGiC" }, "source": [ "## A) Changing to the standard basis\n", " Now that we have learned about matrix multiplication, we can figure out how to change basis using it! We want to write a matrix P so that:\n", "\n", " $$ \\bar{x} = P_{\\beta}[\\bar{x}]_{\\beta}$$ \n", "\n", " In other words, a matrix that transforms from the vector in coordinates relative to $\\beta$ to one with the coordinates relative to the standard basis. $P_{\\beta}$ is called the change-of-coordinates matrix from $\\beta$ to the standard basis.\n", "\n", "What is $P_{\\beta}$ for the basis given above? \n", " \n", "\n", " Hints: \n", " For the specific vector example, \n", "\n", " $$\\begin{bmatrix}\n", "1 \\\\\n", "3\\\\\n", "\\end{bmatrix} = P_{\\beta}\\begin{bmatrix}\n", "2 \\\\\n", "1\\\\\n", "\\end{bmatrix} $$\n", "\n", "You could find it by looking at the computation we did above and thinking how that could have been accomplished by a matrix multiplication or by thinking about what we learned this week: that matrices represent linear transformations and each column corresponds to where the corresponding basis vector ends up after the transformation." ] }, { "cell_type": "markdown", "metadata": { "id": "ina6kOayoFKR" }, "source": [ "**Your math answer**" ] }, { "cell_type": "markdown", "metadata": { "id": "C-Y4wt9xrOdZ" }, "source": [ "Fill in $P_{\\beta}$ below and plot the resulting linear transformation. This visualization highlights something essential: when we are in the basis $\\beta$ the $\\beta$ basis vectors seem like standard basis vectors. The Before transformation is when we are in \"$\\beta-land$\" - it looks like when we are in standard basis vector land. It is only after we transform to standard coordinates, that we see the $\\beta$ vectors defined with respect to them (After transformation)." ] }, { "cell_type": "code", "metadata": { "id": "oWpjDAdfVVVg" }, "source": [ "p_beta = ...\n", "\n", "plot_linear_transformation(p_beta)" ], "execution_count": null, "outputs": [] }, { "cell_type": "markdown", "metadata": { "id": "myTgr0JSl8Ct" }, "source": [ "## B) Changing from the standard basis\n", "\n", "Luckily for us, the matrix $P_{\\beta}$ is always invertible. We know this from the Invertible Matrix Theorem presented in Part 1 because the columns of $P_{\\beta}$ span $R^n$ (since they are a basis for $R^n$). So we can compute the coordinates of a vector with respect to $\\beta$ as:\n", "\n", " $$ [\\bar{x}]_{\\beta} = P^{-1}_{\\beta}\\bar{x}$$ \n", "\n", " We will compute $P^{-1}_{\\beta}$ using code and find the $\\beta$ coordinates of $\\bar{x} = \\begin{bmatrix}\n", "4 \\\\\n", "2\\\\\n", "\\end{bmatrix}$" ] }, { "cell_type": "code", "metadata": { "id": "nIqgnbL7oh96" }, "source": [ "p_beta_inverse = ...\n", "\n", "beta_coordinates = ...\n", "\n", "print(beta_coordinates)" ], "execution_count": null, "outputs": [] }, { "cell_type": "markdown", "metadata": { "id": "gwzcI1LuLuKF" }, "source": [ "## C) Changing basis criterion\n", "\n", "Not every matrix-vector multiplication can be considered a change of basis. What has to be true of the matrix for it to be a change of basis?\n", "\n" ] }, { "cell_type": "markdown", "metadata": { "id": "D3oy2QH0L9xa" }, "source": [ "**Your text answer**" ] }, { "cell_type": "markdown", "metadata": { "id": "jkYM079mr21N" }, "source": [ "## Extra info: Changing between two non-standard bases\n", "\n", "We can also change from basis $\\beta$ to another basis $C$ without first moving to the standard basis. We will not get into this here due to time constraints but you should have all the tools to understand the explanation here if you wish: https://math.berkeley.edu/~arash/54/notes/04_07.pdf" ] } ] }