{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Linear Algebra in NumPy\n",
    "====\n",
    "\n",
    "## Unit 9, Lecture 2\n",
    "\n",
    "*Numerical Methods and Statistics*\n",
    "\n",
    "----\n",
    "\n",
    "#### Prof. Andrew White, March 30, 2020"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "import random\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from math import sqrt, pi, erf\n",
    "import scipy.stats\n",
    "import numpy.linalg"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Working with Matrices in Numpy\n",
    "====\n",
    "\n",
    "We saw earlier in the class how to create numpy matrices. Let's review that and learn about explicit initialization"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Explicit Initialization\n",
    "----\n",
    "\n",
    "You can explitily set the values in your matrix by first creating a list and then converting it into a numpy array"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "As Python list:\n",
      "[[4, 3], [6, 2]]\n",
      "The shape of the array: (2, 2)\n",
      "The numpy matrix/array:\n",
      "[[4 3]\n",
      " [6 2]]\n"
     ]
    }
   ],
   "source": [
    "matrix = [ [4,3], [6, 2] ]\n",
    "print('As Python list:')\n",
    "print(matrix)\n",
    "\n",
    "np_matrix = np.array(matrix)\n",
    "\n",
    "print('The shape of the array:', np.shape(np_matrix))\n",
    "print('The numpy matrix/array:')\n",
    "print(np_matrix)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "You can use multiple lines in python to specify your list. This can make the formatting cleaner"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 4  3]\n",
      " [ 1  2]\n",
      " [-1  4]\n",
      " [ 4  2]]\n"
     ]
    }
   ],
   "source": [
    "np_matrix_2 = np.array([\n",
    "                        [ 4, 3],   \n",
    "                        [ 1, 2],  \n",
    "                        [-1, 4], \n",
    "                        [ 4, 2]   \n",
    "                        ])\n",
    "print(np_matrix_2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Create and Set\n",
    "----\n",
    "\n",
    "You can also create an array and then set the elements."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]\n",
      " [0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]]\n"
     ]
    }
   ],
   "source": [
    "np_matrix_3 = np.zeros( (2, 10) )\n",
    "\n",
    "print(np_matrix_3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[0. 2. 0. 0. 0. 0. 0. 0. 0. 0.]\n",
      " [0. 2. 0. 0. 0. 0. 0. 0. 0. 0.]]\n"
     ]
    }
   ],
   "source": [
    "np_matrix_3[:, 1] = 2\n",
    "print(np_matrix_3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[-1. -1. -1. -1. -1. -1. -1. -1. -1. -1.]\n",
      " [ 0.  2.  0.  0.  0.  0.  0.  0.  0.  0.]]\n"
     ]
    }
   ],
   "source": [
    "np_matrix_3[0, :] = -1\n",
    "print(np_matrix_3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[-1. -1. -1. -1. -1. -1. -1. -1. -1. -1.]\n",
      " [ 0.  2.  0.  0.  0.  0. 43.  0.  0.  0.]]\n"
     ]
    }
   ],
   "source": [
    "np_matrix_3[1, 6] = 43\n",
    "print(np_matrix_3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[-1. -1. -1. -1. -1. -1. -1. -1. -1. -1.]\n",
      " [ 0.  1.  4.  9. 16. 25. 36. 49. 64. 81.]]\n"
     ]
    }
   ],
   "source": [
    "rows, columns = np.shape(np_matrix_3) #get the number of rows and columns\n",
    "for i in range(columns): #Do a for loop over columns\n",
    "    np_matrix_3[1, i] = i ** 2 #Set the value of the 2nd row, ith column to be i^2\n",
    "    \n",
    "print(np_matrix_3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Linear Algebra\n",
    "====\n",
    "\n",
    "The linear algebra routines for python are in the `numpy.linalg` library. [See here](http://docs.scipy.org/doc/numpy/reference/routines.linalg.html)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Matrix Multiplication\n",
    "---\n",
    "\n",
    "Matrix multiplication is done with the `dot` method. Let's compare that with `*`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[2.02700017]\n",
      " [1.20440067]]\n"
     ]
    }
   ],
   "source": [
    "np_matrix_1 = np.random.random( (2, 4) ) #create a random 2 x 4 array\n",
    "np_matrix_2 = np.random.random( (4, 1) ) #create a random 4 x 1 array\n",
    "\n",
    "print(np_matrix_1.dot(np_matrix_2))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "So, dot correctly gives us a 2x1 matrix as expected for the two shapes"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Using the special `@` character:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[2.02700017]\n",
      " [1.20440067]]\n"
     ]
    }
   ],
   "source": [
    "print(np_matrix_1 @ np_matrix_2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "The element-by-element multiplication, `*`, doesn't work on different sized arrays."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "ename": "ValueError",
     "evalue": "operands could not be broadcast together with shapes (2,4) (4,1) ",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mValueError\u001b[0m                                Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-11-f16136fbd191>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[0;32m----> 1\u001b[0;31m \u001b[0mprint\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mnp_matrix_1\u001b[0m \u001b[0;34m*\u001b[0m \u001b[0mnp_matrix_2\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m",
      "\u001b[0;31mValueError\u001b[0m: operands could not be broadcast together with shapes (2,4) (4,1) "
     ]
    }
   ],
   "source": [
    "print(np_matrix_1 * np_matrix_2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Method vs Function\n",
    "----\n",
    "\n",
    "Instead of using `dot` as a method (it comes after a `.`), you can use the `dot` function as well. Let's see an example:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[2.02700017]\n",
      " [1.20440067]]\n",
      "[[2.02700017]\n",
      " [1.20440067]]\n"
     ]
    }
   ],
   "source": [
    "print(np_matrix_1.dot(np_matrix_2))\n",
    "\n",
    "print(np.dot(np_matrix_1, np_matrix_2))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Matrix Rank\n",
    "----\n",
    "\n",
    "The rank of a matrix can be found with singular value decomposition. In numpy, we can do this simply with a call to `linalg.matrix_rank`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1\n"
     ]
    }
   ],
   "source": [
    "import numpy.linalg as linalg\n",
    "\n",
    "matrix =  [ [1, 0], [0, 0] ]\n",
    "np_matrix = np.array(matrix)\n",
    "\n",
    "print(linalg.matrix_rank(np_matrix))\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Matrix Inverse\n",
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "The inverse of a matrix can be found using the `linalg.inverse` command. Consider the following system of equations:\n",
    "\n",
    "$$\\begin{array}{lr}\n",
    "3 x + 2 y + z & = 5\\\\\n",
    "2 x - y & = 4 \\\\\n",
    "x + y - 2z & = 12 \\\\\n",
    "\\end{array}$$\n",
    "\n",
    "We can encode it as a matrix equation:\n",
    "\n",
    "$$\\left[\\begin{array}{lcr}\n",
    "3 & 2 & 1\\\\\n",
    "2 & -1 & 0\\\\\n",
    "1 & 1 & -2\\\\\n",
    "\\end{array}\\right]\n",
    "\\left[\\begin{array}{l}\n",
    "x\\\\\n",
    "y\\\\\n",
    "z\\\\\n",
    "\\end{array}\\right]\n",
    "=\n",
    "\\left[\\begin{array}{l}\n",
    "5\\\\\n",
    "4\\\\\n",
    "12\\\\\n",
    "\\end{array}\\right]$$\n",
    "\n",
    "$$\\mathbf{A}\\mathbf{x} = \\mathbf{b}$$\n",
    "\n",
    "$$\\mathbf{A}^{-1}\\mathbf{b} = \\mathbf{x}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[ 2.47058824  0.94117647 -4.29411765]\n",
      "[ 5.  4. 12.]\n"
     ]
    }
   ],
   "source": [
    "\n",
    "#Enter the data as lists\n",
    "a_matrix = [[3, 2, 1],\n",
    "            [2,-1,0],\n",
    "            [1,1,-2]]\n",
    "b_matrix = [5, 4, 12]\n",
    "\n",
    "#convert them to numpy arrays/matrices\n",
    "np_a_matrix = np.array(a_matrix)\n",
    "np_b_matrix = np.array(b_matrix).transpose()\n",
    "\n",
    "#Solve the problem\n",
    "np_a_inv = linalg.inv(np_a_matrix)\n",
    "np_x_matrix = np_a_inv @ np_b_matrix\n",
    "\n",
    "#print the solution\n",
    "print(np_x_matrix)\n",
    "\n",
    "#check to make sure the answer works\n",
    "print(np_a_matrix @ np_x_matrix)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Computation cost for inverse\n",
    "---\n",
    "\n",
    "Computing a matrix inverse can be VERY expensive for large matrices. Do not exceed about 500 x 500 matrices"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Eigenvectors/Eigenvalues\n",
    "----\n",
    "\n",
    "Before trying to understand what an eigenvector is, let's try to understand their analogue, a stationary point."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "A stationary point of a function $f(x)$ is an $x$ such that:\n",
    "\n",
    "$$x = f(x)$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Consider this function:\n",
    "\n",
    "$$f(x) = x - \\frac{x^2 - 612}{2x}$$\n",
    "\n",
    "If we found a stationary point, that would be mean that\n",
    "\n",
    "$$x = x - \\frac{x^2 - 612}{2x} $$\n",
    "\n",
    "or \n",
    "\n",
    "$$ x^2 = 612 $$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "More generally, you can find a square root of $A$ by finding a stationary point to:\n",
    "\n",
    "$$f(x) = x - \\frac{x^2 - A}{2x} $$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "In this case, you can find the stationary point by just doing $x_{i+1} = f(x_i)$ until you are stationary"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "x = 1\n",
    "for i in range(10):\n",
    "    x = x - (x**2 - 612) / (2 * x)\n",
    "    print(i, x)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Eigenvectors/Eigenvalues\n",
    "----\n",
    "\n",
    "Matrices are analogues of functions. They take in a vector and return a vector.\n",
    "$$\\mathbf{A}\\mathbf{x} = \\mathbf{y}$$\n",
    "\n",
    "Just like stationary points, there is sometimes a special vector which has this property:\n",
    "\n",
    "$$\\mathbf{A}\\mathbf{x} = \\mathbf{x}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Such a vector is called an **eigenvector**. It turns out such vectors are rarely always exists. If we instead allow a scalar, we can find a whole bunch like this:\n",
    "\n",
    "$$\\mathbf{A}\\mathbf{v} =\\lambda\\mathbf{v}$$\n",
    "\n",
    "These are like the stationary points above, except we are getting back our input times a constant. That means it's a particular *direction* that is unchanged, not the value. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Finding Eigenvectors/Eigenvalues\n",
    "---\n",
    "\n",
    "Eigenvalues/eigenvectors can be found easily as well in python, including for complex numbers and sparse matrices. The command `linalg.eigh` will return only the real eigenvalues/eigenvectors. That assumes your matrix is Hermitian, meaning it is symmetric (if your matrix is real numbers). Use `eig` to get general possibly complex eigenvalues Here's an easy example:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Let's consider this matrix: \n",
    "\n",
    "$$\n",
    "A = \\left[\\begin{array}{lr}\n",
    "3 & 1\\\\\n",
    "1 & 3\\\\\n",
    "\\end{array}\\right]$$\n",
    "\n",
    "Imagine it as a geometry operator. It takes in a 2D vector and morphs it into another 2D vector."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "$$\\vec{x} = [1, 0]$$\n",
    "\n",
    "$$A \\,\\vec{x}^T = [3, 1]^T$$\n",
    "\n",
    "Now is there a particular *direction* where $\\mathbf{A}$ cannot affect it?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "A = np.array([[3,1], [1,3]])\n",
    "\n",
    "e_values, e_vectors = np.linalg.eig(A)\n",
    "\n",
    "print(e_vectors)\n",
    "print(e_values)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "So that means $v_1 = [0.7, 0.7]$ and $v_2 = [-0.7, 0.7]$. Let's find out:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "v1 = e_vectors[:,0]\n",
    "v2 = e_vectors[:,1]\n",
    "\n",
    "A @ v1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Yes, that is the same direction! And notice that it's 4 times as much as the input vector, which is what the eigenvalue is telling us."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "A random matrix will almost never be Hermitian, so look out for complex numbers. In engineering, your matrices commonly be Hermitian. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "A = np.random.normal(size=(3,3))\n",
    "e_values, e_vectors = linalg.eig(A)\n",
    "print(e_values)\n",
    "print(e_vectors)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "Notice that there are compelx eigenvalues, so `eigh` was not correct to use"
   ]
  }
 ],
 "metadata": {
  "celltoolbar": "Slideshow",
  "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": 1
}