00 - Linear Algebra.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Session 0: Linear Algebra"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This notebook is a simple translation of the linear algebra revision lecture into Python constructs, using numpy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.random.seed(42)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Please note that Python uses 0-based indexing, as does the rest of this notebook"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Basic Notation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$A \\in \\mathbb{R}^{m x n}$ is a matrix with m rows and n columns"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 0.49671415, -0.1382643 ,  0.64768854],\n",
       "       [ 1.52302986, -0.23415337, -0.23413696],\n",
       "       [ 1.57921282,  0.76743473, -0.46947439],\n",
       "       [ 0.54256004, -0.46341769, -0.46572975],\n",
       "       [ 0.24196227, -1.91328024, -1.72491783]])"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "m = 5\n",
    "n = 3\n",
    "A = np.random.randn(m, n)\n",
    "A"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$x \\in \\mathbb{R}^{n}$ is a vector with n rows"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([-0.56228753, -1.01283112,  0.31424733])"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x = np.random.randn(n)\n",
    "x"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$a_{ij}$ is the entry in the ith row and jth column of A"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.76743472915290878"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "i = 2\n",
    "j = 1\n",
    "a_ij = A[i][j]\n",
    "a_ij"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Matrix Multiplication"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The product of two matrixes $A \\in \\mathbb{R}^{m x n}$ and $B \\in \\mathbb{R}^{n x p}$ is the matrix \n",
    "\n",
    "$C = AB \\in \\mathbb{R}^{m x p}$ where \n",
    "\n",
    "$C_{ij} = \\sum_{k=1}^{n} A_{ik} B_{kj}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[-0.90802408, -1.4123037 ,  1.46564877, -0.2257763 ,  0.0675282 ,\n",
       "        -1.42474819],\n",
       "       [-0.54438272,  0.11092259, -1.15099358,  0.37569802, -0.60063869,\n",
       "        -0.29169375],\n",
       "       [-0.60170661,  1.85227818, -0.01349722, -1.05771093,  0.82254491,\n",
       "        -1.22084365]])"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "p = 6\n",
    "B = np.random.randn(n, p)\n",
    "B"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[-0.76547819,  0.48285148,  0.87840781, -0.84915915,  0.64934202,\n",
       "        -1.45808819],\n",
       "       [-1.11459697, -2.61064038,  2.50489606, -0.18418579,  0.05090089,\n",
       "        -1.81578833],\n",
       "       [-1.56925562, -3.01479942,  1.43759548,  0.42834307, -0.74047335,\n",
       "        -1.90068169],\n",
       "       [ 0.03985168, -1.68032411,  1.33487931,  0.19600514, -0.06809894,\n",
       "        -0.06925218],\n",
       "       [ 1.85974361, -3.74897788,  2.58008658,  1.0510195 , -0.25329297,\n",
       "         2.31921156]])"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "C = np.matmul(A, B)\n",
    "C"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[-0.76547819,  0.48285148,  0.87840781, -0.84915915,  0.64934202,\n",
       "        -1.45808819],\n",
       "       [-1.11459697, -2.61064038,  2.50489606, -0.18418579,  0.05090089,\n",
       "        -1.81578833],\n",
       "       [-1.56925562, -3.01479942,  1.43759548,  0.42834307, -0.74047335,\n",
       "        -1.90068169],\n",
       "       [ 0.03985168, -1.68032411,  1.33487931,  0.19600514, -0.06809894,\n",
       "        -0.06925218],\n",
       "       [ 1.85974361, -3.74897788,  2.58008658,  1.0510195 , -0.25329297,\n",
       "         2.31921156]])"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Sanity check:\n",
    "C_manual = np.zeros((m, p))\n",
    "for i in range(m):\n",
    "    for j in range(p):\n",
    "        for k in range(n):\n",
    "            C_manual[i][j] += A[i][k] * B[k][j]\n",
    "\n",
    "C_manual"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(5, 3)\n",
      "(3, 6)\n",
      "(5, 6)\n"
     ]
    }
   ],
   "source": [
    "print(A.shape)\n",
    "print(B.shape)\n",
    "print(C.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Vector-Vector Products"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Given two vectors $x, y \\in \\mathbb{R}^{n}$, $x^{T}y$, called the inner product or dot product, is a real number given by:\n",
    "\n",
    "$x^{T}y \\in \\mathbb{R} = \\sum_{i=1}^{n} x_{i} y_{i}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[ 0.2088636  -1.95967012 -1.32818605]\n",
      "[ 0.19686124  0.73846658  0.17136828]\n"
     ]
    }
   ],
   "source": [
    "x = np.random.randn(n)\n",
    "y = np.random.randn(n)\n",
    "print(x)\n",
    "print(y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-1.6336427091601982"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.dot(x, y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-1.6336427091601982"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Sanity check\n",
    "dot_prod = 0\n",
    "for i in range(n):\n",
    "    dot_prod += x[i] * y[i]\n",
    "    \n",
    "dot_prod"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Given $x \\in \\mathbb{R}^{m}, y \\in \\mathbb{R}^{n}$\n",
    "\n",
    "$xy^T \\in \\mathbb{R}^{m x n}$ is called the outer product of the vectors."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[-0.11564828 -0.3011037  -1.47852199 -0.71984421 -0.46063877]\n",
      "[ 1.05712223  0.34361829 -1.76304016]\n"
     ]
    }
   ],
   "source": [
    "x = np.random.randn(m)\n",
    "y = np.random.randn(n)\n",
    "print(x)\n",
    "print(y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[-0.12225437, -0.03973886,  0.20389257],\n",
       "       [-0.31830341, -0.10346474,  0.53085791],\n",
       "       [-1.56297846, -0.5080472 ,  2.60669364],\n",
       "       [-0.76096331, -0.24735164,  1.26911425],\n",
       "       [-0.48695148, -0.15828391,  0.81212465]])"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.outer(x, y)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Matrix-Vector Products"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Given a matrix $A \\in \\mathbb{R}^{m x n}$ and a vector $x \\in \\mathbb{R}^{n}$, their product is a vector \n",
    "\n",
    "$y = Ax \\in \\mathbb{R}^{m}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(5, 3)"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A.shape # m x n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 0.32408397, -0.38508228, -0.676922  ])"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x = np.random.randn(A.shape[1])\n",
    "x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(3,)"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x.shape # n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([-0.22421439,  0.74225033,  0.53406958,  0.66955167,  1.98282124])"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "y = np.matmul(A, x)\n",
    "y"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(5,)"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "y.shape # should be m"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Matrix-Matrix Products"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This was already covered above. Some basic properties\n",
    "\n",
    "  * Matrix multiplication is associative, i.e. $ABC = (AB)C = A(BC)$\n",
    "  * Matrix multiplication is distributive, i.e. $A(B + C) = AB + AC$\n",
    "  * Matrix multiplication is **not** commutative, i.e. $AB \\ne BA$ (in general)  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Operations and Properties"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### The Identity Matrix and Diagonal Matrices"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The identity matrix $I_{n} \\in \\mathbb{R}^{n x n}$ is a square matrix with ones on the diagonal and zeros everywhere else"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 1.,  0.,  0.],\n",
       "       [ 0.,  1.,  0.],\n",
       "       [ 0.,  0.,  1.]])"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "I = np.identity(n)\n",
    "I"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It has the property that for all $A \\in \\mathbb{R}^{m x n}$:\n",
    "    \n",
    "$AI_{n} = A = I_{m}A$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 0.49671415, -0.1382643 ,  0.64768854],\n",
       "       [ 1.52302986, -0.23415337, -0.23413696],\n",
       "       [ 1.57921282,  0.76743473, -0.46947439],\n",
       "       [ 0.54256004, -0.46341769, -0.46572975],\n",
       "       [ 0.24196227, -1.91328024, -1.72491783]])"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 0.49671415, -0.1382643 ,  0.64768854],\n",
       "       [ 1.52302986, -0.23415337, -0.23413696],\n",
       "       [ 1.57921282,  0.76743473, -0.46947439],\n",
       "       [ 0.54256004, -0.46341769, -0.46572975],\n",
       "       [ 0.24196227, -1.91328024, -1.72491783]])"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.matmul(A, np.identity(A.shape[1]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 0.49671415, -0.1382643 ,  0.64768854],\n",
       "       [ 1.52302986, -0.23415337, -0.23413696],\n",
       "       [ 1.57921282,  0.76743473, -0.46947439],\n",
       "       [ 0.54256004, -0.46341769, -0.46572975],\n",
       "       [ 0.24196227, -1.91328024, -1.72491783]])"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.matmul(np.identity(A.shape[0]), A)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A diagonal matrix is a matrix where all non-diagonal elements are 0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 0.61167629,  0.        ,  0.        ],\n",
       "       [ 0.        ,  1.03099952,  0.        ],\n",
       "       [ 0.        ,  0.        ,  0.93128012]])"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "D = np.identity(n) * np.random.randn(n)\n",
    "D"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### The Transpose"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The transpose of a matrix results from \"flipping\" the rows and columns"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 0.49671415,  1.52302986,  1.57921282,  0.54256004,  0.24196227],\n",
       "       [-0.1382643 , -0.23415337,  0.76743473, -0.46341769, -1.91328024],\n",
       "       [ 0.64768854, -0.23413696, -0.46947439, -0.46572975, -1.72491783]])"
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A.T"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A square matrix $A \\in \\mathbb{R}^{n x n}$ is symmetric if $A = A^T$, and anti-symmetric if $A = -A^T$\n",
    "\n",
    "For any matrix $A \\in \\mathbb{R}^{n x n}$, \n",
    "\n",
    "$A + A^T$ is symmetric and $A - A^T$ is anti-symmetric."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[-0.83921752, -0.30921238,  0.33126343],\n",
       "       [ 0.97554513, -0.47917424, -0.18565898],\n",
       "       [-1.10633497, -1.19620662,  0.81252582]])"
      ]
     },
     "execution_count": 26,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A = np.random.randn(n, n)\n",
    "A"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[-1.67843505,  0.66633275, -0.77507154],\n",
       "       [ 0.66633275, -0.95834848, -1.3818656 ],\n",
       "       [-0.77507154, -1.3818656 ,  1.62505164]])"
      ]
     },
     "execution_count": 27,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A + A.T"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 0.        , -1.2847575 ,  1.43759841],\n",
       "       [ 1.2847575 ,  0.        ,  1.01054765],\n",
       "       [-1.43759841, -1.01054765,  0.        ]])"
      ]
     },
     "execution_count": 28,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A - A.T"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It follows that any matrix can be represented as a sum of a symmetric and anti-symmetric matrix:\n",
    "\n",
    "$A = \\frac{1}{2}(A + A^T) + \\frac{1}{2}(A - A^T)$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[-0.83921752, -0.30921238,  0.33126343],\n",
       "       [ 0.97554513, -0.47917424, -0.18565898],\n",
       "       [-1.10633497, -1.19620662,  0.81252582]])"
      ]
     },
     "execution_count": 29,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "0.5*(A + A.T) + 0.5*(A - A.T)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### The Trace"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The trace of a square matrix $A \\in \\mathbb{R}^{n x n}$, denoted as tr(A) is the sum of the diagonal elements:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-0.50586593867373053"
      ]
     },
     "execution_count": 30,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A.trace()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-0.50586593867373053"
      ]
     },
     "execution_count": 31,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "tr = 0\n",
    "for i in range(n):\n",
    "    tr += A[i][i]\n",
    "    \n",
    "tr"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Norms"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The norm of a vector $||x||$ is informally a measure of the length of the vector.\n",
    "\n",
    "For example, we have the Euclidean or $l_{2}$ norm:\n",
    "\n",
    "$||x||_{2} = \\sqrt{\\sum_{i=1}^{n} x_{i}^2}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 0.32408397, -0.38508228, -0.676922  ])"
      ]
     },
     "execution_count": 32,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.8435295942689488"
      ]
     },
     "execution_count": 33,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.linalg.norm(x, ord=2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.8435295942689488"
      ]
     },
     "execution_count": 34,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.sqrt(np.dot(x,x.T))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Other examples are $l_1$ norm:\n",
    "\n",
    "$||x||_{1} = \\sum_{i=1}^{n} |x_{i}|$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1.3860882501170702"
      ]
     },
     "execution_count": 35,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.linalg.norm(x, ord=1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "and the infinity norm:\n",
    "\n",
    "$||x||_{\\infty} = max_{i} |x_{i}|$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.67692200030595873"
      ]
     },
     "execution_count": 36,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.linalg.norm(x, ord=np.inf)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "All of the above are from the family of $l_p$ norms, which are parameterized by $p \\geq 1$, and defined as \n",
    "\n",
    "$$||x_p|| = \\left( \\sum_{i=1}^{n} |x_{i}^p| \\right)^\\frac{1}{p}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Another important one is the Frobenius norm, for matrices\n",
    "\n",
    "$||A||_{F} = \\sqrt{tr(A^T A)}$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Linear Independence and Rank"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A set of vectors ${x_1, x_2, x_3...}$ are said to be linearly independent if no vector can be represented as a linear combination of the remaining vectors\n",
    "\n",
    "The column rank of a matrix is the size of the largest subset of columns of A that constitute a linearly independent set"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {},
   "outputs": [],
   "source": [
    "A = np.array([[1, 2, 3], [4, 5, 6], [9, 12, 15]])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 1,  2,  3],\n",
       "       [ 4,  5,  6],\n",
       "       [ 9, 12, 15]])"
      ]
     },
     "execution_count": 38,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2"
      ]
     },
     "execution_count": 39,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.linalg.matrix_rank(A)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### The Inverse"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The inverse of a square matrix $A \\in \\mathbb{R}^{n x n}$ is denoted by $A^{-1}$, and is the unique matrix s.t. \n",
    "\n",
    "$A^{-1} A = I = AA^{-1}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 1.35624003, -0.07201012,  1.0035329 ],\n",
       "       [ 0.36163603, -0.64511975,  0.36139561],\n",
       "       [ 1.53803657, -0.03582604,  1.56464366]])"
      ]
     },
     "execution_count": 40,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A = np.random.randn(n, n)\n",
    "A"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 2.70796683, -0.20849188, -1.68868211],\n",
       "       [ 0.02715442, -1.57233199,  0.3457551 ],\n",
       "       [-2.66129553,  0.16894435,  2.30700572]])"
      ]
     },
     "execution_count": 41,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A_inv = np.linalg.inv(A)\n",
    "A_inv"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[  1.00000000e+00,   5.05780563e-17,   0.00000000e+00],\n",
       "       [ -1.96993814e-16,   1.00000000e+00,  -1.11022302e-16],\n",
       "       [  0.00000000e+00,   0.00000000e+00,   1.00000000e+00]])"
      ]
     },
     "execution_count": 42,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.matmul(A, A_inv)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[  1.00000000e+00,   3.71549568e-17,  -8.88178420e-16],\n",
       "       [ -6.15912127e-17,   1.00000000e+00,  -1.11022302e-16],\n",
       "       [  4.44089210e-16,  -1.38777878e-17,   1.00000000e+00]])"
      ]
     },
     "execution_count": 43,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.matmul(A_inv, A)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We say that $A$ is invertible or non-singular if $A^{-1}$ exists and non-invertible or singular otherwise"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In order for a square matrix A to have an inverse, A must be full rank."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Properties**\n",
    "\n",
    "  * $(A^{-1})^{-1} = A$\n",
    "  * $(AB)^{-1} = B^{-1}A^{-1}$\n",
    "  * $(A^{-1})^T = (A^T)^{-1}$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Orthogonal Matrices"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Two vectors $x, y \\in \\mathbb{R}^{n}$ are orthogonal if $x^T y = 0$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = np.array([5, 0, 0]) # vector along x-axis\n",
    "y = np.array([0, 4, 0]) # vector along y-axis"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[5 0 0]\n",
      "[0 4 0]\n"
     ]
    }
   ],
   "source": [
    "print(x)\n",
    "print(y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0"
      ]
     },
     "execution_count": 46,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.dot(x.T, y)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A vector $x$ is normalized if $||x||_2 = 1$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = np.array([1./np.sqrt(2), 1./np.sqrt(2)])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.99999999999999989"
      ]
     },
     "execution_count": 48,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.linalg.norm(x)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A square matrix U is orthogonal if all its columns are orthogonal to each other and are normalized"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 0.70710678,  0.70710678],\n",
       "       [-0.70710678,  0.70710678]])"
      ]
     },
     "execution_count": 49,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "U = np.array([[1./np.sqrt(2), 1./np.sqrt(2)], [-1./np.sqrt(2), 1./np.sqrt(2)]])\n",
    "U"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1.0\n",
      "1.0\n"
     ]
    }
   ],
   "source": [
    "print(np.linalg.norm(U[:, 0]))\n",
    "print(np.linalg.norm(U[:, 1]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.0"
      ]
     },
     "execution_count": 51,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.dot(U[:, 0], U[:, 1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It follows immediately from the definition of orthogonality and normality that \n",
    "\n",
    "$U^T U = I = U U^T$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 52,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 1.,  0.],\n",
       "       [ 0.,  1.]])"
      ]
     },
     "execution_count": 52,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.matmul(U.T, U)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 53,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 1.,  0.],\n",
       "       [ 0.,  1.]])"
      ]
     },
     "execution_count": 53,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.matmul(U, U.T)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can think about orthogonal matrices as a rotated version of the unit-vectors along the axes in an n-dimensional space.\n",
    "\n",
    "For example, the above U is just 2 vectors along the x and y axes in 2D space, rotated at an angle of 45 degrees."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 54,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 55,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0.5,1,'Visualizing Orthogonal Matrices as Rotated Axes')"
      ]
     },
     "execution_count": 55,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAUwAAAE/CAYAAAAt2PowAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAIABJREFUeJzsnXl4VdW5/z9vEkICiSSG2SBBIIiCoIhDHQutgtqiP+fbq+LVIlatVeutw61VW1urrbW9TrXWqrdWq7YqtVCoA3WqFdBoIpEQIGhMZAiEJCQh0/r9sXbCIWQ45+xzzl6vPZ/nyZNz9rD2d79n7Xev4V1riTGGJEmSJEnSPylBC0iSJEkSLSQdZpIkSZKESdJhJkmSJEmYJB1mkiRJkoRJ0mEmSZIkSZgkHWaSJEmShIkKhykiD4nI9+N8jeUicqn3+RsisiyMc5aIyEXx1OUXESkQESMiaUFrCRcROVFEKoPWASAiN4nII0Hr+CKjKY8G7jBFZKmI3N7D9nki8rmIpBljFhpjfpgoTcaYJ40xJ4Vx3FxjzOPx0CAiOSLyoGeDRhEpFpGLwzivQkS+Eg9NruA9XJtCHzARSRORzSISVmBxuE7ZGPNjY8ylfvQmGu/eOkSkQUTqRWRNOHkn5PyuwkOYx98qIr+PTm34eLq2i8jAeF+rNwJ3mMBjwAUiIt22XwA8aYxpS7ykYBGRdOBlYCxwNDAEuB64U0Su7eUc59/OMaYWmBvy/RRgeywvoNymVcaYLGAf4BrgNyIyKWBNUSMiBcBxgAG+HpgQY0ygf0AmsAM4PmRbLtAMTPO+Pwb8yPs8FHgJ+8BsA94AUrx9BpgQkk7oebneeVuwD9ZLQH7IscuBS73P84E3vc//DTSE/LUCj/V2DvAzL/0NwNyQ9McBrwP1WGd4P/D7XmxyCbAZGNxt+7mehn287xXA94APgV3AU0AH0OQd999AgWeXi4BPgK3AzSFpDgTuBaq8v3uBgSH7/xuo9vZdGmpjrCN/wrPpRuB/Qn6L/uxxMVDq2WM9cFnIvhOByj7yjPGu9WzItueAmwHT3zWAwZ6NOkJ+19HArV46vwfqvPu9NfR3Ao4F3sbmv0+B+SF2/Jln403AQ0Bmf3m2h3v7pZduHbAKOC5k3xHASm/fJuCeXtLYy37Y/HR2yPcvASuwz94K4Eve9juAduzz1wDc15cuYA7Qgn0uGoAPQvLGb7288xnwIyDV25fq2Wqr97tc4f2maX385rcAbwH3AC+FbE8HioCrQtJ+C7jF+54C3ACsA2qAZ4B9vX0Z3m9d4/02K4ARffqroB2mJ/w3wCMh3y8DinpxfD/xMuMA7+84QMJwmHnAmcAgIBt4FnihP4fZTecYrOM4pZdzWoFvej/a5d6xndr+6WWSdOxDV0fvDvNp4PEetqcBbcDJIQ6zyNOVGbLtKyHnFHh2+Q325TQN61wne/tvB94BhgPDsM7ghyEPw+fAwZ7d/o89HeYTwIuePQuAMuCSMO1xKjAeEOAEoBE4LAKHOQXrNHK8v03etlCHGdE1sM6xFTgd+6BlEuIwgf2xzvd8bN7LA6Z7++4FFgH7evb4C/CT/vJsD/f2n166acB1nv0zQvLQBd7nLOCo/hymdx9fx74cDvW27Yt9iV3gXed873te93wdpq4uG4Uc/wLwa+zLaTjwLrtfWAuBj7H5dl/gNfp3mOXAt4AZ3m80ImTfFE//ZOxL8x12O+fveN/zsS+1XwNPhfiZv2DzdqqX9j4aHOax2Ddd50P/FnBNL47vduxDOqGHdHp1mD0cOx3YHvK9K5PQg8PEPjyrgO/1cU55yL5Bnp6R2AetDRgUsv/33TNZyL6XgTt72fc58A3vcwXwX932V9CzwwwtTb8LnOd9Xof3AvC+nwxUeJ8fxXvove8TOm3sZbBdwEEh+y8Dlvdnj17u6wXg6u4PfC/Hdmp4xLvmQuwLYQIhDjPSa2Af/Nd72NbpMG8Enu8hXQF2AuNDth0NbOgvz4bxbGxnd03rdeA2YGg/55yIdZC13m/UDnwnZP8FwLvdzvknu0vLy+nmMPvR1WUj7/sI77qZIdvOB17zPr8KLAzZdxJ9OEysf2jtvG+ss72m2zHXedu3AxNDtpcCs0O+j/LSSgP+C1tAOCTc38OFNkyMMW9iq3XzROQAYCbwh14Ovxv7tlkmIutF5IZwriEig0Tk1yKyUUTqsJkvR0RSw5T5W2CNMeanfRzzeecHY0yj9zELW93bFrINbPWmN7Zif9ju95CGrd5tDTOdHrVhS1pZ3ufR2Op0Jxu9bZ37QtMP/TwUW1rufu5+PV2zmz0Qkbki8o6IbBORWmwb5NAw76WTJ4ALvb8nuu+M8hp92XMM9gXTnWHYF8IqEan1rvU3bztEkGdF5DoRKRWRHV46Q0I0XwIUAh+LyAoROa0PrVXGmBxsG+avgFkh+7r/5rD3bxeJru6MxZakq0Ps8WtsSbPz+qF27q6lOxcBy4wxnfn+D962UB7HFg4WG2PWdtPyfIiOUuwLZAS2xrQUeFpEqkTkLhEZ0JcQJxymR2fmvwBrnE09HWSMqTfGXGeMOQD4GnCtiMz2djdiM24nI0M+XwdMAo40xuwDHO9t797ZtBdeBp+EzbDRUA3sKyKh2sb0cfzLwFwRGdxt+5nYN/c7IdtMt2O6f++PKmym6mR/bxtY3fkh+0I1b8W+qbuf+1l/F/R6Of+EbaIY4T3Yiwnjt+jGG9gXywhse2kk1+jNTn3Z71NsFb87W7FtogcbY3K8vyHGdrr0l2dDNR+HbZM+B8j1NO/o1GyMWWuMOR/reH4KPNdDHtnzZozZ5aU5VURO9zZ3/81hz99uDxv0p6v78Vg77cKWCDvtsY8x5mBvfzV75qX9e9MvIpnedU/wIkY+x3ZiTRORaSGHPoBtJz5ZRI7tpmVuiI4cY0yGMeYzY0yrMeY2Y8xB2Dbd07A+qFdcc5hfwbZ59RqqIyKnicgEr1e9Dvu2aPd2FwH/ISKpIjIH227VSTY2U9eKyL7AD8IRJSJzgW8DpxtjmiK8JwCMMRuxjfW3iki6iByNfXB64/+ASuBZL0ZtgIicjC0p3GqM2dHHuZuAAyKQ9xTwPyIyTESGYhvXO0NEngEuFpHJnrO/JeSe2r39d4hItoiMBa4NObcv0rHtSVuANs/G/YZxdcfYOtbXgK97nyO5xiYgT0SGRHDJJ4GviMg5XhhTnohMN8Z0YJsEfiEiwwFEZD/vN+svz4aSjW262QKkicgt2BIiXjr/KSLDvOvVept7SmcPjDEtwM/Z/fstBgpF5D+8+zgXOAjrcGDvPNSnLu/4AhFJ8a5XDSwDfi4i+4hIioiMF5HO5/EZ4Nsiki8iudhOmd443bvHg7DNaNOxbZVv4Dk3EbkA2/44H/usPi4inTWoh7B5dKx37DARmed9/rKITPVqmXXYAkCf9nTGYRpjKrDtCYOxjee9MRFbAmvAtrs8YIxZ7u27GvsA1QLfwLZZdXIvth1yK7aE9rcwpZ2LrVqVio1raxCRh8I8N5RvYNu1arA9hn/EvoX3wisVfAX7dvwX9se8B9u7fXc/1/kJ1gHWish3w9D1I6wz/xAoBt7ztmGMWYJ10q9hq5T/9M7p1H0Vtu1uPbaE9wdsu2efGGPqsRn7GWyb03/Q92/eV1ofGWM+ivQaxpiPsS+L9Z6tRndPo4c0P8FW66/D9nYXYTvRwJbAyoF3vCafl7G1Eug7z4ayFFiC7TzbiO2pDq26zgE+EpEGbK/1ecaY5v50ezwK7C8iXzPG1GBLU9dh8+N/A6eFVHl/CZwlNubxV2Hoetb7XyMi73mfL8S+tFZj7f8cu5uZfuOl+QE2v/25D90XAb8zxnxijPm88w+4D/iGiOyPfbYvNMY0GGP+gM3Pvwi5l0XY5pB67LN/pLdvpKerDltV/wf9vPBl7xdzkkQgIn8EPjbGhFXSdQERmQyUYMOO/u3iY5MkcaaE+UVHRGZ61ZIUr7lgHnuWgJ1ERM7wmhFyse1mf0k6yyT/rsTEYYrIo2KHpZX0sv9Er3etyPu7pafjvuCMxIZrNGCruZcbY94PVFF4XIZtu1qHbd+5PFg5SZIER0yq5CJyPNYRPGGMmdLD/hOB7xpj+gqDSJIkSRKniUkJ0xjzOrYRPEmSJEm+sCSyDfNoEflA7JRoB/d/eJIkSZK4RaJmY3kPGGuMaRCRU7CdHRN7OlBEFgALADIzM2eMGTOGjIwMWlpa6OjoYPDgwezcuZMBAwYgIrS0tJCZmUlzs42uyMzMpLGxkQEDbMB+a2srgwYNoqnJhlBmZGTQ1NREeno6xhhaW1u70kxJSSE9PZ3m5mYGDhxIe3s7bW1tXftTU1MZMGBA1/62tjba29v32J+WlkZjYyODBw+mtbV1j/1paWmkpqaya9cuZ++pubkZEdnrnjo1u3hPO3fupKWlhWHDhkX0O7lyT3V1dQwYMCAmeS+R99SZfryfp3jcU1FR0VZjTOdIrLCJWViR2OmXXuqpDbOHYyuAw0Pivnpk0qRJZs2aNTHRl0iWL1/OiSeeGLSMqNCo/dFHH2XLli1873vfC1pKVGi0OejVDSAiq4wxh0d6XkKq5CIy0hvlgIgc4V23pr/zBg/uc9SXs8yYMSNoCVGjTbsxhrq6OlJSUti5c2fQcqJCm8070arbD7EKK3oKO4JhkohUisglIrJQRBZ6h5wFlIjIB9iQmvN6GMq2F1ofgFWrVgUtIWq0ad+0aRMtLS00NTWxfv36oOVEhTabd6JVtx9i0obpTQjQ1/77sEOZIiIlRWdcfWZmZtASokab9o8//rirjXD16tVMnTo1YEWRo83mnWjV7QenPVJ6enrQEqKioKAgaAlRo0372rW7Z/LatKnHCa6cR5vNO9Gq2w9OO8zOXi1tlJaWBi0hajRp7+jooKGhoet7a2srdXV1ASqKDk02D0Wrbj847TAHDgxscThfaH7zatJeVVVFS0tL1/edO3eybl1P8/u6jSabh6JVtx+cdpjt7f1O9eck9fX1QUuIGk3aS0tL96iFGGP46KO9ZnpzHk02D0Wrbj84vYxoW5vOSXFqavqNmHIWTdqnTp1KVlYW7733Hlu3buWkk05i1Ki9VvZwHk02D0Wrbj847TCTcZiJR5P2kSNHMnLkSOrr62ltbeXoo48OWlJUaLJ5KFp1+8HpKnkyDjPxaNW+a1ePk9erQKvNter2g9MOMzU13AUd3SIrK6v/gxxFq3atMbug1+ZadfvB6VzWOeBfG6NH97s8jLNo1a715Qp6ba5Vtx+cdpha4zDLysqClhA1WrW3trYGLSFqtNpcq24/OO0wtcZhjh/f09LVOtCqXWttBPTaXKtuPzjtMJNhRYlHq3atMbug1+ZadfvBaYep9SGora0NWkLUaNXe0dERtISo0Wpzrbr94LTDTMZhJh6t2rU234Bem2vV7QenHWYyDjPxaNWejMNMPFp1+8Fph6k1VCQnJydoCVGjVbvmOEytNteq2w9O57K0NKdHbvZKXl5e0BKiRqt2rS9X0Gtzrbr94LTD1FrN0jjFWCdatWuOw9Rqc626/eC0w8zIyAhaQlQUFhYGLSFqtGrXHIep1eZadfvBaYeptdRQVVUVtISo0apdawga6LW5Vt1+cNphan0IQpdN0IZW7ZrjMLXaXKtuPzjtMJNxmIlHq/ZkHGbi0arbD047zGQcZuLRql1rByHotblW3X5w2mEmw4oSj1btybCixKNVtx+cdphaH4Ls7OygJUSNVu0iErSEqNFqc626/eC0w9RazaqoqAhaQtRo1a51ZivQa3Otuv3gtMPUGoc5efLkoCVEjVbt6enpQUuIGq0216rbD047zJaWlqAlRIXmN69W7VpjdkGvzbXq9oPTDlNrbF1TU1PQEqJGq3ZjTNASokarzbXq9oPTDjMZh5l4tGpPxmEmHq26/eC0w0zGYSYerdq1dhCCXptr1e0Hpx2m1gkVhg8fHrSEqNGqXWsIGui1uVbdfnDaYWqNrdPcY6tVu9a8AnptrlW3H5x2mFp7ySsrK4OWEDVatWuOw9Rqc626/eC0w8zMzAxaQlRMmTIlaAlRo1W75tKOVptr1e0Hpx1mc3Nz0BKioqysLGgJUaNVu+Y4TK0216rbDzFxmCLyqIhsFpGSXvaLiPxKRMpF5EMROSwW13UVrfN4gl7tmuMwtdpcq24/xKqE+Rgwp4/9c4GJ3t8C4MFwEtVaJZ86dWrQEqJGq3bNVXKtNteq2w8xcZjGmNeBbX0cMg94wljeAXJEZFR/6TY2NsZCXsIpKioKWkLUaNWutYMQ9Npcq24/JKoNcz/g05Dvld62PtEahzlqVL/vAmdxXXtjI2zcCOXl9vvWrbB5M3R0pPLhh1BfD21tsGED1NWBhpq66zbvDa26/ZCoGXp7CpLrMSuLyAJstZ0RI0awfPlyJk+eTEVFBU1NTcyYMYNVq1YxfPhw0tPTqaysZMqUKZSVldHe3s7UqVMpKirq+jGrq6uZPn06xcXFpKamUlhYSElJCfn5+bS0tLB58+auNDMzMykoKKC0tJSCggLq6+upqanp2p+VlcXo0aMpKytj/Pjx1NTUUFtb27U/JyeHvLw8iouLyc7OpqqqioaGhq79eXl5ZGdnU1FR4ew9rV+/nurq6r3uad26dRQWFsb1nvbbbzrl5cUMGZLKG28UUltbwpYt+Ywe3cKsWZt56KEZbN26itTUTLZvL+Cuu0pZv76AFSu2MHRoC5deWs8VV6wiPz+L668fzdChZWzYMJ7zz69h1qxaliyZAawiOzuH/Pw8CgrWMXZsITU1wf5OxcXF1NbWxiTvJeJ36rynbdu2kZeXF/fnKR73FDXGmJj8AQVASS/7fg2cH/J9DTCqvzQLCwuNRl577bWgJURNIrR3dBizcaMxO3caU1ZmzKmnGjNqlDG5ucb8/vf2mBtvNObuu+33FSvstsZGe253li5dau68884er1NXZ0xNjf2+aJExv/iFMf/938bcdJPd9t3vGrPvvsZ85SvGXH+9MTt2GLNrlzHt7XG48V7Qml+06jbGGGClicLPJaqEuQi4UkSeBo4Edhhjqvs7adCgQXEXFg+mT58etISoiYf2tjZIS4NXXoGf/hTeew8GDICXXoIDDoBLLoH77oOxY6FzwM6Pf7x3On31AfbU6SMCoZOCf+1re593991w3XXw/vtWV2YmPPssLFwI06fDYYfBlVfChAkR3nQEaM0vWnX7IVZhRU8B/wQmiUiliFwiIgtFZKF3yGJgPVAO/Ab4Vjjpap0+qri4OGgJURMr7Vu2wOOPw9lnw/DhUFMDY8bAd74DxcVQXQ0zZkBuLpxxBhQU7HaW0eCn02fkSJg7F26+2Try//gPqKiAH/wARo2C1FSree5cuP9+24YaS7TmF626/RCTEqYx5vx+9hvgilhcSwOaJ4KIVrsx8MEHsHw5XH01/OEP8MYbcNpptvSYl2f/Cgtjq7eTWI8l33dfmD3b/gE0NMB//ZctFd96KzzwAJx+Orz7Lhx1lHWq0aI1v2jV7Qenl2XUukRFYby8QgKIRvuTT8L//I91GqedBi0t1mlefXUcBPZCvCMqsrJsafnss6G93f5VVcEVV9j/F11kPxcURJ621vyiVbcfnB4aqbVKXlLS44AnFYSjvaMDli2Diy+G1lYYPx7+9CdYuxbuvReCmMs3kXGYqamQnm7bXIuK4O23ISUF1q2DnTvhmWfsSyNctOYXrbr94LTD1Dp6Iz8/P2gJUdOf9hdfhIkT4XvfgyOPtM7zqKNs50iQM6wFuYb9hAm2M2v2bNi0CR58EPbfH2680ZY++0NrftGq2w9OO0yjIeq4BzSPOulJ+/vvwwUXwI4dtlf7qadsj/LChcGUJnvClbxywAHw2mu2LXfXLutAP//cRgj0htb8olW3H5x2mFpnoNm8eXPQEqImVPuGDXDOOXDKKbYEmZYGU6fCEUcEW5rsCdcmgjjwQLjnHjj0UOswFy6EWbPgn//c+1it+UWrbj843emTXAQt8cyYMYMNG2DQINs+OWMG/O534PpP4fIiaNOnw+rVNszqggvgzTdh2LDdPeta84tW3X5wuoSZXAQtsVRVwe23r+Lww2HFChsC9L3vue8swf1F0AYMgEsvhbIyG/d53XW29P7xx3rzi1bdfnDaYaakOC2vVzROS9fYaKvaAwZksmaNDQ/ShJY1fTqz9B132GaO446Dl17Sl19AZz73i9MeSWsveUE0wXgBsWgRXHONrYKXlcHVVxcwdGjQqiJH28xWgwfDDTfYUKzjjisA4NVXbdSBFjTl81jhtMPUukRFaWlp0BL6ZfNmOO88WzU8/XS7bdAgHdp7QmuPbU4ODBlSSmOjDUOaPXv31HWuozWv+MFph+lyQ35faHjzPvWUjRX84AM44YTd2zVo74kg4zD9UlBQwKBBNgD+a1+zca3vvBO0qv7Rmlf84HQucy1UJFzq6+uDltAjn34Kl19uS5W9DVt0VXt/uBKHGQ2dNk9NhWuvtZOR5OfbEKSsLBvK5SJa84ofnC5hal1ruqamJmgJe2AMPPKI7WQ48kg45pjej3VNe7hofbnC3jYfN872qn/6qY3dvO02O3bdNbTmFT847TCTcZixobUVXn/djj75/vftOOjecE17uGhtvoHebX7OOXas+uuv24k9XENrXvGD0w4zGYfpj08/tR07u3bBE0/AwQf3f44r2iPF9TjMvujL5vvtB0uX2hddYyN89FEChfWD1rziB6cdptb59rKysoKWwOuv2+r3YYfZdrBwcUF7NGiN2YX+bZ6WZh3n++/DiSfC008nRld/aM0rfnC600dbbF0no0ePDvT6mzbZWcMfewxOOimyc4PWHi1aX64Qvs2POQb+/nfbKfT++3DnncGO6deaV/zg9GtZaxxmWVlZINdtarLr0YwYYYPQI3WWEJx2v2idqAUis/n06bBype05F7Ht00GhNa/4wWmHqbUhf/z48Qm/5pYtNp7yz3+2o0WiXT8uCO2xQGttBCK3eV4e/Od/2gmLDzrIvhyDQGte8YPTDjMZVhQemzZZZ/nVr9q1dPw052kNFfkihRWFy/jxdnjliSfCq69+HltRYaA1r/jB6TZMrQ9BbW1twq5lDAwZArfcYnvE/ZJI7bGkQ9Mg7G74sfkll8DAgc3Mnj0KgOeee44zzzwzVtL6RGte8YPTJcxkHGbffPwxHH+8dZqxcJagN7ZOa/MN+Lf5f/5nBi+++BIAZ511VsKWv9WaV/zgtMNMxmH2zpo1dqKGSy+FWM6ypTW27osahxkO7e3tzJu3ez6+WbNP8SspLLTmFT847TC1hork5OTENf2WFjtf5R132OVdY0m8tccLzXGYfmze1tbWNfHI+vXruf32F6itHcrrr8dKXe9ozSt+cLoNU+sMNHl5eXFLe8cO22b5xht25u5YE0/t8UTryxWit3lLS0tXU8Qnn3zCmDFj+P73x/GlL82jtNQ218QTrXnFD06/lrVWs9atWxeXdLdvt7Oi//Of8XGWED/t8UZzHGY0Nm9ubu5yltXV1YwZM6Zr3+zZcNlldq34eM6tqTWv+MFph5mRkRG0hKgoLCyMeZqtrXDWWXDqqXD00TFPvot4aE8EmuMwI7V5Y2Nj1/IQW7ZsYWQvb88tW+z8mvHqzNaaV/zgtMPUWmqoqqqKeZoPPGA7d+6+O+ZJ70E8tCcCrSFoEJnN6+vru6JHtm3bxtA+1hNZuNCO9jrnnPiMCNKaV/zgtMPU+hA0NDTENL2mJvjWt+CPf9y9NGu8iLX2RKE5DjNcm9fW1rLPPvsAsGPHDnJzc/s95+c/h2nToK7Ol8Qe0ZpX/OC0w0zGYcJf/2pXFkxLS8xyt1pj677ocZg1NTVdDrKhoaHLcfZHWtruWsmf/xy1xB7Rmlf84LTD/HePwywuhosvhvvuS9ysNFpj67R2EEL/Nt+0aVNX1buxsTGqgkR9vZ2E+G9/i0pij2jNK35w2mH+O4cVdXTA+efDL39pF8VKFFpDRb6oYUWfffZZV6dOc3Nz1GuBFxTAc8/BhRdCrBZ71JpX/OC0R9L6EGRnZ/tOIyUFXnnFTtWWSGKhPQgkyIkhfdKbzTdu3Ni1MmNLS4vvSIBjjoGf/hRKSmDyZF9JAXrzih+cLmFqrWZVVFT4Ov/3v4cf/CDxzhL8aw8KrTNbQc82Ly8v73KWra2tMQubuvhiOPvs2JQyteYVPzjtMLXGYU728fr+9FO45hqYNy+GgiLAj/YgSe9rZTfH6W7zjz/+mIkTJwJ7Dn2MFdu32ynhVqzwl47WvOIHpx1mS0tL0BKiIto3rzF2uq5vf9uuxRMEWksNWmN2YU+bf/jhh12OqL29PS7NUrm5tm38wgttyFq0aM0rfnDaYWqNrWuKMhd2dNjRPDfeGGNBERCt9qAxxgQtIWo6bb5y5UqmTZsG2LwfzwlFzjsPDjkEHnoo+jS05hU/xOQXEZE5IrJGRMpF5IYe9s8XkS0iUuT9XRpOuv9OcZjl5fDaa7BggY2dCwqtsXXa4zDffvttZs6cSVZWFh0dHQnpxHrkEbjqKluziQatecUPvh2miKQC9wNzgYOA80XkoB4O/aMxZrr390g4af+7xGG2t9tp2lxYc1prbJ3WDkKAP/3pTxxzzDGMHDWKurq6hPX4Z2fb+MzjjotuJJDWvOKHWJQwjwDKjTHrjTEtwNNATLostE6oMHz48IiOf/hhW6q86qo4CYqASLW7gtYQtKVLl3LXXXcxYcJEqquqEh4elZsLkybZqIxI0ZpX/BALh7kf8GnI90pvW3fOFJEPReQ5ERnTw/690BpbF2mP7fPPw733+lu8LFZo7W3WmFcWLVrEnDlzGDNmDGvXBrdk7Y9/DP/3f7BhQ2Tnac0rfohFa1lPObV7q8hfgKeMMbtEZCHwODCrx8REFgALAIYNG8by5cuZPHkyFRUVNDU1MWPGDFatWsXw4cNJT0+nsrKSKVOmUFZWRnt7O1OnTqWoqIhRo+yiUNXV1UyfPp3i4mJSU1MpLCykpKSE/Px8Wlpa2Lx5c1eamZmZFBQUUFpaSkFBAfX19dTU1HTtz8rKYvT/Wbx5AAAgAElEQVTo0ZSVlTF+/Hhqamqora3t2p+Tk0NeXh5vvvkmgwYNoqqqioaGhq79eXl5ZGdnU1FR0XVPjY1NPPfcDN57bxWrVwd/TytXrqSysnKve1q3bh2FhYVh3VOif6ctW7bQ1NREfX19RL9TkPe0fPlynnzySWbNms3555/Pv/71r5jkvWjuqbR0FXfcMZyKinQ2bgz/nioqKsjJyYn78xSP3ylqjDG+/oCjgaUh328Ebuzj+FRgRzhpT5s2zWhky5YtYR332WfGzJxpTEtLnAVFQLjaXWLp0qXmnnvuCVpG2DzxxBMGMHPmzDHGuGPz2lpjNmwI/3hXdEcDsNJE4e9iUQlcAUwUkXEikg6cBywKPUBERoV8/ToQ1jiD5ubmGMhLPGVl4VWvbr3VrifuUlNtuNpdQ0sc5sMPP8yFF17ImWeeyZIlSwB3bP6nP8H8+eH3mruiO5H4dpjGmDbgSmAp1hE+Y4z5SERuF5Gve4d9W0Q+EpEPgG8D8/1e12XCmcdz9Wp44QW46aYECIoArXOQGgVxmPfeey+XXXYZF154Ec8991zXdldsfuGFdpb2xYvDO94V3YkkJhF/xpjFwOJu224J+XwjtqoeEdHOzBI04bSRtLbC//6v7aV0CV/tOwHiegfET37yE2666SYWLlzIgw8+uMc+V2yelmYn57jxRjjllP6nFHRFdyJxoF+2dxobG4OWEBVFRUV97l+/HiZMgHPPTZCgCOhPu6u4PIz2lltu4aabbuKaa67Zy1mCWzY/9VQbtRFO0IFLuhOF0w5TaxxmZ49ibyxcGPvZr2NFf9pdxdU4zOuvv54f/vCH3Hzzzdxzzz09HuOSzUVgzBjbvt5fjdsl3YnCaYf5ReTjj+GDD+zCVEm+2FxxxRX87Gc/44477uBHP/pR0HLCJj3dtmN6fVJJQnDaYWrp+exOdXV1r/vuuw+++U1wdehzX9pdxrUOiPnzL+aBBx7gnnvu4aZ+evZctPlVV9k29r5wUXe8cXrG9UGDBgUtISqmT5/e675vfhN6WUbaCfrS7jIudfqcc845PPvsszz44IMsXLiw3+NdtPk558D3vw+bN0NvIyBd1B1vnC5hap0+qri4uMft//qXdZZBzKQeLr1pdx1XOn1OPfVUnn32WX73u9+F5SzBTZsPHAhr1/buLMFN3fHGaYeplZ46IDo6bJyb67G+rnae9IcLY8lPPPHLLF68mKeffpr58+eHfZ6rNk9LszWi3mYyclV3PHHaYWpdoqKwsHCvbX//O2RmwrHHBiAoAnrSroGgIyoOP3wm//jHcl544QXOjTBezFWbi0BtLTz+eM/7XdUdT5x2mFqr5CUlJXttu/9+25DuQEGoT3rSroEgq+STDjyQVatWsmTJEuZFsRiTyzb/9rdtR2VPix+4rDteON3p41JDfiTk5+fvte2BB0DDMs49addAEGvYG2PIHzOGqs8+45VXXmHWrB4n4OoXl21+7LE2LnPDBhg/fs99LuuOF047TA3jg3uie2nnnXdgyBDQkL9c6TyJlETnFWMMubm57NixgzfffJNjjjkm6rRctrkIvPxyz/tc1h0vnK6Sa43D3Lx58x7fb7sNtNReumvXQiLjMI0xDBgwgB07dvDuu+/6cpbgvs1ra+0SKt3fSa7rjgdOO8wvwiJodXXw1lswZ06AgiJA68JWiVoErXM1x/b2doqKipg5c6bvNF23+ZAh8I9/7L3mlOu644HTDvOLsAja3/5m24GyswMUFAFaF7ZKxCJoHR0dXaE0H330UdeSuH5x3eYicPrpdjrCUFzXHQ+cbsOM57rM8SR0Wro5c+DwwwMUEyFap9SLdxxme3t7V8dSWVkZEydOjFnaGmx++um24zIUDbpjjdMeSWsveUFBAQAtLfDmm3DAAcHqiYRO7dqIZxxma2trl7PcsGFDTJ0l6LD5CSfAM8/suU2D7ljjtMPUukRFaaldgWP5clA0SQ2wW7s24tVj29LS0vXirqysjIuT0GBzEVi0CEImilehO9Y47TAT1ZAfazofqhdfhCjimANFa6khHnGYzc3NXXnw888/Z7/9elo92j9abN7eDr/+9e7vWnTHEqcdpmtTdoVLfX09AG+/DV/+csBiIqRTuzZiHYe5c+fOrja6LVu2MCKOM6ZosfkJJ9gJZDofSy26Y4nTDrOtrS1oCVFRU1MDwD33QIw6UhNGp3ZtxPLlWl9fT1ZWFgDbt29n6NChMUu7J7TYfN997WxbGzbY71p0xxKne8k1x2Fu2waTJ7s7UXBvaI2ti1XzTW1tLbneynR1dXVkJyAeTJPNS0rsjOygS3escLqEqTkO88kn7boo2tAaWxeLOMytW7d2OcudO3cmxFmCLptXVNjOH9ClO1Y47TC1zreXlZXFihW64i876ayKasNvzO7nn3/OsGHDADtLViJn+9dk8+pq+MlP7GdNumOF0w4z6DkOo2X06NGsXKnTYY4ePTpoCVHh5+VaWVnZtQLirl27Ej4PqyabH3YYfPghtLbq0h0rnHaYWuMw16wpY/58OPjgoJVETpnrU8L3QrQTtVRUVDBmzBhgz5jLRKLJ5tnZMHasHVeuSXescLrTR2sc5rhx49WFE3Uyvvukh0qIpjaydu3arlnD29raAmsC0mbz55+H/feHrVt16Y4FTpcwtYYVPfNMDd/7XtAqokNrqEikYUWrV6/ucpbt7e2Btpdrs/m++8L69fp0xwKnHabWwPXa2lriHLoXN2pra4OWEBUdPa2h0AtFRUUc7LWXtLe3Bz7Jizabv/023HSTPt2xwGmHqTUO89NPZ6C1PVxrbF24zTcrVqzg0EMPJTU1tWtuy6DRZvPRo6GqSp/uWBB8bukDrXGYBxywikmTglYRHVpj68KJw3zrrbc44ogjyMnJobW11YmleUGfzUePhs8+06c7FjjtMLXGYZ5xRo7KkCKAnJycoCVERX8lxVdffZVjjz2W/DFj2LZtmzPOEvTZfMQIWyXXpjsWOO0wg1gJ0C/GwN1356F1fag8DUtb9kBfL9clS5Ywe/ZsJh14IJ9+8olTzhL02TwtDa64AvbdV5fuWOC0w0zEsgOxpq4Otm5dh9K5j1m3bl3QEqKitzjMF154gVNOOYUZMw7nY0fnb9Ro82OOgTff1KfbL047zESPuIgFVVVQV1cYtIyo6Qy10UZPcZhPP/00Z5xxBieccCIrV64IQFV4aLR5Xh6kpOjT7RenHabGZXY3b4bx46uClhE1VVU6tXcPQXvsscc4//zzOeXUU1m+/LWAVIWHRpvbnnJ9uv3idCOhxjjME06Ajo6GoGVETUODTu2hcZgPPfQQl19+Oeeccw5//OMfA1QVHhptfswxsGOHPt1+iUkJU0TmiMgaESkXkRt62D9QRP7o7f+XiBSEk67GOMzSUmhs1BufpjW2rjMO8xe/+AWXX3458+fPV+EsQafNL7wQ5s/Xp9svvh2miKQC9wNzgYOA80XkoG6HXQJsN8ZMAH4B/DSctHuKw2xqamLNmjUsWbKEhx56iOrqal/6Y83bb8PLL+uNT9MaW7dr1y5+/OMfc+2113LFFVfwu9/9LmhJYaPR5n/+M/z2t27prqur4+GHH+bFF1+kpKQkLktoxKJKfgRQboxZDyAiTwPzgNUhx8wDbvU+PwfcJyJi+lmIJS0tjYaGBjZu3Mj69etZv3499fX1pKWlsWvXLlJTU9m0aVPX1Fwu0NYGxugNt9AW4tJJU1MTt912G9/97ne5++67g5YTERptvn49bNnilu4dO3awdetWqqurWb16Ne3t7WRmZlJQUMCECRMYO3as79hR8bt4lIicBcwxxlzqfb8AONIYc2XIMSXeMZXe93XeMVv7Sjs/P98sWLCAAQMG9LiMakpKSkRjiJMkSfLvSXp6Om1tbQwcOJCZM2cya9asVcaYiIeXxKKE2VMUcHcvHM4x9kCRBcACgNzcXFJTU7tmk+mtE2jYsGFdc2dmZGTQ1NTUFWbS2tpKZmbmXvvT09Pp6Oigra2NQYMG0djYSEpKCunp6TQ3N/e7f+DAgbS1tdHe3t61PzU1FWPS2LmzjtzcIbS0tNDR0dG1Py0tjZSUFFpaWsjIyOhzf0+aE3FPDQ0NDBgwYI976izR96c5qHvasWMHra2t7KirY0x+fli/k0v3tH37djIzM33nvUTe09atzbS372LEiCFxfZ4iuafa2toeZzgTETo6OhARMjIyfPWNxMJhVgJjQr7nA93jDTqPqRSRNGAIsK2nxIwxDwMPA0ybNs3cdNNN1NbWsnHjRsrLy6moqKC5uZnU1FTa2tqYN28e06dPj8FtxIaGBvjss01MmhS/ZVnjyaZNm+K6pGw8WLZsGStXruSX997LmP33Z2NFhXOjefpCo82Ngc8/38SoUe7o/vTTT3n88ce7hsmmpqay//77M2HCBAoKChg6dKjvfBELh7kCmCgi44DPgPOA/+h2zCLgIuCfwFnAq/21X4KdAVtEyM3NJTc3t8sx1tXV8cknn7Bx40bGjRsXg1uIHc88Ax98UMEvf+lORoqEiooKdQ8v2Pbu119/neOPP55hw4axZcsWNU5To82LiqC8vIKzz3ZH97BhwzjssMMYNWoUY8eOJTc3N+Z5wLfDNMa0iciVwFIgFXjUGPORiNwOrDTGLAJ+C/yfiJRjS5bnhZN2b+2T++yzD1OmTGHKlCl+5cccO/y9KWgZUdPUpFO7MYbjjjuOd955h6OOOoqMjAyam5tVOE2NNn/pJWhra+Lss4NWspuMjAxOOeWUuF4jJoHrxpjFwOJu224J+dwMRGxajXGYGRmwaZPe+DSNMYGwOw7zyCOP5L333uOwww4jJSXFiQmC+0OjzZuaIDNTn26/OJ2TNM6HOXUqHHaYW/FpkaAxJhD2nKjl0EMPpaSkBKBromCX0Wjzzz+HwYP16faL00MjNS6zO3kyGDM8aBlRM3y4Tu3dp3c7+OCDWbNmDZMmTerqIHR1flWNNr/ySti5U59uvzhdwtTQ/tSdlhb49rfT8RneGhhBLDMbC3rKK4WFhaxfvx6wnUKuLqqn0eYjRsCIEfp0+8Vph9lTsLrrpKfDgAGVbN8etJLoqKysDFpCVPTmDMeNG8cnn3wC0OsAiKDRaPNDDoG1a/Xp9ovTDjMzMzNoCVGxffsUtM585WLkQTj0VUobM2ZM15wDAwcO7ArKdgVtNm9uhvp6mDlTl+5Y4LTDdC1jh8vMmWUo7K8CoKysLGgJUdHf3KkjR45k8+bNgH0RNzY2JkJWWGizeXU1jBoF5eW6dMcCpx2mVs45p50jjwxaRXRonIMUbBxmfwwbNoxt2+wAs8GDB8dlNpto0GbzAQPg8sv16Y4FTjtMvVXyqbz0UtAqomPq1KlBS4iKcDtOcnNz2bFjB2AHQHR+DhJtNs/Phxtu0Kc7FjjtMF2qNkVCWVkRzz8ftIroKCoqClpCVETSmbPPPvt0zXKek5NDTU1NvGSFhTab338//OIX+nTHAqcdpsY4TIADDxyFwlhkAKfmFo2ESGMsBw8e3PVCHjp0KJs2bYqHrLDQZvPXXrNhRdp0xwKnHaZWCgpg7Vo7fCyJu4ROYzZy5Mh/y0W9omHlSjg84pkkvxg47TA1rhoJsHVrNRs3gsYmWNeW/AiXaDsgBg4c2FWd32+//di4cWMsZYWFJps3Ndk2zAkTdOmOFU47zEGDBgUtISqmT5/Opk2wwt2lsHvFpblFI8HPaJkBAwZ0vZwLCgooLy+Playw0GTzzEx4801ISdGlO1Y47TA1TnsFUFxczDvvwK9+FbSSyCkuLg5aQlT4HcETOnRy4sSJlJaWxkJWWGiy+RNP0BUBokl3rHDaYWolNTWVmTNtW482XJ2goj9iMe9A6DIoBx10EB9++KHvNMO9rhaee86O9AFdumOF0w4zIyMjaAlRUVhYyEEHwSefQF1d0Goio7CwMGgJURGriIrQhfWmTZvGygS89TTZfOVKmDnTftakO1Y47TC1VslLSkpIS4MlS+xkHJronEdSG7GcVKNz0azs7GxmzpzJ22+/HbO0e0KLzb2Rpey/v/2vRXcscdphapz2CiA/Px+woRfe7GJq6NSujbS02E7tKiLs2LGDkaNGccwxx/Daa6/FNP1QtNh8+HDYuBE6Wz+06I4lTjtMv2umB0VnaaekBM48M2AxEeLi9GfhEI+8IiJUV1UxYcJEZs2axdKlS2N+DdBj8wce2F3KBD26Y4nTDlNrHGbnrDiHHw47dsCaNQELioDNoU+EIuI5EcTatWVMn34oc+bMYdGiRTFPX4PNd+6048dDI/006I41TjtMjYugwe5FrVJSYN48ePHFgAVFgMYFuWD3Imjx4v333+PYY49j3rx5PPPMMzFNW4PNly2DI4+E3Nzd2zTojjVOO0yNi6DBnotaXXUVzJ4doJgI0bggF+y5CFq8eOON15kzZw7nnnsuTzzxRMzS1WDzxYvh9NP33KZBd6xxehE015dH7Y3QaekOOgg2bbLhRfvsE6CoMNE6pV6i1n9asmQJZ511FhdddBHNzc0sWLDAd5oabH7ffdC91UOD7ljjtEfS2kteUFCwx/fvfAeefjoYLZHSXbsWEjmz1XPPPccFF1zIZZddxi9/+Uvf6blu8+JiO0NR95HKruuOB047TK1LVHQfVnf66fDCCwGJiZBEDgmMJYnusX3iicdZuHAh3/nOd7jzzjt9peW6zX/7255HrbmuOx44XSWPd0N+vOj+5p07Fy67DLZsgWHDgtEULlpLDbGOwwyHBx98kEGDBnHjjTfS1NTEbbfdFlU6Ltu8rQ3+9Cf429/23uey7njhtMPUumZI97Vi9tkHHnnE9pq7jivr3ERKUDG7P//5z8nMzOT222+nsbGRu+++O+I0XLb5smUwdiwcfPDe+1zWHS+cdpi9rTXtOj0teXDWWbB6NQwZAgEUhsIm6OUaoiXIl+uPfvQjBg0axM0330xjYyP3339/ROe7bPO5c+l1QT+XdccLp8s82uMwu7NgAcQh7jmmaI2tC7r55qabbuKee+7hgQce4OKLL47oXFdtvnq1nc4tL6/n/a7qjidOO8wvQhxmKFdeCf/7vwkWEyFaY+sSEYfZH9dccw0PPvggjz32GOeee27Y57lq83vvtTNu9YaruuOJw5VDvfPtZWVl9bj9zDPhuutsmIarK5T2pt11XInZXbhwIZmZmcyfP5+Ghgb++te/9nuOizbftg2efRY+/rj3Y1zUHW+cdphaV40cPXp0j9sHDIA//3n39Fgu0pt213Hp5XrRRRcxcOBAzj//fE444UT+8Y/lfR7vos2XLoXTTrOrQ/aGi7rjjRuv5V7QGodZVlbW674jj4SiIti+PYGCIqAv7S7j2kQt5513Hs8//zyvv/4PDj98Zp/Humjz88+HRx/t+xgXdccbpx1m0A350TJ+/Pg+9//mN/DrXydITIT0p91VXKyNnH766SxevJhVq1Zy4OTJvR7nms0XL7Z5tD+TuqY7ETjtML9IYUWhfP/78POfg4tRGVpDRVyN2Z07dy6vvPIKaz7+mPwxY3qMF3XJ5q2tcO21MGZM/8e6pDtROO0wXX0I+qO2trbP/ZMmwdln2wlZXaM/7a7SuQ6Pi8yaNYs333yTzyoryc3N3ctpumTz3/7Wrjt+8sn9H+uS7kThq9NHRPYF/ggUABXAOcaYvVrnRKQd6FyT8xNjzNfDSf+LFocZyp13gotrvGmNrXO9+eaYY47h3Xff5Ygjjuha0rdzhiWXbF5eDnfdtXsZir5wSXei8FvCvAF4xRgzEXjF+94TTcaY6d5fWM4SvnhxmKHss49d7+fmmxMgKAK0xta5EIfZHzNnzqSoqIiOjo49Vqd0xeYNDfCzn8Fhh4V3vCu6E4lfhzkPeNz7/Dhweh/HRoxLoSKRkJOTE9ZxY8bAY4+5tX55uNpdw5U4zP6YNm0aH330EWDzd0dHhxM2r66GwkKIZHi4C7oTjd9cNsIYUw3g/R/ey3EZIrJSRN4RkbCdahAz0MSCvN7GknVj8GC49Va4/npwZb23cLW7hqaX60EHHcTatWsBq9sFx/ODH8A3vgHZ2eGfozWv+KFfjyQiLwMje9gVSWVyf2NMlYgcALwqIsXGmHW9XG8BsABg2LBhLF++nMmTJ1NRUUFTUxMzZsxg1apVDB8+nPT0dCorK5kyZQplZWW0t7czdepUioqKGDVqFADV1dVMnz6d4uJiUlNTKSwspKSkhPz8fFpaWti8eXNXmpmZmRQUFFBaWkpBQQH19fXU1NR07c/KymL06NGUlZUxfvx4ampqqK2t7dqfk5NDXl4er7zyCieddBJVVVU0NDR07c/LyyM7O5uKioquexo/vonRo2ewZMkqCgqCv6e33nqLkSNH7nVP69ato7CwMKx7SvTvtGXLFhobG6mvr4/odwr6np5//nkuvPBCbrvtNr75zW8yfvx433kvmnv6299W8eGHw3nooXSWLw//nioqKpg7d27cn6d4/E7RIn6mxRKRNcCJxphqERkFLDfGTOrnnMeAl4wxz/WX/rRp08wHH3wQtb6gqKqqingURGOj/d99VutEE432oFm2bBklJSVce+21QUuJmMrKSo444giqq6vZtWtXIKsMbNlih+vOmhXZeRrzSicissoYc3ik5/mtki8CLvI+XwTstT6iiOSKyEDv81DgGGB1OIm7NnojXKqqqiI+55Zb7DKmQRONdhfQGoKWn5/Pk08+Cdie/kSPbvvTn+x6U5E6S9CbV/zg12HeCXxVRNYCX/W+IyKHi8gj3jGTgZUi8gHwGnCnMSYsh6n1IWhoaIj4nJtvtuPMX3klDoIiIBrtLuByHGZ/iAhbtmwB7MJiiYoOKSuzKwFEazqtecUPvnpVjDE1wF6LyBpjVgKXep/fBqJqNPgix2F2JzfXzsp+6aWwZg0Etf6b1tg61+Mw+2LGjBlkZ2ezfft2cnNzycrKoq6ujuxIemAipK0NLrwQbrsNJk6MLg2tecUPTsdifJHjMHtizhx4+eXgnCXoja3TEIfZG502z8nJoa6uDoB99tknriNpyspg3Di4/PLo09CaV/zgtMP8oocV9cT48fCrX8GLe7UGJwatoSKawoq6E2rz7Ozsrqpubm4uW7dujfn1Nm+GAw+Ep57yt86U1rziB6cdptaHwG9V6rDDYOFCm7ETTTyrgfFEwhnL5yjdbT548GCampoAG1q3adOmmF1r1y6YPdvOd+kXrXnFD047TK3VrIqKCl/nH3ssXHABfPOb0TfIR4tf7UGhdWYr6NnmGRkZXfl/5MiRVFZWxuRaN90EEybY5h+/aM0rfnDaYWa4ODtFGEzuY+7DcPnhD+1kw4n2A7HQHgRBxC/Git5snp6e3hVaN2bMGN8O6o034IUX7FyXsSiQa80rfnDaYba0tAQtISpi8eYdONCWBtauteFGiUJrqUFrzC70bfPOmY0Axo0b1zWkMlJaW23N5Z13YOjQqJLYC615xQ9OO0ytsXWd7U+xwBjbnvnWWzFLsk9iqT2R+BmxFjT92Tw1NbUrJrmwsJDVq8MKY+5i/Xo4+GCorYVhw6KWuRda84ofnHaY/05xmL0xZQo8/jicdRYk4oWuNbZOexxmf6SkpHQ5zYMPPphwhwzv2GEXM7v6ahvrG0u05hU/OO0w/93iMHtj7lzbppmIXnOtsXVaOwghfJt3zqGZlpbG9OnTWbFiRb/nXHGFHfZ4xRV+Ve6N1rziB6cdposLW4XD8OG9zXIXPZdeCjNm2Ab7eHYExUN7ItAaggaR2VxEaGlpIScnhyOOOII333yz12M7OuCnP4V7742Fyr3Rmlf84LTD1BpbF68e244OeOYZuOqq+M2fqbW3WWtegchtLiJs27aN/fLzOe6443j11Vf3Ouaee2w1fL/9IF7jP7TmFT847TC19pLHKmauOwMGwHPPwfvv24chHk4zXtrjjeY4zGhsLiJUfvopkw48kNmzZ7NkyZKufb/6Fdx3n52YOp5ozSt+cNphZmZmBi0hKqZMmRK3tIcMsaM0Ro2Kj8OMp/Z4orm048fmH5eWcvjhMznllFN44YUXuOCCG7j6auG112D//WMosge05hU/OO0wEz03YKwoKyuLa/pDhsCNN8Lq1XaN81g6znhrjxea4zD92nzFinc58cQvc8YZZ/D73/+UU087k7FjYySuD7TmFT847TC1kqh5PPfbD/7+d1iwAGJ1Sa1zkGqOw4yFzVtad8+g+NCDcerl6YbWvOIHpx2m1iq5nzVDIiE31zrM9evhyitjk2aitMcazVVyPzY3Bq6/vpm33/oVXz35dNasWUN+fn4M1fWO1rziB6cdZmPnQjfKKCoqSti1srPhr3+Fb3/bhhv5DV1NpPZYorWDEPzZfNs2KC/PYPPmDpb97XkKCwtjqKxvtOYVPzjtMLXGYXausJcoMjJg8mR49ln40pdsiTNaEq09VmiOw4zG5jU1th17n33g+edh2LDEh1VpzSt+cNphJomM886zU8IdfTQsWxa0miTxoqgIZs607daKw09V4rTD1NrzWV1dHch1RWxb5rPP2qm8oiEo7X7R3AERic0rKuCrX4Wf/ATuuit+QenhoDWv+MFphzko6EW6o2T69OmBXv/44+3Y8w8/tAtdRbK4X9Dao0Vzp084Nm9rg3ffhYICW8I899z46+oPrXnFD047TK3TRxUXFwctAYDCQkhNte2a69aFd44r2iNFc6dPfzbfutXOkP6jH9le8f32S5CwftCaV/zgtMPUiisdEBkZ8Oijdj7N884LL8DdFe2RonkseV82/+AD2145Y4bt3HHpNrXmFT847TC1LlGRyNCO/hCBb30L/vEP20lw/fXw2We9H++S9kjQGlEBPdt8xw6oqoKRI+F//9fOOuSaf9KaV/zgtMPUWiUvKSkJWsJeDH9JByoAAA/1SURBVBpkS5iDB8Ohh8Ijj/Rc4nRRezhorpJ3t/lLL9mJo//8Zxgxwk4A7CJa84ofnF74W2tDfqJGWkTKgAFw661w5plwySW2jfP44/c8xlXt/aF1DXvY0+bf+pYNCXviCfjylwMUFQZa84ofnC5hah0f7HppZ+pUuxjW8cfDQw/ZuRM7o3Jc194bWvMKWJu/8oot8c+fb6MbXHeWoDev+MFph6k1DnNzItaS8EmK98t/5SuwaJFd0vfll3Vo7wmtcZjvvQf33LOZ73zHjt454gjbfKIBrXnFD047zOQiaPFnwgR49VX47ndh8WKrva9OIVfRuAjaP/5h2ycPOWQGq1bFbvnbRKEpn8cKpx1mchG0xJCSYsOO7rkHVqxYxUknwSmn2NKPFrQsgrZ+PVx0EfzhD3ad8PJyOOSQVWhsrteWz2OB0w4zJcVpeb2idVo6gMGDM3nvPeswTzsNHnwwaEXh4XocZnMzXH65rXKPGwennmrDhAYN0ptftOr2g9MeSWsveUFBQdASoqagoICBA+2Y9LVr4ZxzYMMGO8Tyn/+M3+JrfnE1DnP1atvrPXAgHHQQrFljIxWGDNl9jNb8olW3H5x2mFqXqCgtLQ1aQtSEah88GPLyYN99bc/6hRfC9OluVtVd67F9+mkbhTB7NpSW2gEEV11l7dkdrflFq24/OO0wNTbkg+43b0/ahwyxI4TWrLHtnGPGwMqVdiq5lSsTr7EnXIjDLCuDF16wn0tL7cqen3xi//eF1vyiVbcfnHaYWkNF6uvrg5YQNX1pT0mxJaZhw2DsWDjgADj7bDvOOej1sIKMw1y0CGbNguOOg48+sttuu80OEAinpUBrftGq2w9OO0yta03X1NQELSFqwtU+bJid8XvdOrjjDsjPh7/8xfb83nknlJQktr0zkS/XDRvs+O7bbrPfy8rgsstsafLmmyNPT2t+0arbD047zGQcZuKJVHtKip16bNAgO7Ht//wPVFbC174Gb75pg7GXLLG9xPEkns03bW12Igyw0QNHHWXbcQ891G777nft/JTRStCaX7Tq9oMvhykiZ4vIRyLSISKH93HcHBFZIyLlInJDuOkn4zATjx/tGRnWed53n403PPZY62h+/GMYPtyGKf3rX7bk+dlnsS2BxjIOs7Ow+rOf2SaHESPs+u9gh5JWV8Pvfgdf/3psrqc1v2jV7Qe/JcwS4P8Br/d2gIikAvcDc4GDgPNF5KBwEtc6315WVlbQEqImVtpF7N/UqXa5jA0b7DjpYcNg0ybb2z5yJMyda5fUAPj88+idaDQxu8ZARwds3w633w6nn27bZm+5xe5vbYV58+zY7t/+1m7bf//dw0pjhdb8olW3H3x1LRpjSqHfoOEjgHJjzHrv2KeBecDq/tJ3NbauP0aPHh20hKiJl/a8PDjrrN3fN2+2pcz33rNOFOCkk2zJdPRoOPFEePhh60zXr4dRo+zfiSdaJ9fcbNezSUuzji8lJZWWFlt97pzG7pNPbBtrdbX9u/pq20zwwx/u3vbII3DyydDUBN/4Bvz85zawHGwbbSLQml+06vZDImIx9gM+DfleCRzZ28EisgBYADBs2DCWL1/O5MmTqaiooKmpiRkzZrBq1SqGDx9Oeno6lZWVTJkyhbKyMtrb25k6dSpFRUVdS4BWV1czffp0iouLSU1NpbCwkJKSEvLz82lpaWHz5s1daWZmZlJQUEBpaSkFBQXU19dTU1PTtT8rK4vRo0dTVlbG+PHjqampoba2tmt/Tk4OeXl5vPLKK5x00klUVVXR0NDQtT8vL4/s7GwqKiqcvae33nqLkSNH7nVP69ato7CwMKb3VF5exKRJ9p6WL6/m9dens2JFMTt3pjJsWCHLl5fQ3JzP9u0tfPbZZsrLZyCyivLyTO6+u4Dx40vZsKGAyy7bwo4djYwcWc+MGasYMSKLu+4azSOPlFFUNJ7Ro2vIza1l27YZ1NWt4r/+K4f8/Dx27lzHtGmFfPxxFSefvPueGhsT+zstW7aMyZMnxyTvxeN36u2eKioqmDt3btyfp3jcU7RIf+EYIvIyMLKHXTcbY170jlkOfNcYs1dUnoicDZxsjLnU+34BcIQx5qr+xB1yyCHmww8/7PcmXOPTTz9lzJgxQcuICo3aly1bRklJCddee23QUqJCo81Br24AEVlljOm136U3+m2NMcZ8xRgzpYe/F8O8RiUQatV8oCqcE5NhRYlHq3atMbug1+ZadfshEWFFK4CJIjJORNKB84BF4Zyo9SGora0NWkLUaNXe0dERtISo0Wpzrbr94Des6AwRqQSOBv4qIku97aNFZDGAMaYNuBJYCpQCzxhjPgon/WQcZuLRql3rMFrQa3Otuv3gy2EaY543xuQbYwYaY0YYY072tlcZY04JOW6xMabQGDPeGHNHuOkn4zATj1btWubD7AmtNteq2w9Oj/TRGoeZk5MTtISo0apd69ypoNfmWnX7welc5sIMNNGQ19McXkrQql3ryxX02lyrbj847TC1VrPWrVsXtISo0apd64J5oNfmWnX7wWmHmZGREbSEqCgsLAxaQtRo1a51VBjotblW3X5w2mFqLTVUVYUVZuokWrVrDUEDvTbXqtsPTjtMrQ9BQ0ND0BKiRqt2zXGYWm2uVbcfnHaYyTjMxKNVezIOM/Fo1e0Hpx1mMg4z8WjVrrWDEPTaXKtuPzjtMJNhRYlHq/ZkWFHi0arbD047TK0PQXZ2dtASokar9n7mZHUarTbXqtsPTjtMrdWsioqKoCVEjVbtWme2Ar0216rbD047TK1xmJMnTw5aQtRo1Z6enh60hKjRanOtuv3gtMNsaWkJWkJUaH7zatWuNWYX9Npcq24/OO0wtcbWNTU1BS0harRq72/lAJfRanOtuv3gtMNMxmEmHq3ak3GYiUerbj847TCTcZiJR6t2rR2EoNfmWnX7wWmHqXVCheHDhwctIWq0atcaggZ6ba5Vtx+cdphaY+s099hq1a41r4Bem2vV7QenHabWXvLKysqgJUSNVu2a4zC12lyrbj847TAzMzODlhAVU6ZMCVpC1GjVrrm0o9XmWnX7wWmH2dzcHLSEqCgrKwtaQtRo1a45DlOrzbXq9oPTDlMrWufxBL3aNcdharW5Vt1+cNphaq2ST506NWgJUaNVu+YquVaba9XtB6cdZmNjY9ASoqKoqChoCVGjVbvWDkLQa3Otuv3gtMPUGoc5atSooCVEjVbtmuMwtdpcq24/OO0wkyRJksQlnHaYWns+q6urg5YQNVq1a+6A0Gpzrbr94LTDHDRoUNASomL69OlBS4gardo1d/potblW3X5w2mFqnT6quLg4aAlRo1W75k4frTbXqtsPTjtMrWjugNCqXfNYcq0216rbD04vy6h1iYrCwsKgJUSNJu0fffQRRUVFfPLJJ7S0tPDkk08yfvx4jjrqqKClRYQmm4eiVbcfnC5haq2Sl5SUBC0hajRpr62tpby8vKs6Xl5ernJCCE02D0Wrbj847TC1NuTn5+cHLSFqNGkvLCzcY1b+tLQ0Dj744AAVRYcmm4eiVbcfnHaYWscHa+6A0KR96NChpKXtblXKzMxk3LhxASqKDk02D0Wrbj847TC1xmFu3rw5aAlRo0m7iJCTk9P1PT09XWW7tyabh6JVtx+cdpjJRdASjzbtoWtja60iarN5J1p1+8GXwxSRs0XkIxHpEJHD+ziuQkSKRaRIRFaGm35yEbTEo037xIkTGTx4MCkpKSrbL0GfzTvRqtsPfsOKSoD/B/w6jGO/bIzZGkniKSlOF4B7Reu0dKBPe25uLgMGDKCtrY2xY8cGLScqtNm8E626/eDLYRpjSiF+QcNae8kLCgqClhA12rSLCPvuuy/btm1L5pcEo1W3HxJVhDPAMhFZJSILwj1J6xIVpaWlQUuIGo3aDznkELKysoKWETUabQ56dfuh3xKmiLwMjOxh183GmBfDvM4xxpgqERkO/F1EPjbGvN7L9RYAnU51l4hojI4dCkTU/OAQWrUPvfTSSzXqBsU2R6dugEnRnNSvwzTGfCWahLulUeX93ywizwNHAD06TGPMw8DDACKy0hjTa2eSq2jVDXq1a9UNerVr1Q1WezTnxb1KLiKDRSS78zNwErazKEmSJElU4Tes6AwRqQSOBv4qIku97aNFZLF32AjgTRH5AHgX+Ksx5m9+rpskSZIkQeC3l/x54PketlcBp3if1wPTorzEw9GrCxStukGvdq26Qa92rbohSu2idbx2kiRJkiQanZHhSZIkSRIAzjjMeA+zjCcRaJ8jImtEpFxEbkikxt4QkX1F5O8istb7n9vLce2ezYtEZFGidYbo6NOGIjJQRP7o7f+XiBQkXuXehKF7vohsCbHxpUHo7I6IPCoim3sL7xPLr7z7+lBEDku0xt4IQ/uJIrIjxOa39JuoMcaJP2AyNjZqOXB4H8dVAEOD1hupdiAVWAccAKQDHwAHOaD9LuAG7/MNwE97Oa7BAa392hD4FvCQ9/k84I9KdM8H7gtaaw/ajwcOA0p62X8KsAQQ4CjgX0FrjkD7icBLkaTpTAnTGFNqjFkTtI5oCFP7EUC5MWa9MaYFeBqYF391/TIPeNz7/DhweoBa+iMcG4bez3PAbAl+wR9Xf/t+MXaAybY+DpkHPGEs7wA5IjIqMer6JgztEeOMw4yAqIZZOsB+wKch3yu9bUEzwhhTDeD9H97LcRkislJE3hGRoJxqODbsOsYY0wbsAPISoq53wv3tz/Sqtc+JyJjESPONq/k6XI4WkQ9EZImI9DvdVUIXQUv0MMtYEgPtPZVyEhKi0Jf2CJLZ37P7AcCrIlJsjFkXG4VhE44NA7NzH4Sj6S/AU8aYXSKyEFtKnhV3Zf5x0d7h8h4w1hjTICKnAC8AE/s6IaEO0yR4mGUsiYH2SiC01JAPVPlMMyz60i4im0RklDGm2qtK9TiNdojd14vIcuBQbLtcIgnHhp3HVIpIGjCEGFfLoqBf3caYmpCvvwF+mgBdsSCwfO0XY0xdyOfFIvKAiAw1fUxDqapKrnyY5QpgooiME5F0bIdEYL3NISwCLvI+XwTsVVoWkVwRGeh9HgocA6xOmMLdhGPD0Ps5C3jVeC38AdKv7m7tfl8HtEwFtAi40OstPwrY0dnE4zoiMrKzfVtEjsD6w5o+Twq6Jyukx+oM7NtqF7AJWOptHw0s9j4fgO1h/AD4CFsdVqHd7O5RLMOWzFzRnge8Aqz1/u/rbT8ceMT7/CWg2LN7MXBJgHr3siFwO/B173MG8CxQjh2Ke0DQNg5T90+8PP0B8BpwYNCaPV1PAdVAq5fHLwEWAgu9/QLc791XMX1EuDio/coQm78DfKm/NJMjfZIkSZIkTFRVyZMkSZIkSJIOM0mSJEnCJOkw/387dSAAAAAAIMjfeoUBCiKASZgAkzABJmECTMIEmIQJMAVgOg+zltrUxQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x119cd7978>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure(figsize=(5, 5))\n",
    "ax = plt.axes()\n",
    "ax.grid(color='gray', linestyle='--', linewidth=0.5)\n",
    "\n",
    "plt.xlim((-1.5, 1.5))\n",
    "plt.ylim((-1.5, 1.5))\n",
    "\n",
    "# The U matrix in 2D space\n",
    "ax.arrow(0, 0, *U[:, 0], head_width=0.05, head_length=0.01)\n",
    "ax.arrow(0, 0, *U[:, 1], head_width=0.05, head_length=0.01)\n",
    "\n",
    "# x-axis\n",
    "ax.arrow(0, 0, 1.4, 0, head_width=0.05, head_length=0.1, linewidth=1, color='gray')\n",
    "ax.arrow(0, 0, -1.4, 0, head_width=0.05, head_length=0.1, linewidth=1, color='gray')\n",
    "\n",
    "# y-axis\n",
    "ax.arrow(0, 0, 0, 1.4, head_width=0.05, head_length=0.1, linewidth=0.7, color='gray')\n",
    "ax.arrow(0, 0, 0, -1.4, head_width=0.05, head_length=0.1, linewidth=0.7, color='gray')\n",
    "\n",
    "# Unit circle\n",
    "circ = plt.Circle((0, 0), radius=1, edgecolor='b', facecolor='None', linestyle='--')\n",
    "ax.add_patch(circ)\n",
    "\n",
    "plt.title(\"Visualizing Orthogonal Matrices as Rotated Axes\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Determinant"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The determinant of a square matrix $A \\in \\mathbb{R}^{n x n}$ is a function $\\mathbb{R}^{n x n} \\rightarrow \\mathbb{R}$, denoted as det A or $|A|$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 56,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 1.35624003, -0.07201012,  1.0035329 ],\n",
       "       [ 0.36163603, -0.64511975,  0.36139561],\n",
       "       [ 1.53803657, -0.03582604,  1.56464366]])"
      ]
     },
     "execution_count": 56,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 57,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-0.36796431478202435"
      ]
     },
     "execution_count": 57,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.linalg.det(A)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Eigenvalues and Eigenvectors"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Given a square matrix A, we say that $\\lambda \\in \\mathbb{C}$ is an eigenvalue of A and $x \\in \\mathbb{C}^{n}$ is the corresponding eigenvector if\n",
    "\n",
    "$Ax = \\lambda x, x \\ne 0$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Intuitively, this definition means that multiplying A by the vector x results in a new vector\n",
    "that points in the same direction as x, but scaled by a factor $\\lambda$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Can be re-written as \n",
    "\n",
    "$(\\lambda I - A)x = 0, x \\ne 0$\n",
    "\n",
    "$|\\lambda I - A| = 0$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 58,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 0.70710678,  0.70710678],\n",
       "       [-0.70710678,  0.70710678]])"
      ]
     },
     "execution_count": 58,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "U"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 59,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[ 0.70710678+0.70710678j  0.70710678-0.70710678j]\n",
      "[[ 0.70710678+0.j          0.70710678-0.j        ]\n",
      " [ 0.00000000+0.70710678j  0.00000000-0.70710678j]]\n"
     ]
    }
   ],
   "source": [
    "# Compute the eigenvalues and right eigenvectors of a square array.\n",
    "eigenvalues, eigenvectors = np.linalg.eig(U)\n",
    "print(eigenvalues)\n",
    "print(eigenvectors)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 60,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 0.70710678+0.70710678j,  0.70710678-0.70710678j])"
      ]
     },
     "execution_count": 60,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Compute the eigenvalues of a general matrix.\n",
    "np.linalg.eigvals(U)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Properties:**\n",
    "\n",
    "  1. The trace of A is equal to the sum of its eigenvalues\n",
    "  2. The determinant of A is equal to the product of its eigenvalues\n",
    "  3. The rank of A is equal to the number of non-zero eigenvalues of A."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 61,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1.41421356237\n",
      "(1.41421356237+0j)\n"
     ]
    }
   ],
   "source": [
    "print(np.trace(U))\n",
    "print(sum(eigenvalues))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 62,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1.0\n",
      "(1+0j)\n"
     ]
    }
   ],
   "source": [
    "print(np.linalg.det(U))\n",
    "print(np.product(eigenvalues))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 63,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2\n",
      "2\n"
     ]
    }
   ],
   "source": [
    "print(np.linalg.matrix_rank(U))\n",
    "print(sum(eigenvalues != 0))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Eigenvalues and Eigenvectors of Symmetric Matrices"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Two remarkable properties come about when we look at the eigenvalues and eigenvectors\n",
    "of a symmetric matrix $A \\in \\mathbb{S}^n$.\n",
    "\n",
    "1. All eigenvalues of A are real\n",
    "2. The eigenvectors of A are orthonormal"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 64,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[-5.23949021,  0.52289515, -0.13262482],\n",
       "       [ 0.52289515,  0.18352155, -1.63045634],\n",
       "       [-0.13262482, -1.63045634,  2.95578809]])"
      ]
     },
     "execution_count": 64,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A = np.random.randn(n, n)\n",
    "A = A + A.T # Use the fact that A + A.T is symmetric\n",
    "A"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 65,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[-5.28951678 -0.53329967  3.72263588]\n",
      "[[-0.99537851  0.08816279 -0.03806489]\n",
      " [ 0.0959833   0.90114531 -0.42275802]\n",
      " [ 0.00296953  0.42445784  0.90544283]]\n"
     ]
    }
   ],
   "source": [
    "eig_vals, eig_vectors = np.linalg.eig(A)\n",
    "print(eig_vals)\n",
    "print(eig_vectors)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can check for orthogonality of the matrix by computing the dot product of every column with every other column. Here is a nifty way to do it"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 66,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[  1.00000000e+00,   1.02656300e-16,  -3.33934269e-17],\n",
       "       [  1.02656300e-16,   1.00000000e+00,  -2.77555756e-16],\n",
       "       [ -3.33934269e-17,  -2.77555756e-16,   1.00000000e+00]])"
      ]
     },
     "execution_count": 66,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "orth_check = np.matmul(eig_vectors.T, eig_vectors)\n",
    "orth_check"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 67,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 1.,  0., -0.],\n",
       "       [ 0.,  1., -0.],\n",
       "       [-0., -0.,  1.]])"
      ]
     },
     "execution_count": 67,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Round the really small ones, at a significance level of 13 decimal digits\n",
    "np.around(orth_check, decimals=13)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Typically, the matrix of eigenvectors is represented using the letter U"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 68,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[-0.99537851,  0.08816279, -0.03806489],\n",
       "       [ 0.0959833 ,  0.90114531, -0.42275802],\n",
       "       [ 0.00296953,  0.42445784,  0.90544283]])"
      ]
     },
     "execution_count": 68,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "U = eig_vectors\n",
    "U"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can therefore represent A as $A = U \\Lambda U^T$, remembering from above that the inverse of an orthogonal matrix is just its transpose"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The Gradient"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The gradient of f w.r.t. $A \\in \\mathbb{R}^{m x n}$ is the matrix of partial derivatives, of dimensions $m x n$.\n",
    "\n",
    "The gradient of f w.r.t. $x \\in \\mathbb{R}^{n}$ is the matrix of partial derivatives, of dimensions $n$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 69,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sympy import symbols, diff"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 70,
   "metadata": {},
   "outputs": [],
   "source": [
    "x1, x2, x3 = symbols('x1 x2 x3', real=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 71,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2*x1"
      ]
     },
     "execution_count": 71,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "diff(x1**2 + x2**3 + x3**4, x1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 72,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3*x2**2"
      ]
     },
     "execution_count": 72,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "diff(x1**2 + x2**3 + x3**4, x2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 73,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "4*x3**3"
      ]
     },
     "execution_count": 73,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "diff(x1**2 + x2**3 + x3**4, x3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It is very important to remember that the gradient of a function is only defined if the function\n",
    "is real-valued, that is, if it returns a scalar value."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 74,
   "metadata": {},
   "outputs": [],
   "source": [
    "def gradient(f, x):\n",
    "    \"\"\"Calculates the gradient of a function f w.r.t x\"\"\"\n",
    "    \n",
    "    syms = [x1, x2, x3]\n",
    "    sub_dict = dict(zip(syms, x))\n",
    "    return np.array([diff(f, s).subs(sub_dict) for s in syms])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 75,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([2, 4, 6], dtype=object)"
      ]
     },
     "execution_count": 75,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "inp_vector = np.array([1, 2, 3])\n",
    "grads = gradient(x1**2 + x2**2 + x3**2, inp_vector)\n",
    "grads"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As you can see, this is \"symply\" (eh, eh, get it?) the vector `2x`, which is the derivative of `x**2` (input)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 76,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([2, 4, 6])"
      ]
     },
     "execution_count": 76,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "2*inp_vector"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### The Least Squares Equation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If we have $Ax = b$, where $A \\in \\mathbb{R}^{m x n}$ and $b \\in \\mathbb{R}^{m}$ s.t. $b \\notin R(A)$, then we want to find a vector $x$ such that $Ax$ is as close to b, as measured by the $l_2$ norm. It can be derived that the equation to solve this is:\n",
    "\n",
    "$$x = (A^T A)^{-1} A^Tb$$"
   ]
  }
 ],
 "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.6.4"
  },
  "varInspector": {
   "cols": {
    "lenName": 16,
    "lenType": 16,
    "lenVar": 40
   },
   "kernels_config": {
    "python": {
     "delete_cmd_postfix": "",
     "delete_cmd_prefix": "del ",
     "library": "var_list.py",
     "varRefreshCmd": "print(var_dic_list())"
    },
    "r": {
     "delete_cmd_postfix": ") ",
     "delete_cmd_prefix": "rm(",
     "library": "var_list.r",
     "varRefreshCmd": "cat(var_dic_list()) "
    }
   },
   "types_to_exclude": [
    "module",
    "function",
    "builtin_function_or_method",
    "instance",
    "_Feature"
   ],
   "window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}