slides/lecture15-midterm-review.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "i_f5u2x9nn6I",
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "<left><img width=25% src=\"img/cornell_tech2.svg\"></left>\n",
    "\n",
    "# Mid-Semester Review\n",
    "\n",
    "### Applied Machine Learning\n",
    "\n",
    "__Volodymyr Kuleshov__<br>Cornell Tech"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Announcements\n",
    "\n",
    "* The prelim is next __Wednesday__, in class, closed-book.\n",
    "    * Please make sure to have a pen or pencil\n",
    "    * Calculator is not needed (you can leave derivations unsimplified), but you can bring one if you want.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Part 1: Review of Supervised Learning\n",
    "\n",
    "We start with an overview of the supervised learning algorithms seen in class."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Supervised Machine Learning\n",
    "\n",
    "To apply supervised learning, we define a dataset and a learning algorithm.\n",
    "\n",
    "$$ \\underbrace{\\text{Dataset}}_\\text{Features, Attributes, Targets} + \\underbrace{\\text{Learning Algorithm}}_\\text{Model Class + Objective + Optimizer } \\to \\text{Predictive Model} $$\n",
    "\n",
    "The output is a predictive model that maps inputs to targets. For instance, it can predict targets on new inputs."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Linear Regression\n",
    "\n",
    "In linear regression, we fit a model\n",
    "$$ f_\\theta(x) := \\theta^\\top \\phi(x) $$\n",
    "that is linear in $\\theta$. \n",
    "\n",
    "The features $\\phi(x) : \\mathbb{R} \\to \\mathbb{R}^p$ may be non-linear in $x$ (e.g., polynomial features), allowing us to fit complex functions."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "We define the least squares objective for the model as:\n",
    "$$ J(\\theta) = \\frac{1}{2} \\sum_{i=1}^n (y^{(i)} - \\theta^\\top x^{(i)})^2 = \\frac{1}{2} (X \\theta - y)^\\top  (X \\theta - y) $$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "We can set the gradient to zero to obtain the *normal equations*:\n",
    "$$ (X^\\top X ) \\theta = X^\\top y. $$\n",
    "\n",
    "Hence, the value $\\theta^*$ that minimizes this objective is given by:\n",
    "$$ \\theta^* = (X^\\top X)^{-1} X^\\top y.$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "<center><img width=100% src=\"img/algorithms3.png\"></center>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Overfitting\n",
    "\n",
    "Overfitting is one of the most common failure modes of machine learning.\n",
    "* A very expressive model (a high degree polynomial) fits the training dataset perfectly.\n",
    "* The model also makes highly incorrect predictions outside the training set, and doesn't generalize."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "We can measure overfitting and underfitting by estimating accuracy on held out data and comparing it to the training data.\n",
    "* If training performance is high but holdout performance is low, we are overfitting.\n",
    "* If training performance is low and holdout performance is low, we are underfitting."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "We will see many ways of dealing with overfitting, but here are some ideas:\n",
    "* Use a simpler model family (linear models vs. neural nets)\n",
    "* Keep the same model, but collect more training data\n",
    "* Modify the training process to penalize overly complex models."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Regularization\n",
    "\n",
    "The idea of regularization is to penalize complex models that may overfit the data."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "Regularized least squares optimizes the following objective (__Ridge__).\n",
    "$$ J(\\theta) = \\frac{1}{2n} \\sum_{i=1}^n \\left( y^{(i)} - \\theta^\\top \\phi(x^{(i)}) \\right)^2 + \\frac{\\lambda}{2} \\cdot ||\\theta||_2^2. $$\n",
    "If we use the L1 norm, we have the __LASSO__."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "<center><img width=100% src=\"img/algorithms4.png\"></center>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Regression vs. Classification\n",
    "\n",
    "Consider a training dataset ${D} = \\{(x^{(1)}, y^{(1)}), (x^{(2)}, y^{(2)}), \\ldots, (x^{(n)}, y^{(n)})\\}$.\n",
    "\n",
    "We distinguish between two types of supervised learning problems depending on the targets $y^{(i)}$. \n",
    "\n",
    "1. __Regression__: The target variable $y \\in {Y}$ is continuous:  ${Y} \\subseteq \\mathbb{R}$.\n",
    "2. __Classification__: The target variable $y$ is discrete and takes on one of $K$ possible values:  ${Y} = \\{y_1, y_2, \\ldots y_K\\}$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "<center><img width=100% src=\"img/algorithms5.png\"></center>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Parametric vs. Non-Parametric Models\n",
    "\n",
    "Nearest neighbors is an example of a *non-parametric* model.\n",
    "* A parametric model has a finite set of parameters $\\theta \\in \\Theta$ whose dimensionality is constant with data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "* In a non-parametric model, the function $f$ uses the entire training dataset to make predictions, and the complexity of the model increases with dataset size."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "* Non-parametric models have the advantage of not losing any information at training time. \n",
    "* However, they are also computationally less tractable and may easily overfit the training set."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Probabilistic Models\n",
    "\n",
    "A probabilistic discriminative model outputs a vector of class probabilities\n",
    "\n",
    "$$\n",
    "\\underbrace{\\left[ \n",
    "\\begin{array}{c}\n",
    "x_1 \\\\\n",
    "x_2 \\\\\n",
    "x_3 \\\\\n",
    "x_4 \\\\\n",
    "\\end{array}\n",
    "\\right]}_\\text{input}\n",
    "{\\to}\\text{model $P_\\theta(y|x)$}{\\to}\n",
    "\\underbrace{\\left[ \n",
    "\\begin{array}{c}\n",
    "0.1 \\\\\n",
    "0.7 \\\\\n",
    "0.2 \\\\\n",
    "\\end{array}\n",
    "\\right]}_\\text{output}\n",
    "$$\n",
    "\n",
    "The model takes $x$ to be a fixed input."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "For example, a logistic (softmax) model outputs a probability distribution \n",
    "\n",
    "$$\n",
    "P_\\theta(y| x) \n",
    "= \\left[ \\begin{array}{c}\n",
    "P_\\theta(y=0| x) \\\\\n",
    "P_\\theta(y=1| x) \\\\\n",
    "\\end{array} \\right]\n",
    "= \\left[ \\begin{array}{c}\n",
    "1-\\sigma(\\theta^\\top x) \\\\\n",
    "\\sigma(\\theta^\\top x) \\\\\n",
    "\\end{array} \\right]\n",
    "$$\n",
    "\n",
    "where $\\theta^\\top x$ is a linear model and\n",
    "$$\\sigma(z) = \\frac{1}{1 + \\exp(-z)}$$\n",
    "is the *sigmoid* or *logistic* function."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Generative Models: Intuition\n",
    "\n",
    "Another approach to classification is to use *generative* models.\n",
    "\n",
    "* A generative approach first builds a model of $x$ for each class:\n",
    "$$ P_\\theta(x | y=k) \\; \\text{for each class $k$}.$$\n",
    "$P_\\theta(x | y=k)$ *scores* each $x$ according to how well it matches class $k$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "* A class probability $P_\\theta(y=k)$ encoding our prior beliefs.\n",
    "These are often just the % of each class in the data."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "In the context of spam classification, we would fit two models on a corpus of emails $x$ with spam/non-spam labels $y$:\n",
    "\n",
    "\\begin{align*}\n",
    "P_\\theta(x|y=\\text{0}) && \\text{and} && P_\\theta(x|y=\\text{1})\n",
    "\\end{align*}\n",
    "\n",
    "* $P_\\theta(x | y=1)$ *scores* each $x$ based on how much it looks like spam.\n",
    "* $P_\\theta(x | y=0)$ *scores* each $x$ based on whether it looks like non-spam.\n",
    "\n",
    "We also choose a prior $P(y)$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "In the context of spam classification, given a new $x'$, we would compare the probabilities of both models:\n",
    "\n",
    "\\begin{align*}\n",
    "P_\\theta(x'|y=\\text{0})P_\\theta(y=\\text{0}) && \\text{vs.} && P_\\theta(x'|y=\\text{1})P_\\theta(y=\\text{1})\n",
    "\\end{align*}\n",
    "\n",
    "We output the class that's more likely to have generated $x'$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "<center><img width=100% src=\"img/algorithms6.png\"></center>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "# Generative vs. Discriminative Approaches\n",
    "\n",
    "What are the pros and cons of generative and discriminative methods?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "* If we only care about prediction, we don't need a model of $P(x)$. It's simpler to only model $P(y|x)$ (what we care about).\n",
    "    * In practice, discriminative models are often be more accurate."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "* If we care about other tasks (generation, dealing with missing values, etc.) or if we know the true model is generative, we want to use the generative approach."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Text Classification & Bag of Words\n",
    "\n",
    "Perhaps the most widely used approach for representing text documents is called \"bag of words\".\n",
    "\n",
    "We start by defining a vocabulary $V$ containing all the possible words we are interested in, e.g.:\n",
    "$$ V = \\{\\text{church}, \\text{doctor}, \\text{fervently}, \\text{purple}, \\text{slow}, ...\\} $$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "A bag of words representation of a document $x$ is a function $\\phi(x) \\to \\{0,1\\}^{|V|}$ that outputs a feature vector\n",
    "$$\n",
    "\\phi(x) = \\left( \n",
    "\\begin{array}{c}\n",
    "0 \\\\\n",
    "1 \\\\\n",
    "0 \\\\\n",
    "\\vdots \\\\\n",
    "0 \\\\\n",
    "\\vdots \\\\\n",
    "\\end{array}\n",
    "\\right)\n",
    "\\begin{array}{l}\n",
    "\\;\\text{church} \\\\\n",
    "\\;\\text{doctor} \\\\\n",
    "\\;\\text{fervently} \\\\\n",
    "\\\\\n",
    "\\;\\text{purple} \\\\\n",
    "\\\\\n",
    "\\end{array}\n",
    "$$\n",
    "of dimension $V$. The $j$-th component $\\phi(x)_j$ equals $1$ if $x$ contains the $j$-th word in $V$ and $0$ otherwise."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Bernoulli Naive Bayes Model\n",
    "\n",
    "The *Bernoulli Naive Bayes* model $P_\\theta(x,y)$ is defined for *binary data* $x \\in \\{0,1\\}^d$ (e.g., bag-of-words documents)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "The $\\theta$ contains prior parameters $\\vec\\phi = (\\phi_1,...,\\phi_K)$ and $K$ sets of per-class parameters $\\psi_k = (\\psi_{1k},...,\\psi_{dk})$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "The probability of the data $x' = [0,1,0,...]$ for each class equals\n",
    "\\begin{align*}\n",
    "P_\\theta(x=x'|y=k) \n",
    "& = \\prod_{j=1}^d P_\\theta(x_j = x_j' \\mid y=k) %= \\prod_{j=1}^d \\psi_{kj}^{x_j} \\cdot (1-\\psi_{kj})^{(1-x_j)}\\\\\n",
    "\\end{align*}\n",
    "\n",
    "where each $P_\\theta(x_j \\mid y=k)$ \n",
    "is Bernoulli: \n",
    "\n",
    "\\begin{align*}\n",
    "P_\\theta(x_j = 1\\mid y=k) & = \\psi_{jk} \\\\\n",
    "P_\\theta(x_j = 0\\mid y=k) & = 1- \\psi_{jk}\n",
    "\\end{align*}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "The probability over $y$ is Categorical:\n",
    "$P_\\theta(y=k) = \\phi_k$. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "For example, suppose that $x' = [0, 1, 0, 0, 1]$. Then we have:\n",
    "\n",
    "\\begin{align*}\n",
    "P_\\theta(x=x'|y=k) \n",
    "& = \\prod_{j=1}^d P_\\theta(x_j = x_j' \\mid y=k) \\\\\n",
    "& = (1-\\psi_{1k}) \\cdot \\psi_{2k} \\cdot (1-\\psi_{3k}) \\cdot (1-\\psi_{4k}) \\cdot \\psi_{5k} \\\\\n",
    "\\end{align*}\n",
    "\n",
    "i.e., we multiply the probabilities of each word being present or absent."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Maximum Likelihood Learning\n",
    "\n",
    "We can learn a generative model $P_\\theta(x, y)$ by maximizing the *likelihood*:\n",
    "\n",
    "$$ \\max_\\theta \\frac{1}{n}\\sum_{i=1}^n \\log P_\\theta({x}^{(i)}, y^{(i)}). $$\n",
    "\n",
    "This seeks to find parameters $\\theta$ such that the model assigns high probability to the training data."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "Given a dataset $\\mathcal{D} = \\{(x^{(i)}, y^{(i)})\\mid i=1,2,\\ldots,n\\}$, we want to optimize the log-likelihood $\\ell(\\theta) = \\log L(\\theta)$:\n",
    "\\begin{align*}\n",
    "\\ell & = \\sum_{i=1}^n \\log P_\\theta(x^{(i)}, y^{(i)}) = \\sum_{i=1}^n \\sum_{j=1}^d \\log P_\\theta(x^{(i)}_j | y^{(i)}) + \\sum_{i=1}^n \\log P_\\theta(y^{(i)}) \\\\\n",
    "& = \\sum_{k=1}^K \\sum_{j=1}^d \\underbrace{\\sum_{i :y^{(i)} =k} \\log P(x^{(i)}_j | y^{(i)} ; \\psi_{jk})}_\\text{all the terms that involve $\\psi_{jk}$} + \\underbrace{\\sum_{i=1}^n \\log P(y^{(i)} ; \\vec \\phi)}_\\text{all the terms that involve $\\vec \\phi$}.\n",
    "\\end{align*}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "Let's first consider the optimization over $\\vec \\phi = (\\phi_1, \\phi_2, \\ldots, \\phi_K)$. \n",
    "$$ \\max_{\\vec \\phi} \\sum_{i=1}^n  \\log P_\\theta(y=y^{(i)} ; \\vec \\phi). $$\n",
    "* We have $n$ datapoints, each having one of $K$ classes\n",
    "* We want to learn the most likely class probabilities $\\phi_k$ that generated this data\n",
    "\n",
    "What is the maximum likelihood $\\phi$ in this case?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "Intuitively, the maximum likelihood class probabilities $\\phi$ should just be the class proportions that we see in the data."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "Let's calculate this formally. Our objective $J(\\vec \\phi)$ equals\n",
    "\\begin{align*}\n",
    "J(\\vec\\phi) & = \\sum_{i=1}^n  \\log P_\\theta(y^{(i)} ; \\vec \\phi) = \\sum_{i=1}^n  \\log \\phi_{y^{(i)}}\n",
    "\\end{align*}\n",
    "where $\\phi_k$ are the probabilities of each class."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "Normally, the parameters of the categorical distribution are class probabilities $\\phi_k$, that sum to one.\n",
    "\n",
    "However, let's assume for this derivation, that they don't have to sum to one. Anything that we will derive will still be true if they do sum to one, so we don't lose anything by making this assumption.\n",
    "\\begin{align*}\n",
    "J(\\vec\\phi) & = \\sum_{i=1}^n  \\log P_\\theta(y^{(i)} ; \\vec \\phi) = \\sum_{i=1}^n  \\log \\left( \\frac{\\phi_{y^{(i)}}}{\\sum_{k=1}^K \\phi_k}\\right)\n",
    "\\end{align*}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "We further expand this using properties of the log:\n",
    "\\begin{align*}\n",
    "J(\\vec\\phi) & = \\sum_{i=1}^n  \\log \\left( \\frac{\\phi_{y^{(i)}}}{\\sum_{k=1}^K \\phi_k}\\right) \\\\\n",
    "& = \\sum_{i=1}^n \\log \\phi_{y^{(i)}} - n \\cdot \\log \\sum_{k=1}^K \\phi_k \\\\ \n",
    "& = \\sum_{k=1}^K \\sum_{i : y^{(i)} = k} \\log \\phi_k - n \\cdot \\log \\sum_{k=1}^K \\phi_k\n",
    "\\end{align*}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "Let's take the partial derivative with respect to $\\phi_k$:\n",
    "\n",
    "\\begin{align*}\n",
    "\\frac{\\partial}{\\partial_{\\phi_k}} J(\\vec\\phi) \n",
    "& = \\frac{\\partial}{\\partial_{\\phi_k}} \\left (\\sum_{l=1}^K \\sum_{i : y^{(i)} = l} \\log \\phi_l - n \\cdot \\log \\sum_{l=1}^K \\phi_l \\right) \\\\\n",
    "& = \\frac{\\partial}{\\partial_{\\phi_k}} \\sum_{i : y^{(i)} = k} \\log \\phi_k - \\frac{\\partial}{\\partial_{\\phi_k}} n \\cdot \\log \\sum_{l=1}^K \\phi_l \\\\\n",
    "& = n_k \\cdot \\frac{\\partial}{\\partial_{\\phi_k}} \\log \\phi_k - n \\cdot \\frac{\\partial}{\\partial_{\\phi_k}}  \\log \\sum_{l=1}^K \\phi_l \\\\\n",
    "& =  \\frac{n_k}{\\phi_k} - \\frac{n}{\\sum_{l=1}^K \\phi_l}\n",
    "\\end{align*}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "Taking the derivative and setting it to zero, we obtain \n",
    "$$ \\frac{\\phi_k}{\\sum_l \\phi_l} = \\frac{n_k}{n}$$\n",
    "for each $k$, where $n_k = |\\{i : y^{(i)} = k\\}|$ is the number of training targets with class $k$.\n",
    "\n",
    "Thus, the optimal $\\phi_k$ is just the proportion of data points with class $k$ in the training set!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Max-Margin Classification\n",
    "\n",
    "Maximizing margins is a principle for choosing good decision boundaries."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "* Margin: distance from the decision boundary to closest datapoint.\n",
    "<center><img width=40% src=\"img/margin.png\"></center>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "* We want the margin to be as large as possible."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Support Vector Machine: Primal Form\n",
    "\n",
    "We saw that maximizing the margin of a linear model amounts to solving the following optimization problem.\n",
    "\n",
    "\\begin{align*}\n",
    "\\min_{\\theta,\\theta_0} \\frac{1}{2}||\\theta||^2 \\; &  \\\\\n",
    "\\text{subject to } \\; & y^{(i)}((x^{(i)})^\\top\\theta+\\theta_0)\\geq 1 \\; \\text{for all $i$}\n",
    "\\end{align*}\n",
    "\n",
    "We minimize the L2 norm among all correct decision boundaries."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Support Vector Machine: Dual Form\n",
    "\n",
    "The solution to this problem also happens to be given by the following *dual* problem:\n",
    "\\begin{align*}\n",
    "\\max_{\\lambda} & \\sum_{i=1}^n \\lambda_i - \\frac{1}{2} \\sum_{i=1}^n \\sum_{k=1}^n \\lambda_i \\lambda_k y^{(i)} y^{(k)} (x^{(i)})^\\top x^{(k)}  \\\\\n",
    "\\text{subject to } \\; & \\sum_{i=1}^n \\lambda_i y^{(i)} = 0\n",
    "\\end{align*}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "Also, the predictions at a new point $x'$ are given by \n",
    "$$y' = \\sum_{i=1}^n \\lambda_i^* y^{(i)} K(x^{(i)}, x').$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "* It's better to use the primal when we have $n > d$\n",
    "* Use the dual when we have $d > n$\n",
    "* The dual allows us to use kernels"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Kernelized SVM\n",
    "\n",
    "We can replace each dot product via a kernel function\n",
    "\n",
    "\\begin{align*}\n",
    "\\max_{\\lambda} & \\sum_{i=1}^n \\lambda_i - \\frac{1}{2} \\sum_{i=1}^n \\sum_{k=1}^n \\lambda_i \\lambda_k y^{(i)} y^{(k)} K(x^{(i)}, x^{(k)})  \\\\\n",
    "\\text{subject to } \\; & \\sum_{i=1}^n \\lambda_i y^{(i)} = 0 \\\\\n",
    "& C \\geq \\lambda_i \\geq 0 \\; \\text{for all $i$}\n",
    "\\end{align*}\n",
    "\n",
    "Kernels can compute the dot product more efficiently."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "How to choose features or kernels?\n",
    "* Most often you want to start with a kernel (don't worry about the features)\n",
    "* There are conditions that guarantee that a kernel corresponds to some dot product"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Are kernelized SVMs: \n",
    "\n",
    "* For classification or regression? \n",
    "* Discriminative or generative? \n",
    "* Parametric or non-parametric?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Consider an update rule for kernelized k-NN:\n",
    "$$y' = \\sum_{i=1}^n y^{(i)} K(x^{(i)}, x').$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "This is practically the same as kernelized SVM:\n",
    "$$y' = \\sum_{i=1}^n \\lambda_i^* y^{(i)} K(x^{(i)}, x').$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "You can view the SVM as learning an optimal weighting of data in kernelized k-NN (hence, it's non-parametric)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "<center><img width=100% src=\"img/algorithms8.png\"></center>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "<left><img width=25% src=\"img/cornell_tech2.svg\"></left>\n",
    "# Part 2: Unsupervised Learning\n",
    "\n",
    "Next, we review algorithms for unsupervised learning."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "# Unsupervised Learning\n",
    "\n",
    "We have a dataset *without* labels. Our goal is to learn something interesting about the structure of the data:\n",
    "* __Clusters__ hidden in the dataset.\n",
    "* A __low-dimensional representation__ of the data.\n",
    "* Recover the __probability density__ that generated the data."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "# Density Estimation\n",
    "\n",
    "The problem of density estimation is to approximate the data distribution $P_\\text{data}$ with the model $P$.\n",
    "$$ P \\approx P_\\text{data}. $$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "It's also a general learning task. We can solve many downstream tasks using a good model $P$:\n",
    "* Outlier and novelty detection\n",
    "* Generating new samples $x$\n",
    "* Visualizing and understanding the structure of $P_\\text{data}$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "Let's look at an example in the context of the 1D points we have seen earlier.\n",
    "\n",
    "We will fit a model of the form\n",
    "$$P(x) = \\sum_{i=1}^n K(x, x^{(i)})$$\n",
    "with a Gaussian kernel $K(x,z; \\delta) \\propto \\exp(-||x-z||^2/2\\delta^2)$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[ 1.62434536 -0.61175641 -0.52817175 -1.07296862  0.86540763 -2.3015387\n",
      "  6.74481176  4.2387931   5.3190391   4.75062962  6.46210794  2.93985929\n",
      "  4.6775828   4.61594565  6.13376944  3.90010873  4.82757179  4.12214158\n",
      "  5.04221375  5.58281521]\n"
     ]
    }
   ],
   "source": [
    "# https://scikit-learn.org/stable/auto_examples/neighbors/plot_kde_1d.html\n",
    "import numpy as np\n",
    "np.random.seed(1)\n",
    "\n",
    "N = 20 # number of points\n",
    "# concat samples from two Gaussians:\n",
    "X = np.concatenate((\n",
    "    np.random.normal(0, 1, int(0.3 * N)), \n",
    "    np.random.normal(5, 1, int(0.7 * N))\n",
    "))[:, np.newaxis]\n",
    "bins = np.linspace(-5, 10, 10) # locations of the bins\n",
    "\n",
    "# print out X\n",
    "print(X.flatten())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(-0.02, 0.32)"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAskAAAD4CAYAAAAejHvMAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nOzdd3zV1f3H8dfJDlkEEiCQsPceAVy4WkVFcaE46wBHrR2/1rpq1dpaq7bVDrUOcODAXal71YUywh4yw0hCICEhe957z++PExRTlCBJvvcm7+fjcR93fe+9nwv3m/u+53uGsdYiIiIiIiJfC/O6ABERERGRYKOQLCIiIiLSiEKyiIiIiEgjCskiIiIiIo0oJIuIiIiINBLhdQGNpaSk2N69e3tdhoiIiIi0cUuWLNltrU3d331BF5J79+5NVlaW12WIiIiISBtnjNn2bfepu4WIiIiISCMKySIiIiIijSgki4iIiIg0opAsIiIiItKIQrKIiIiISCMKySIiIiIijSgki4iIiIg0opAsIiIiItJIk0KyMeYkY8x6Y8wmY8yN+7n/amPMKmPMcmPMZ8aYofvcd1PD49YbYyY3Z/EiIiIiIi3hgCHZGBMOPACcDAwFzt83BDd41lo7wlo7GrgH+GvDY4cC5wHDgJOABxueT0REREQkaDWlJXkCsMlam22trQPmAqfvu4G1tmyfq3GAbbh8OjDXWltrrd0CbGp4PhERERGRoBXRhG16ADn7XM8FJjbeyBjzE+CXQBRw/D6PXdDosT3289grgSsBevbs2ZS6RURERERaTLMN3LPWPmCt7QfcANxykI99xFqbaa3NTE1Nba6SRERERES+l6aE5DwgY5/r6Q23fZu5wBnf87EiIiIiIp5rSkheDAwwxvQxxkThBuLN23cDY8yAfa5OATY2XJ4HnGeMiTbG9AEGAIsOvWwRERERkZZzwD7J1lqfMeZa4B0gHJhtrV1jjLkDyLLWzgOuNcb8EKgH9gCXNDx2jTHmBWAt4AN+Yq31t9B7ERERERFpFsZae+CtWlFmZqbNysryugwRERERaeOMMUustZn7u08r7omIiIiINKKQLCIiIiLSiEKyiIiIiEgjCskiIiIiIo0oJIuIiIiINKKQLCIiIiLSiEKyiIiIiEgjCskiIiIiIo0oJIuIiIiINKKQLCIiIiLSiEKyiIiIiEgjCskiIiIiIo0oJIuIiIiINKKQLCIiIiLSiEKyiIiIiEgjCskiIiIiIo0oJIuIiIiINKKQLCIiIiLSiEKyiIiIiEgjCskiIiIiIo0oJIuIiIiINKKQLCIiIiLSiEKyiIiIiEgjTQrJxpiTjDHrjTGbjDE37uf+Xxpj1hpjVhpjPjDG9NrnPr8xZnnDaV5zFi8iIiIi0hIiDrSBMSYceAA4AcgFFhtj5llr1+6z2TIg01pbZYz5MXAPML3hvmpr7ehmrltEREREpMU0pSV5ArDJWpttra0D5gKn77uBtfa/1tqqhqsLgPTmLVNEREREpPU0JST3AHL2uZ7bcNu3mQG8tc/1GGNMljFmgTHmjP09wBhzZcM2WYWFhU0oSURERESk5Rywu8XBMMZcBGQCx+xzcy9rbZ4xpi/woTFmlbV2876Ps9Y+AjwCkJmZaZuzJhERERGRg9WUluQ8IGOf6+kNt32DMeaHwG+Aqdba2r23W2vzGs6zgY+AMYdQr4iIiIhIi2tKSF4MDDDG9DHGRAHnAd+YpcIYMwZ4GBeQC/a5PdkYE91wOQU4Eth3wJ+IiIiISNA5YHcLa63PGHMt8A4QDsy21q4xxtwBZFlr5wH3AvHAi8YYgO3W2qnAEOBhY0wAF8j/1GhWDBERERGRoGOsDa4uwJmZmTYrK8vrMkRERESkjTPGLLHWZu7vPq24JyIiIiLSiEKyiIiIiEgjCskiIiIiIo0oJIuIiIiINKKQLCIiIiLSiEKyiIiIiEgjCskiIiIiIo0ccDERERERaXnWWgoratldXketz48/YEmKjSQ5LorOcVE0LNYlIq1EIVlERMQD9f4AX2wuYv7m3SzILmbTrnIq6/z73TYxJoLBaYmMyejIMYNSyezViagIHQwWaUkKySIiIq1o465ynlm4nXkrdlBcWUdkuGF0RkfOycygd+cOdEuKIToynDBjKK2up6iilk0FFXyZX8bs+Vt4+JNsEmIiOG1Ud6ZnZjAyPUmtzCItQCFZRESkFazOK+UfH27knTW7iIoI44ShXTlzdA+O6N+ZDlFN+zqurPUxf9Nu3l69k1eW5vLswu2M7dmRa4/vz3GDuigsizQjY631uoZvyMzMtFlZWV6XISIi0ix2V9Ry79vreT4rh8SYCC49ojeXHtmHTnFRh/S8ZTX1vLo0j0c+ySavpJqR6Un89tShjO/dqZkqF2n7jDFLrLWZ+71PIVlERKT5WWt5dVket89bQ1Wdn8uP6sO1x/cnMSayWV+n3h/g1WV5/PXdDewsq2HKyDRuO3UoXRJjmvV1RNqi7wrJ6m4hIiLSzEqq6rj51VW8uWon43snc9dZI+nfJb5FXisyPIxzMzM4bWR3Hv5kMw99tJnPNu7mttOGcuaYHuqCIfI9qSVZRESkGa3dUcaVc7LYVVbD/50wkKuO7kd4WOsF1ezCCq5/aSVZ2/Zw/OAu3DNtJCnx0a32+iKh5LtakjV/jIiISDN5Y2U+Zz/0OfX+AC9cdTjXHNu/VQMyQN/UeJ6/6nB+e+pQ5m/azSl/+5SF2UWtWoNIW6CQLCIi0gwen7+Fnzy7lKHdE/nPT49iTM9kz2oJDzPMOKoPr15zJHHREVzw2EIe/GgTgUBwHT0WCWYKySIiIofAWstf3l3P7/6zlsnDuvLMzIl0SQiOQXNDuycy79ojOWlYN+55ez0/fW4Z1d+yYImIfJMG7omIiHxPgYDl1nmreXrBdqZnZnDnmcOJCA+u9qeEmEj+ecEYRn6SxJ/eXkfunioe/VGmZr8QOYDg2pNFRERChLWW2/+zhqcXbOeqo/vyp7NHBF1A3ssYw1XH9OPhi8axYVcFpz8wn3U7y7wuSySoBefeLCIiEsSstfzp7XU89cU2rjy6LzeePDgkplo7cVg3Xvrx4VgL5/7rC7K2FntdkkjQUkgWERE5SH//YBMPf5zNxYf14qYQCch7DeuexEs/PpyU+GgumrWQ/64r8LokkaCkkCwiInIQnl24nfve38C0cen8buqwkArIe6Und+CFqw+nf5d4rngqi38vy/O6JJGg06SQbIw5yRiz3hizyRhz437u/6UxZq0xZqUx5gNjTK997rvEGLOx4XRJcxYvIiLSmj5aX8BvX1vNcYNS+dNZIwhr5TmQm1NKfDTPXXEYmb2T+cXzy3l24XavSxIJKgcMycaYcOAB4GRgKHC+MWZoo82WAZnW2pHAS8A9DY/tBNwGTAQmALcZY7ybOFJEROR7WrujjJ88s5RBXRP45wVjg3aQ3sFIiInkicsmcNygVG5+dRXPLNzmdUkiQaMpe/gEYJO1NttaWwfMBU7fdwNr7X+ttVUNVxcA6Q2XJwPvWWuLrbV7gPeAk5qndBERkdZRUFbD5U8sJjE2ktmXjicuuu3MoBoTGc6/Lh7H8YO78JtXVzNngYKyCDQtJPcAcva5nttw27eZAbx1MI81xlxpjMkyxmQVFhY2oSQREZHWUecLcM0zSymtrmfWJePpltT25heOjgjnoYvG8sMhXfjtv1cz54utXpck4rlmPVZkjLkIyATuPZjHWWsfsdZmWmszU1NTm7MkERGRQ3LnG2vJ2raHe6aNZGj3RK/LaTHREeE8cOFYfjikK799bQ1Pfr7V65JEPNWUkJwHZOxzPb3htm8wxvwQ+A0w1VpbezCPFRERCUYvL8nlyS+2ccWkPpw2qrvX5bS46IhwHrxwLCcM7cpt89bw/GIN5pP2qykheTEwwBjTxxgTBZwHzNt3A2PMGOBhXEDed8LFd4ATjTHJDQP2Tmy4TUREJKitzivl5ldXcVjfTtxw0mCvy2k1URFh/POCMRwzMJUbX1nFvBU7vC5JxBMHDMnWWh9wLS7cfgm8YK1dY4y5wxgztWGze4F44EVjzHJjzLyGxxYDv8cF7cXAHQ23iYiIBK3SqnqufnoJneKi2sxMFgcjOiKcf100jvG9O/HL55fz/tpdXpck0uqMtdbrGr4hMzPTZmVleV2GiIi0U9ZarnlmKe+t3cWLVx/OmJ7td+bS8pp6LnpsIV/uLOfxS8dzZP8Ur0sSaVbGmCXW2sz93de+fhqLiIgcwLOLtvPW6p38evKgdh2Qwc2j/OTlE+ibEsfMJ7PI2qqDwdJ+KCSLiIg0WL+znDv+s5ZJA1K4YlJfr8sJCh07RDFnxkS6JcVw2eOLWZ1X6nVJIq1CIVlERASorvPz0+eWkhATwV/OHRXSS043t9SEaJ6eOZHE2EgunrWQjbvKvS5JpMUpJIuIiAB/eGMtG3ZV8NdzR9Mloe0tGHKoenSM5ZmZEwkPC+OiWQvJKa468INEQphCsoiItHufbizkmYXbufLovhw9UItafZveKXE8PXMCNfUBLpq1kIKyGq9LEmkxCskiItKuVdT6uPHlVfRNjeOXJwz0upygN7hbIo9fNp7C8lounrWIkqo6r0sSaREKySIi0q7d/dY6dpRWc++0kcREhntdTkgY2zOZR3+UyZbdlVz6+GIqa31elyTS7BSSRUSk3fpicxFzFmzj8iP7MK5XJ6/LCSlH9k/hHxeMYVVeKVfOyaKm3u91SSLNSiFZRETapao6Hze8vJJenTtw3YmDvC4nJE0e1o17p41k/qYifvrcMnz+gNcliTQbhWQREWmX/vzOBrYXV3H32SOJjVI3i+/rrLHp/G7qMN5bu4vrX1pJIBBcK/mKfF8RXhcgIiLS2pZsK+bxz7dw8WG9OKxvZ6/LCXmXHNGb8pp6/vzuBhJiIrh96jCM0TzTEtoUkkVEpF2pqffz65dW0j0plhtOHux1OW3GT47rT1mNj0c+ySYxNpJfqQuLhDiFZBERaVfuf38j2YWVzJkxgfhofQ02F2MMN508mLLqev7x4SYSYiK48uh+Xpcl8r3pr4OIiLQbK3JKeOSTzZw3PoNJA7RoSHMzxnDnmSMor/XxxzfXkRATyfkTenpdlsj3opAsIiLtQq3Pz69fWkGXhBhunjLE63LarPAww33njqay1sfNr64iPjqC00Z197oskYOm2S1ERKRdeODDTWzYVcEfzxpOYkyk1+W0aVERYTx04TjG9+rE/z2/nP+uL/C6JJGDppAsIiJt3podpTz40WbOGtOD4wd39bqcdiE2KpzHLs1kcFoCP356CYu2FHtdkshBUUgWEZE2rd4f4PqXVtKxQxS3njbU63LalcSYSJ68bAI9OsYy44nFrM4r9bokkSZTSBYRkTbt4Y83s2ZHGX84YzgdO0R5XU670zk+mqdnTiQxNpIfzV7EpoIKr0sSaRKFZBERabM27Crn7x9sYsrINE4a3s3rctqttKRYnpk5kTBjuOixheQUV3ldksgBKSSLiEib5PMH+PVLK4mPieCOqcO8Lqfd650Sx5wZE6iq83HxrIUUlNd4XZLId1JIFhGRNmn2/C2syCnh9qnD6Bwf7XU5AgxJS+SJyydQUF7Lj2YtorSq3uuSRL6VQrKIiLQ52YUV/OXdDZw4tCunjUzzuhzZx9ieyTxycSbZhZVc+sQiKmt9Xpcksl9NCsnGmJOMMeuNMZuMMTfu5/6jjTFLjTE+Y8y0Rvf5jTHLG07zmqtwERGR/QkELNe/tJKYyHD+cMZwjDFelySNHDUghX9cMIaVuaVc9sRiquoUlCX4HDAkG2PCgQeAk4GhwPnGmMZz6GwHLgWe3c9TVFtrRzecph5ivSIiIt/pyS+2krVtD7eeOpQuiTFelyPfYvKwbtw/fTRZW4u57HEFZQk+TWlJngBsstZmW2vrgLnA6ftuYK3daq1dCQRaoEYREZEm2V5UxT1vr+fYQamcNbaH1+XIAZw2qjv3nzeGxVuLufyJxVTX+b0uSeQrTQnJPYCcfa7nNtzWVDHGmCxjzAJjzBn728AYc2XDNlmFhYUH8dQiIiJOIGC54eWVRIQZ7jprhLpZhIipo7pz3/TRLNqioCzBpTUG7vWy1mYCFwD3G2P6Nd7AWvuItTbTWpuZmpraCiWJiEhb89zi7XyRXcTNU4aQlhTrdTlyEE4f3YO/njuahVuKmPGkgrIEh6aE5DwgY5/r6Q23NYm1Nq/hPBv4CBhzEPWJiIgcUF5JNXe9uY4j+3fmvPEZB36ABJ0zxvTgL+eO4ovsImY+paAs3mtKSF4MDDDG9DHGRAHnAU2apcIYk2yMiW64nAIcCaz9vsWKiIg0Zq3lpldWEbCWP501Ut0sQtiZY9L5yzmj+HxzEVc8lUVNvYKyeOeAIdla6wOuBd4BvgResNauMcbcYYyZCmCMGW+MyQXOAR42xqxpePgQIMsYswL4L/Ana61CsoiINJuXluTyyYZCbjx5MBmdOnhdjhyis8am8+dpo5i/eTeXPb5Y8yiLZ4y11usaviEzM9NmZWV5XYaIiISAXWU1nPDXjxncLZG5Vx5GWJhakduKV5fl8qsXVjC2ZzKzLxtPYkyk1yVJG2SMWdIwdu5/aMU9EREJSdZafvPqKmp9Ae6eNlIBuY05c0w6/7xgLMtzSrj4sYWUVNV5XZK0MwrJIiISkuat2MH7Xxbw68mD6JMS53U50gJOGZHGvy4ax5f55Zz/6EKKKmq9LknaEYVkEREJOYXltdw+bw1jenbksiP7eF2OtKAfDu3Ko5dkkl1YwXmPLKCgrMbrkqSdUEgWEZGQYq3lln+vorLOz73TRhKubhZt3jEDU3nisgnklVQz/ZEF7Cip9rokaQcUkkVEJKTMW7GDd9bs4pcnDKR/lwSvy5FWcni/zsyZMYHd5bWc+/AX5BRXeV2StHEKySIiEjIKymu4bd4aRmd05IpJfb0uR1rZuF6deOaKiZTX+Dj34S/YXFjhdUnShikki4hISLDWcsurq6mq8/Pnc0apm0U7NTK9I3OvPIx6f4Bz//UFq/NKvS5J2iiFZBERCQnzVuzg3bW7+NUJA+nfJd7rcsRDQ9ISeeGqw4mOCOP8RxeQtbXY65KkDVJIFhGRoLe3m8WYnh2ZqW4WAvRNjefFHx9Banw0F81ayMcbCr0uSdoYhWQREQlqbtEQ183i3mnqZiFf69ExlheuPpy+KfHMfHIxb63K97okaUMUkkVEJKj9e3ke763dxXUnqpuF/K+U+Gieu/IwRqV35CfPLuWFrByvS5I2QiFZRESCVu6eKm799xoyeyUz4yh1s5D9S4qN5KkZEziyfwrXv7SS2Z9t8bokaQMUkkVEJCj5A5ZfvbACC9w3fbS6Wch36hAVwWOXZHLy8G7c8fpa7n9/A9Zar8uSEKaQLCIiQenRT7NZuKWY204bSkanDl6XIyEgOiKcf5w/hnPGpXP/+xv5/etfEggoKMv3E+F1ASIiIo2t2VHKX95dz8nDuzFtXLrX5UgIiQgP4+6zR5IQE8ns+Vsor6nnrrNGEBGudkE5OArJIiISVGrq/fxi7nKSO0TxxzNHYIy6WcjBCQsz/PbUISTFRnLf+xuoqPVx/3mjiY4I97o0CSH6WSUiIkHl7rfXsbGggnvPGUVyXJTX5UiIMsbw8x8O4NZTh/LW6p3MfDKLqjqf12VJCFFIFhGRoPHfdQU8Pn8rlx7Rm2MGpnpdjrQBlx/Vh3umjWT+pt1cPGsRpdX1XpckIUIhWUREgkJ+aTW/fGE5Q9ISufHkwV6XI23IuZkZPHDBWFbmlnD+IwvYXVHrdUkSAhSSRUTEcz5/gJ89t4w6X4AHLhhDTKT6jkrzOnlEGrMuGU/27gqmP/wF+aXVXpckQU4hWUREPHff+xtYvHUPfzxrBH1TtaqetIyjB6YyZ8ZECspqOedfX7C9qMrrkiSIKSSLiIinPtlQyIMfbea88RmcPrqH1+VIGze+dyeeveIwKmt9TPvX52zcVe51SRKkFJJFRMQzu8pq+L/nlzOwSwK3nTbM63KknRiRnsTzVx2OBaY/soDVeaVelyRBqEkh2RhzkjFmvTFmkzHmxv3cf7QxZqkxxmeMmdbovkuMMRsbTpc0V+EiIhLa6v0BfvrcMqrq/Dxw4Rhio9QPWVrPwK4JvHjV4cRGhnP+IwvI2lrsdUkSZA4Yko0x4cADwMnAUOB8Y8zQRpttBy4Fnm302E7AbcBEYAJwmzEm+dDLFhGRUHfnG1+yaEsxd501gv5dErwuR9qh3ilxvHj14aQmRHPxrEV8tnG31yVJEGlKS/IEYJO1NttaWwfMBU7fdwNr7VZr7Uog0Oixk4H3rLXF1to9wHvASc1Qt4iIhLCXl+TyxOdbmXFUH84Yo37I4p3uHWN5/qrD6dW5A5c/sZh31+z0uiQJEk0JyT2AnH2u5zbc1hRNeqwx5kpjTJYxJquwsLCJTy0iIqFoVW4pN726iiP6deYmzYcsQSA1IZq5Vx7G0O6J/PiZpby2PM/rkiQIBMXAPWvtI9baTGttZmqqVlgSEWmrdlfUctWcLFLjo/nH+WOICA+KryEROnaI4umZExnfO5lfPL+cZxdu97ok8VhT/jrlARn7XE9vuK0pDuWxIiLShtT6/Pz46SUUVdbx8MXj6Bwf7XVJIt8QHx3BE5dN4NiBqdz86ioe/STb65LEQ00JyYuBAcaYPsaYKOA8YF4Tn/8d4ERjTHLDgL0TG24TEZF2xFrL9S+tZPHWPfz5nFEM75HkdUki+xUTGc7DF2cyZUQad775Jfe9twFrrddliQciDrSBtdZnjLkWF27DgdnW2jXGmDuALGvtPGPMeOBVIBk4zRjzO2vtMGttsTHm97igDXCHtVZzrIiItDP3vb+R15bv4NeTB3HaqO5elyPynaIiwvj7+WPoEBXO3z7YSFWdj5tPGYIxxuvSpBUdMCQDWGvfBN5sdNut+1xejOtKsb/HzgZmH0KNIiISwl5eksvfP9jIuZnpXHNsP6/LEWmS8DDD3WePpENUOI9+uoXqej93TB1OWJiCcnvRpJAsIiLyfXy+aTc3vrKSI/t35s4zR6glTkJKWJjh9qnDiIkK5+GPs6muC3D32SM04LSdUEgWEZEWsSKnhCueyqJPShwPXjiOSAULCUHGGG48aTBxURH89b0N1NT7uW/6aKIi9Hlu6xSSRUSk2W3cVc4ljy+iU3wUc2ZMJCk20uuSRL43Yww/+8EAYiPDufPNL6mp9/PAhWOJidRS6m2ZfgaJiEizyimu4qJZC4kMD+PpGRPpmhjjdUkizeKKo/vy+zOG88G6AmY+mUVVnc/rkqQFKSSLNJG1lpp6P5W1Pnz+xiuwiwjArrIaLpq1kJr6AHNmTKBX5zivSxJpVhcf1os/nzOKzzfv5pLZiyivqfe6JGkh6m4hggvAeSXVrMsvZ8vuSrJ3V5K7p4rC8lqKKusoraqnrlEwjggzxEaGk5IQTWp8NKmJ0fTpHMeArvH0S3Wn2CgdipP2Y0dJNRc8uoDC8lqenjmRwd0SvS5JpEVMG5dOTGQYv5i7nAsfW8hTl0+gY4cor8uSZqaQLO2Szx9gWU4JC7OLWJ5TwvKcEnZX1H11f6e4KDKSY0lPjmV0Rkc6dogiOiKMqIgwIsIMtb4AtT4/lbV+iirrKCirYU1eKW+v3ok/4CadDzMwuFsi43olM65XMpm9k0lP7uDVWxZpUTnFVVzw2AL2VNYzZ8YExvZM9rokkRZ16sjuxESEc80zSznvkQXMmTGR1AStItmWmGBbRSYzM9NmZWV5XYa0QQVlNbz/ZQEfbyjg801FlNe6vmT9UuMYnZHM6IwkhvVIom9K3PduEaj1+dlWVMWmggrW5ZexdHsJy7bvobLOD0DflDiOGZTKMQNTOaxvZw36kDZhe1EV5z+6gLKaeubMmMjojI5elyTSaj7buJsrnsoirWMMz8ycSFpSrNclyUEwxiyx1mbu9z6FZGnLCspreHv1Tl5fmc/ircVYCz06xnL0wBSOHpDKEf1SSOrQsqPu/QHLhl3lfLG5iI83FLIgu4haX4AOUeEcP7gLU0akceygLuqaISFp465yLp61iBqfn6dnTNRy09IuLd5azGWPLyY5LpJnZx5GRicdNQwVCsnSrtT5Arz/5S6eX5zDpxsLCVgY0CWeKSPTOGVEGgO6xHu6oEFNvZ8F2UW8u3YX76zeSVFl3VeB+fTRPTh2UKrmk5WQsHhrMTOeWExURDhPXT6Bod3VB1narxU5Jfxo9iJiI8N55oqJ9EuN97okaQKFZGkXNu4q57lFOfx7eR7FlXWkJcVw9th0po7uzsCuCV6Xt18+f4BFW4p5fVX+V4E5JT6KM0b3YFpmugY+SdB6e3U+P5u7nPSOsTx5+QS1nIkAX+aXcfGshQAavBoiFJKlzbLW8tmm3Tz66RY+2VBIZLjhxKHdOCcznUkDUgkPC50lcOv9AT5eX8hLS3L5YN0u6v2W4T0SmTY2ndNH9yA5TiOnxXvWWh79NJu73lrH6IyOzLpkPJ302RT5yqaCCi58bAG1vgBPXT6Bkenqox/MFJKlzanzBZi3YgePfZrNup3lpCZEc+kRvTl/Qs828YVdXFnHvOV5vLgklzU7yogKD+OEoV05d3wGR/VPCanw3yTl5bB6NZSVfX3b4MHQq5d3Ncn/qKn3c9Mrq3h1WR5TRqTx53NGhVZfemthyxbYtMldBggLgxNO8Lau76uyElatgtLSr28bMAD69vWuJgHcYNYLHltAaVU9j182nszenbwuSb6FQrK0GSVVdTyzcDtPfr6VgvJaBnVNYOakPkwd3Z3oiBD6sj4Ia3eU8eKSHP69LI89VfWkJcUwbVw654zLoGfnED3EHQjAnj3QuTMUFMB5530dWva65ho45xzYsQNuvx3GjIHMTHfysE95e7WjpJqrn17CytxSrjtxID85rr+nffubbPdu9zkzBu6/H1577Zv3R0XBO++4yw89BPn5MHo0HHUUdOnS+vV+F2uhuNi9nz174Oyz/3e/mTkTLrwQCgvhN79x72XcOBg/3v0gkFaTX1rNhY8uJJvvpgQAACAASURBVL+0hscuyeTI/ilelyT7oZAsIW9bUSWzP9vCC1m5VNf7mTQghSsm9WXSgJTQ+KJuBrU+Px98WcDzi3P4ZGMh1sLhfTtz7vh0ThqWFhotej4ffPghPP00DBkCN93kvuSffx6GDoXU1K+3TUyEuDjYuBEeeADWrHGP79ULLroIjjsOwkPgPbcB763dxa9fWoHPb7lv+mhOGNrV65IOLD8fnn0W3n4bZs2Cnj3d0Yq8PHeUIqrhiJMx0K2bu/zYY/DBB7Bzp/tsTZ7sAmf37t69DwC/H/77X7ff9OkDt93mbn/xRddy3HWf/4+EBIiPh+xs+Oc/3Xuur4f0dLff/OAHEKElElpLYXktFz22kC1FlfzrorEcPzgE9p12RiFZQtaSbXt47NNs3lmzk/Aww+mjezBzUp92PxhiR0k1ryzN5YWsXLYXV5EQHcHU0d05NzODkelJwffDwedzYeXZZ1146dcPLr8cjjii6c9RWwvz57ugkJsLzz3nWtSkxdT6/Nz91npmz9/C8B6J/PP8sfROCfJlpnNzYc4ceP9913I6ZQpceil0PIh+oXl58PLL8PrrcNJJ8Mtftli538nvh3ffhWeecTX17u32m0mTmv4cdXXw+eduv9m61T1XVwW11rSnso4fzV7El/ll/P38MZwyIs3rkmQfCskSUvwBy7trdvLop9ks3V5CUmwkF07sySVH9KZrYozX5QWVQMCycEsxL2bl8ObqfGrqAwzulsA5mRmcMbo7neODZPWnxx+Hp55yLXgXXwyHH/79u0xY61rJ+vVzl2fNci1+GRnNW3M7tyq3lOteXMH6XeVcekRvbjplcPB3aaqqgunTXcvp1Knu8qH8kCoqcp/TTp3cEY1Fi1zXoNY6gvHMM651e8AAt98cdVTz7Dfg9snjjnPBW1pcWU09lz++mKXb93DvtFGcPS7d65KkgUKyhITKWh8vZuUwe/5WthdX0bNTB2Yc1YdzMtPpEKXDgwdSVlPPf1bs4IXFOazILSUy3HDMwC5MGdmNHwzpSmJMyy6a8j8CATegKDkZKipg+XI48sjm7U+8axdccYVrZb7qKjjzTPVXPkR1vgD/+HAjD360mZT4KO46a0TwHyIuKnJB1hhYuBAGDnSfu+b02GMutA4d6roJpbdQyLHW9Tfu1MmF/qws13LcnJ/roiLXd7myEmbMcH3/1V+5xVXV+bjiqSzmbyrizjOHc+FEDUwOBgrJEtR2ldXwxOdbeXbhdkqr6xnXK5krJvXhhKHd2t4sDq1k3c4yXszK5Y2V+ewsqyEqPIxJA1I4ZUQaPxzalaTYFg7MO3fCXXe5L/mHHmrZPpBFRfDnP8OCBW5w3403Bt+AqxDx6cZCbpu3huzCSs4em86tpw5t8RUpD4m1rkvEgw/Cdde5/rYt6cMP3eC/ujq4+mo4/fTmDa8FBfCnP7nBeY8+CpEt+G9fUgJ/+Qt89hmMHOn2mzR1A2hpNfV+rnlmKR+uK+CWKUOYOUkzkXhNIVmC0todZcz6bAvzVuThD1hOGt6NmZP6MrZnM7cAtWOBgGVZTglvrsrnrVX57CitISLMMK5XMscMSuWYgakMTUts3j7MK1bAb3/r+lP+9KeuK0RLt+5aC2++6Qb4paa6Q8lqGWuy3D1V3PnGl7y1eie9Onfg9qnDOG5QkP/QqK+Hv/7V9XUfNw6uv751fhzt3g333AOLF7tgOXly8zzvmjVwyy3uqMg117i+1K2x37z3Hvztb5CU5LpEaVBfi6vzBfjF88t4c9VOfnXCQK49PkRmimmjFJIlaAQClg/XFTDrsy18kV1Eh6hwzs3M4PIj+4TudGYhwlrL8pwS3l27i4/XF7I2381JnJoQzZH9OjOudyfG905mYJcEwr5vC/5bb7ng0r07/PGP0KNHM76DJsjJcV08hg9v3dcNUbsravnnh5t4duF2wsLg2uP6M3NSX2Iig7zvcWkp3HorrFwJl1ziTq0ZMqx1AwOPO655QuX778Pdd7uQf9ddbiaO1rRjh5sybtSo1n3ddsznD3D9Syt5ZVkePz62H9dPHqSg7BGFZPFcZa2Pl5bk8vj8LWwtqqJ7UgyXHNGb88b3DO7DuW1YQVkNn2zczUfrC1i4pZjC8loAEmIiGNszmeE9EhmSlsjgbon0SYk7cNcXn88dgk5OdlNUxce3wrv4Dk895Q6Lz5ihfsqNFJbX8sTnW3h8/lZqfQHOGZfOz34wgO4dY70urWk+/RT+8Ae44QY4/nhva9mzB373O/j5z930bAfL74ef/ARiY93zJHo8c8/cue49XXWVjsa0sEDAcstrq3l24XYuPaI3t5469Ps3UMj3ppAsnskrqeapz7fy3KLtlNX4GJ3RkRlH9eGk4d2IDNcf4GBhrSWnuJrFW4vJ2raHpdv2sKmwAn/A/X2Ijgijf5d4enbqQM9OHchoOKUlxdApwpIcF014TLTr55iQ4P38xda6Fu3XX4djjnGLKrRk/84QkV1YwWOfbeGlJbnU+wOcMiKNX54wkH6pHv+gaarycvf5AtcXPRimAMzOdl09qqrgjjvcYjdNUV/vBrdGN+w38fHed3Ww1s2t/MorbnrGW2919UmLsdbyhze+ZNZnW5iemcEfzxqhsTit7JBDsjHmJOBvQDjwmLX2T43ujwaeAsYBRcB0a+1WY0xv4EtgfcOmC6y1V3/Xaykkh75AwPLJxkKeXrCdD9ftwhjDScO7MeOoPupvHEJqfX42FVSwLr+cdTvL2FhQQU5xFTl7qqnzBQCIqa/h+o+fojI6licmX07nhGiSO0QRFx1BbFQ4cVHhdIiKoENUOB2iwokIDyMizBBmDBHhhvAwQ0SYITwsrOHcEBnecD3cELn3PNwQERbWcP/X90VFhNEhOpy4qIj//WKx1i228NBDLrj8/vcQ0/6mEKyp9/PW6nzmLsph4ZZioiLCmDYunSsm9aVPsM95vK/5810XnrvucgPNgklhoeufvH27C5YHmse4psb12w8LcwP1gu1Ix6uvwj/+4f6d//hH6KCucC3JWst9723g7x9uYsrINP567qjgn26xDTmkkGyMCQc2ACcAucBi4Hxr7dp9trkGGGmtvdoYcx5wprV2ekNIft1a2+QOggrJoWt3RS0vZOXw3KLt5BRXkxIfxbmZGVwwsSfpyfoj21YEApaC8lpyt+0k5Q+3Ep29mc/OuYJlA8dRVFHLnqp6quv8VNb5qKr1U1Xno6rOjy/QsketYiPDiYuOID7ancdFR5AYE8nhGxfzw9efpLz/YDbceAedEmLoHBdNaoI7tcVWm1qfn883F/Hump28vjKf8hofvTp3YPr4DM4Zl0FqQoi1Dr7/vgvHAwe6QXN7W5ODSXm5mxpu7VrXFeTbFsqpqHCBeu1a+PWv4eSTW7fOpvrgAxeQBw50rcteHx1qBx7+eDN3vbWOiX068ciPMlt+FiIBDj0kHw7cbq2d3HD9JgBr7V37bPNOwzZfGGMigJ1AKtALheQ2LRCwLNpazDMLt/P26nzq/ZbD+nbiwom9mDysG1ER6lLRJhUVuSm38vLg9tubtHJenS+AP2DxBfae26/P/e52X8Di87vb6wMBfH6Lzx+gPuDO997vCwSo91vqfAGq6nxU1PqorPVRUeun8qvLPkqr6ymurKPvuqXE1FTzUb9v/h0MDzN0SYimW1IMaUkxdEuMded7ryfF0DUxJiS6BuXuqWJhdjEfbyjkv+sKKK/1ERcVzglDuzJ9fE8m9ukUmv0dX3vNzb4werQLn8HcqlldDQ8/7FbF21/f4pISF4y3bnUtyUcf3eolHpTPP3fT0Z16qteVtBuvLc/juhdX0CcljicumxA64wRC2KGG5GnASdbamQ3XLwYmWmuv3Web1Q3b5DZc3wxMBOKBNbiW6DLgFmvtp/t5jSuBKwF69uw5btu2bQf9JqV15RRX8crSPF5e2rAsckwE08alc+HEnvTvEoStPNJ8rHVTVG3bBnfe6eYmDnLWWspqfBRX1lGzYCE7k7qwIyqBnaU15JfWNJxXk19aQ1Wd/xuPNQZS46NJS4ohLSmWbkkxdO8YQ7ekWLp7FKTLaupZl1/Ol/llrMgtYWF2MXkl1QB0jovih0O6Mnl4V47olxL8M1V8l6wsFyqPOMINBo2K8rqipquvh08++XruZmvd4L4NG1zXn/Hjva3vYC1b5uZR7tbN60ravM837eaqOUvoEB3OE5dNYEiax4M52zgvQ3I5EG+tLTLGjAP+DQyz1pZ92+upJTl4Vdb6eGv1Tl5aksOC7GKMgSP6dWbauHROGpZGbFQIfxnLwdmwwY3KHzLE60oOTm0tnH++m0ngvvv+Z15day3ltb6vwnN+STU7y2rIL6khv8xdzy+toaLW943HfRWkO8aSlhhD18RoOnaIomOHyIZTFB1jI0mMjSQmMpyo8DCiI8OICg8jzJhvtKL7/AHKanyUVtdRUuVawneU1JC7p4q8kmq2FVV9FYgBOsVFMaF3Jw7r24mJfTszqOshTOEXbPbOfz15sveD2g7Wq6/C3/8Ol17qpqgDN8ivshJGjPC0tINWXw8XXuj6UN93nxYdaQXrdpZx6ezFVNT6ePjicRzZP8Xrktosz7pb2EZPboz5CLjOWvutKVghObj4A5aF2UW8siyPN1flU1Xnp3fnDkwbl86ZY9PpoUNB7Ud+PnzxBZx1lteVHJovv3RdRZKS3Opp32MBivKaenaW1rCjtIadpdXsKKlpuF7NztIadpXVUFbjO/ATHYQuCdGkJ8eSntyBQd0SGJrmpujrmhjd9uZXffll13ocymEsEIB774V589z8w//8Z2hPqbZxI/zqV+4H5v33h/b/TYjIL63m0tmL2VxYwe/PGM75E1p5/ux24lBDcgSuu8QPgDzcwL0LrLVr9tnmJ8CIfQbunWWtPdcYkwoUW2v9xpi+wKcN2xV/2+spJHvPWsvS7SX8Z8UO3liVT2F5LfHREZw6Mo1p49IZ1yu57X0py3fLy4P/+z83Kv/JJ91cyKFs3Tp3GD8+3rWMtcAhZH/AUlZdz56qOkqq6ymtqqe0up5an586X4DahpO19qtZPyLCDBHhYSTERHzV+pzcIYouidGh3W2iqax1n68nn4TzznNz9Yay/Hw3MC8nx/VB/vnPg28mi4OxaZMLytHRbr9p7cWC2qGymnqufXYZn2wo5NIjenPLlCFEhMAYiVDSHFPAnQLcj5sCbra19k5jzB1AlrV2njEmBpgDjAGKgfOstdnGmLOBO4B6IADcZq39z3e9lkKyN6y1rNlRxn9W7uD1FfnklVQTFRHG8YO6cNqo7hw/uIu6U7RXubkuINfVubmH+/XzuqLmsWGDa1E+4ww30Eq8Za1bTnzOHDjpJPcjJpRbXvPz4Re/cPMnDx3qjmDMng0pIX7YPDsbfvlL93909XfO6CrNxOcP8Mc31zF7/hYmDUjhn+eP1SJczUiLici32lRQzrwV+by+YgfZuyuJCDNMGpDCaaO6c8LQriTEaEds17Zvd1+IPp8LyH37el1R88rPd63IxriQFsqtfKHMWnjsMXj2WTjlFPfjJZT/L/Ly3H5TXQ1/+Qv07w87d37dRSHUP2s7d7puSmFhof9eQsjzi7dzy79Xk5HcgccuyaRvqCwCFOS+KySH2EgIaQ5bdlfy5qp8Xl+Zz5f5ZRgDh/XpzMxJfTl5eDeS40JoBLm0rI0b3Zfg/fdD795eV9P89oaW/Hw3vdhNN0F6urc1tUd1dbB0KZx2mjtqEeqha9Omr4+89O/vbtv7WXvqKSgrc0tRh+r73Ns9qaDArTJ43XVt8+9DkJk+vid9UuK5+uklnPHAfP5+/hiOHXTwYyqk6dSS3E40DsYAY3p2ZOqo7kwZkUaXxPa3Epl8h/r6r5dxrqoK7rlpm8OWLa7lb+/o/Z4aINMqrHVHKSIjXatrTEzoBkc48H5jLTz4ILz0kuvm87Ofhfb7zclxP2r8ftdi3taONAWpnOIqrngqi/W7yvnp8QP4+Q8GtMlFkVqLulu0U9mFFby5Kp83Vu38KhiP7dmRU0akccqINE1SLvuXne1aVG+4AcaO9bqa1rNtm/vCB9cCqJaxlmWtm/Fh2za3NHOoTfHW2LZtcP317sfWxInfvp218K9/wQsvwNSprt9yKAfl3Fz3Hurrv+5aIi2uus7PLf9ezctLczmqfwp/O280neNDbCXNIKGQ3I4oGMshWbvWLZkbE9M+R6/v7YNdX+/ev1rGWkYg4ALVm2/Cuee6AWChHBQ3bnQDDcPDm3YkYt8+2Cef7MJ1KNvbB7uqyv3AHDDA64raBWstL2Tl8NvX1tCpQxQPXDiGcb06eV1WyFGf5DZubzB+fWU+63aWAzCuVzK/PXUoJw/vpmAsTbN0KfzmN9Cpkwsw7XFlrZ493RLIDz/cPt9/a/D53EqNH30EP/qRW2wjlAPyypXuyEtCgttvmvLD0hiYOdNNpZaa2vI1trQePdx+8+CD2m9akTGG6eN7Mqx7Etc8s5TpDy/g5z8YwI+P7adp4pqJWpJD1ObCCt5cmc8bq74ZjE8ZkaZgLAdv82bXmpeR4RZA6NzZ64qCQ20trF8PI0d6XUnbcddd8O67bmnzc87xuppDs20bXHmlC4Z//vOhBd5Vq1xXhdg28Le7rs5NeTdqlNeVtBul1fXc8u/V/GfFDsb1Sua+c0fTs3MbH0vSTNTdoo34tmA8ZUQaJ4/oRlpSG/jjKt6w1h36nTrVtYiJ8+CD8MorcPPNcPzxXlfTNmRnu+4Jkyd7XcmhsxbmznVdJjp2/P7PU1rqFk/p29f1zw71fXBvV5Ibbmgb/88h5LXledzy79UEApbbpg7jnHHpWvzrABSSQ9j+gnHm3hZjBWM5FNa6UfZHHaUlZr9NZaULyCtXutXfpk8P7a4BXtm503WvOO88rys5dNbCq6/ChAnNO13gp5+66dTS0uDuu0N7n6yudisMLlniFuq56CLtN60or6SaXz6/nIVbijlhaFf+cMZwumoGq2+lkBxi9tfHOLNXMlNGpnHy8DS6JenDLofI53P9J99+Gy680PWPlP2rrXWtex99BKee6pYWDvWZGFrTunWuz67P51bUC+UV53w++Pvf4T//cQMOf/zj5n3+lSvhllvcVIR33gnDhjXv87em+nrXdeu991xr8nXXab9pRf6AZdZn2fzl3Q1EhYdx4ymDOX98T8I0Vdz/UEgOAVt3V/LGqnzeWJnP2oZZKRSMpUWUl8Ntt8GyZXDJJe6kVp7vZi3MmgVvvOEG9XXRBP5N8sknLux17ux+aITy/NOVlXD77ZCV5X5YzpjRMvtNbq6bYebII5s/hLc2a93iKa+84qa9C+XW8RC1dXclN72yii+yi5jQuxN/PGs4/buEeHeeZqaQHKS2F1W5YLxqB6vzvp6ubcrI7pyirhTSEgoK3FRVO3a48xNP9Lqi0FJaCklJbgqz3bsVlr/Lyy/DAw/A0KFuNcND6bPrtaIit79s3w6/+pXrg9ySysshLs61KO/a5T5nofxDdu9+Y637G9S1q9cVtSvWWl5cksudb3xJZa2PS47ozc9+MICk2EivSwsKmgIuiOQUVzXMY5zPytxSAEZndOSWKUM4eUQaPTQrhbSkpCQ3Av9Xv9KMDd9HUpI7nzsXnnnGzW97zDHe1hSsUlPhuOPcv1F0iC9ykJjouon85CcwblzLv97egXuVlXDttW5A3803f/35CzV76375ZXdE5rrr4Ac/8LamdsQYw7mZGRw/uAt/eXc9s+dv4d/L8rhu8iDOzczQan3fQS3JrSCnuIq3V+/k9VX5rMgpAWBkehJTGhb4yOikaVqkBdXVuUA3fXrbX166tezaBb/7nZvm6qyz3PR5kWqVYe1at1RxW5jRwOeDp5+Gs8/2brYJa+H1110/6ORkuPVWGD7cm1qaw+7dbnDiqlVw2mnuB0BUlNdVtTur80r53X/WsHjrHgZ1TeBXJw7khKFd2+0sGOpu0cqstWzYVcE7a3byzpqdrNnhulIM75HIlBHdmTIiTfMXSuvIzXVfShs3ugFBar1pPj6f65/80ksweLALMO21z+XemVIeftj9Gzz+eGgP0srPd/vNunWuJbylu1ccyMaNbhzBrl1wxRWhPcuKz+dak+fOdfNC33Zb884SIk1ireWNVfn89d0NZO+uZFRGR3594iCO7N+53YVlheRWEAhYluWU8G5DMN5aVAW4PsaTh3Vj8rBu9E6J87hKaTd8PnjhBXjiCbfE9I03whFHeF1V2/Tpp24U/733wqBBXlfT+rZvd+999Wo32OyGG0J3nt9A4OsuAeHh7r0cfbTXVTmVlXDPPW62lbvuCt2QvNcXX7ip7v7wh9BuHQ9xPn+AV5bmcf/7G9hRWsPYnh25+ph+/HBI13YzE4ZCcgupqvPxxeYiPlxXwHtrd1FQXktEmOHwfp2ZPKwbJw7tShfNTShe+Mc/3IjySZPclGVaQa9lVVd/vVLa3LluDt2+fb2tqTWUlLi5j6Oi3KHzE04I7fD28MPu/++II+AXvwi+JaOthZoa91nbudOtXHjBBaHbar/vfvPiizB6NAwY4G1N7VStz8/zi3N45JNscvdU0y81jquO7sfU0d2JiQz3urwWpZDcjLbsruS/6wr47/oCFm4pps4XoENUOMcMTGXysG4cN7iLRoyKN8rL3ZdOly7uC3TDhuBpBWsvysvhRz9y5+ec46YKi4/3uqrmt3kz9OvnLr/3HmRmuj6zoaiiwrXSdu3q+syuWgXHHhv8Yf/55920ar17fz2gMNhr/jZVVW6/2bPH9QG/6CI3WFJanc8f4I1V+fzr42y+zC8juUMk52RmcMGEnm32aLhC8iEorapnwZYiPt+0m482FLKtoRtFv9Q4jhvUhWMHdWF8n2SiI9r2Ly0JYjU17hDxc8+5lpg//MHritq3sjJ46CF45x0XkC+4AM48M/RneABYvx4efdStpPbQQ64vdqiqrXUr5z37rHsf99zjdUUHb8EC+Nvf3I/iMWPgyitD9/+kvBweecTNRR4bC+efD9Omue5i0uqstXyxuYinF27j3TW78AUsR/VP4ayxPThxWDfio0P06MV+KCQfhPKaehZvLeaLzUV8vrmItfllWAsxkWEc3rczxw3uwrEDu2jgnXivpgbeegvmzHEtMEcc4VbO69PH68oEXGvro4/CihUuiIVqSytAdrb7nH30kZvO66KLYOrU0JyZoLbW/YB56ik3//HEiW6/6d/f68q+n/p6twLgnDluYO6113pd0aHZssX1CV+82M0uEmxdXtqhXWU1PL84h+cX55BXUk10RBg/HNqVqaO6c+yg1JBvJFRI/hbWWnL3VLMsp4Rl2/ewdNseVu8owx+wREWEMbZnRw7vm8Lh/TozKiMp5D8I0sbMmQOzZ7v5jq+8MrSXsG3LCgvdF721bkGK/v3hjDOgWzevK2uaujp3CNzvd0sxn3tuaE8luLebwrBhbr9pK/OFV1W5gYfx8S5gPvWUm55w0qTQ7LO8d78Bt6x5Robbb7p397audsxay9Lte3ht+Q5eX5lPcWUdcVHhHNk/hR8M6cJxg7qE5DgshWTcf25heS1r88tYs6OM5TklLNtewu6KWgBiI8MZkZ7ExD6dOLxfZ8b2TG7zndUlhFRXu0Or770Hp5wCRx3lBk3l5rov+1Dti9ie1NS4pZk//dQF5sMPdwPdDjssuA4pb93qWoxXr3azVhjjljDv3z/0Zq2oqYGFC91+c/zx7lRW5lorR45su/vN/PlutcP8fDdod8oU18+6d+/Qe891dW6/+eQT9yNg4kS3Uuhhh3096E9aXb0/wPxNu3l37S7+u66A/NIawE11e2S/FCb06URmr04kdQj+MVrtLiRX1PrYUljJxoJyvswv48t8d15UWffVNn1T4hjdsyNjeiYzJqMjg7slEBEedqjlizQfa+HDD+Hjj90XfV0ddOoEV12l5aRDWUGBOzz+xhuum8zNN7uwXFXlvvS9CDH5+fD22+6ztm2bq2HcOPjtb0NvAJW1LuR//LH7YVlb65bEnjnThcX2IhCARYvcLDdZWW5A73PPuf/bqqrQOxpQVATz5rn9Zu8y4aec4hoQYmJCL/y3IdZa1u0s58N1BXy8vpDlOSXU+QMYA4O6JjC+dydGpCcxrHsiA7okEBURXFmrzYVkay3FlXXsKKkhr6SanOIqsndXkl1YwZbdlRSU1361bVREGAO7xjOkWyJD0hIZ2j2RId0SQ+LXjbRz1sKll7qR98cc407Dh0NYcP2Bke8pEICVK2HgQBdYli6FsWO9qeXTT92iDqNHu8/ZpEnuB1mouuIKF6SOPtq1oI4c2b73mz17YMcOd9TJ73czeIwe7XVV308g4I5y9OnjjmwsXeoGLSokB42aej/Lc0pYvKWYRVuLWbptD5V1fgCiwsMY1C2BQd0S6JsaR7/UePqlxtGzU5xn4fmQQ7Ix5iTgb0A48Ji19k+N7o8GngLGAUXAdGvt1ob7bgJmAH7gZ9bad77rtUaPHWfnvvlfiirq2F1RS1FFHUWVtRSU1ZJfWsOOkmrySqqp9QW+8bhOcVH0SYmjT0ocfVPj6JsSR9/UePqmxKmFWEJXYSGkpOgLoD2orIQ4j6ZYqqtzrx/Kgwv3tXu3C/ntORh/m/p6FzTbwmwr4O1+I00SCFi2FlWyekcZa/JKWZVXysaCCgr3adAMDzOkJ8eSlhRD96RY0jrGkJYUS/eOMXRLjCUlPoqOHaJaJEgfUkg2xoQDG4ATgFxgMXC+tXbtPttcA4y01l5tjDkPONNaO90YMxR4DpgAdAfeBwZaa/3f9nrRaQNs2iX3f+O28DBDSnwU3TvG0r1jLD06xtI9Kear6+nJsXTsEIKjrEVERETaobKaerYUVpK9u4Lswkq27K4kv7SG/JJqdpXX9MNsRAAACJhJREFU4g/8bz5NiI4gOS6K5LgoOnWIJLlDFHHREXSIDic+KoIO0RHER4fTISqCuOhw4qIiiIkMJyoijMjwMKLCwxoum69ui42K+NaQ3JQhrxOATdbabABjzFzgdGDtPtucDtzecPkl4J/GLf59OjDXWlsLbDHGbGp4vi++7cW6d4zlnxeMpXN8FCnxUXSOiyYpNrLdLI8oIiIi0tYlxkQyKqMjozI6/s99/oCbbGFHaTU7S2sorqxjT2UdxVV7z+vZXVHHxoIKKmt9VNb5qWvUw6A5NCUk9wBy9rmeC0z8tm2stT5jTCnQueH2BY0e26PxCxhjrgSuBOjZsydTRqY1tX4RERERaUPCwwzdkmLoltT0mX/q/QGqav1U1vm+Cs6VtT5q6l2ArvMHqPMFqPdb6nx+d+4PcO3d3/6cQTF5orX2EeAR+P/27i9ErrMO4/j36Ww1ZgbTheqF2WAiBCU0NduEUl3wZkvo1pLceBFF8c+V4J8qBbEN8cJbpSoogrR606hgrJAVo1upV4HWxuw2axsjsUKbWLFeqGWFlI0/L87ZZnMys5mZndn33TnPBwZ2zp7znue87zln3pl5z5ziwr3EcczMzMxsE7m1cQvbtt7S8w8zrHX7nW5GQF8Gdqx6PlFOazuPpDFgG8UFfN0sa2ZmZmaWlW4+SX4O2C1pF0UH9wjwsco8J4FPUow1/gjwdESEpJPAjyU9SnHh3m7g94MKb927evUqp06dYn5+nsnJSWZmZmg06nGzlGFve7X8gwcPMjc3V8u67kbO+2LO2VIYhfrodRs6zT/Iulhd1p3lHf/OnTvX9vzRy/mk34wbvVy/BrG+9ZQxCseD9SgibvoA7qf4hYu/AEfLaV8HDpV/bwF+Blyk6AS/Z9WyR8vlLgAzN1vX/v37wwZreXk5pqeno9VqhaRotVoxPT0dy8vLqaMN3bC3vVp+s9mM8fHxWtZ1N3LeF3POlsIo1Eev29Bp/itXrgysLqrraDQa0Wg0AohWqxXj4+PRbDZ7Pp/0214bvVy/BrG+9ZQxCseDtQeciU79307/SPVwJ3nwZmdno9VqBfDmo9VqxezsbOpoQzfsbW9XfvVRl7ruRs77Ys7ZUhiF+uh1GzrNf+zYsYHVRTfnjH7OJ/2210Yv169BrG89ZYzC8WDtrdVJ9i+t18D8/DxLS0vXTVtaWmJhYSFRoo0z7G1vV35VXeq6GznvizlnS2EU6qPXbeg0/+nTpwdWF92cM9bSab39ttdGL9evQaxvPWWMwvFgvXMnuQYmJydpVu5I1Gw22bdZb0vag2Fve7vyq+pS193IeV/MOVsKo1AfvW5Dp/mnpqYGVhfdnDPW0mm9/bbXRi/Xr0Gsbz1ljMLxYH3o9BFzqoeHWwxencdSpRqTvDKmsE513Y2c98Wcs6UwCvWx2cYkV88fvZxPPCZ5uGWMwvFg7bHGcIub3pZ6ox04cCDOnDmTOsbIWbkqd2FhgX379tXqqtxhb3u1/JWr0etY193IeV/MOVsKo1AfvW5Dp/kHWRery9q7dy8Ai4uLbc8fvZxP+s240cv1axDrW08Zo3A82I0kdbwttTvJZmZmZlZLa3WSPSbZzMzMzKzCnWQzMzMzs4rshltIep3ixiOWn9uBf6YOYW25bfLltsmb2ydfbpt8jVLbvDsi3tHuH93clnqjXeg0NsTSknTGbZMnt02+3DZ5c/vky22Tr7q0jYdbmJmZmZlVuJNsZmZmZlaRYyf5B6kDWEdum3y5bfLltsmb2ydfbpt81aJtsrtwz8zMzMwstRw/STYzMzMzS8qdZDMzMzOziqw7yZIekhSSbk+dxQqSviHpT5LOSfqFpNtSZ6o7SfdJuiDpoqSvps5jBUk7JP1O0ouSXpD0YOpMdj1JDUnzkn6ZOotdI+k2SSfK15rzkj6QOpMVJH25PJ/9UdJPJG1JnWmYsu0kS9oBHAReTp3FrvMUcEdE3An8GXg4cZ5ak9QAvgfMAHuAj0rakzaVlZaBhyJiD3AP8Dm3TXYeBM6nDmE3+A7w64h4H/B+3EZZkLQd+CJwICLuABrAkbSphivbTjLwLeArgK8szEhEzEXEcvn0GWAiZR7jbuBiRLwUEW8APwUOJ85kQES8GhFny79fp3ih3542la2QNAF8GHgsdRa7RtI24EPA4wAR8UZE/CttKltlDHibpDFgK/C3xHmGKstOsqTDwOWIeD51FlvTZ4BTqUPU3HbglVXPL+GOWHYk7QQmgWfTJrFVvk3xQcz/Ugex6+wCXgN+VA6FeUxSM3Uog4i4DHyT4hv+V4F/R8Rc2lTDlayTLOm35ZiW6uMw8AjwtVTZ6u4mbbMyz1GKr5OPp0tqlj9JLeDnwJci4j+p8xhIegD4R0T8IXUWu8EYcBfw/YiYBJYAX2uRAUnjFN9U7gLeBTQlfTxtquEaS7XiiLi33XRJeyka4HlJUHydf1bS3RHx9w2MWFud2maFpE8BDwDT4R/aTu0ysGPV84lymmVA0q0UHeTjEfFk6jz2pingkKT7gS3A2yU9EREj/YK/SVwCLkXEyrcuJ3AnORf3An+NiNcAJD0JfBB4ImmqIcpuuEVELEbEOyNiZ0TspDhg7nIHOQ+S7qP4ivJQRPw3dR7jOWC3pF2S3kJxEcXJxJkMUPEu/3HgfEQ8mjqPXRMRD0fERPkacwR42h3kPJSv9a9Iem85aRp4MWEku+Zl4B5JW8vz2zQjflFlsk+SbdP6LvBW4Knyk/5nIuKzaSPVV0QsS/o88BuKK41/GBEvJI5lhSngE8CipIVy2iMR8auEmcw2gy8Ax8s3/i8Bn06cx4CIeFbSCeAsxXDLeUb89tS+LbWZmZmZWUV2wy3MzMzMzFJzJ9nMzMzMrMKdZDMzMzOzCneSzczMzMwq3Ek2MzMzM6twJ9nMzMzMrMKdZDMzMzOziv8Dpy2EIrGFXgIAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 864x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "from sklearn.neighbors import KernelDensity\n",
    "from matplotlib import pyplot as plt\n",
    "\n",
    "kde = KernelDensity(kernel='gaussian', bandwidth=0.75).fit(X) # fit a KDE model\n",
    "x_ticks = np.linspace(-5, 10, 1000)[:, np.newaxis] # choose 1000 points on x-axis\n",
    "log_density = kde.score_samples(x_ticks) # compute density at 1000 points\n",
    "gaussian_kernel = lambda z : lambda x: np.exp(-np.abs(x-z)**2/(0.75**2)) # gaussian kernel\n",
    "kernel_linspace = lambda x : np.linspace(x-1.2,x+1.2,30)\n",
    "\n",
    "plt.figure(figsize=(12,4))\n",
    "plt.plot(x_ticks[:, 0], np.exp(log_density)) # plot the density estimate\n",
    "plt.plot(X[:, 0], np.full(X.shape[0], -0.01), '.k', markersize=10) # plot the points in X\n",
    "plt.plot(kernel_linspace(4), 0.07*gaussian_kernel(4)(kernel_linspace(4)), '--', color='r', alpha=0.75)\n",
    "plt.plot(kernel_linspace(5), 0.07*gaussian_kernel(5)(kernel_linspace(5)), '--', color='r', alpha=0.75)\n",
    "plt.plot(kernel_linspace(1), 0.07*gaussian_kernel(1)(kernel_linspace(1)), '--', color='r', alpha=0.75)\n",
    "plt.xlim(-4, 9)\n",
    "plt.ylim(-0.02, 0.32)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "Pros of KDE:\n",
    "* Can approximate any data distribution arbitrarily well.\n",
    "\n",
    "Cons:\n",
    "* Need to store entire dataset to make queries, which is computationally prohibitive.\n",
    "* Number of data needed scale exponentially with dimension (\"curse of dimensionality\")."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "# Clustering\n",
    "\n",
    "Clustering is the problem of identifying distinct components in the data.\n",
    "* A cluster $C_k \\subseteq \\mathcal{X}$ can be thought of as a subset of the space $\\mathcal{X}$.\n",
    "* Datapoints in a cluster are more similar to each other than to points in other clusters"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "We can perform clustering via density estimation with a GMM model.\n",
    "$$P_\\theta (x,z) = P_\\theta (x | z) P_\\theta (z)$$\n",
    "* $z \\in \\mathcal{Z} = \\{1,2,\\ldots,K\\}$ is discrete and follows a categorical distribution $P_\\theta(z=k) = \\phi_k$.\n",
    "* $x \\in \\mathbb{R}$ is continuous; conditioned on $z=k$, it follows a Normal distribution $P_\\theta(x | z=k) = \\mathcal{N}(\\mu_k, \\Sigma_k)$.\n",
    "\n",
    "The parameters $\\theta$ are the $\\mu_k, \\Sigma_k, \\phi_k$ for all $k=1,2,\\ldots,K$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Intuitively, a GMM represents well the two clusters in the geyser dataset:\n",
    "\n",
    "Raw data | Single Gaussian | Mixture of Gaussians\n",
    "--|--|---\n",
    "<img width=90% src=\"img/oldfaithful_v2.png\"> | <img width=90% src=\"img/oldfSingle_v2.png\"> | <img width=90% src=\"img/oldfMOG_v2.png\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "<center><img width=40% src=\"img/algorithms11u.png\"></center>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "<center><img width=40% src=\"img/algorithms12u.png\"></center>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Linear Dimensionality Reduction\n",
    "\n",
    "Suppose ${X} = \\mathbb{R}^d$ and ${Z} = \\mathbb{R}^p$ for some $p < d$. The transformation \n",
    "$$f_\\theta : {X} \\to {Z}$$\n",
    "is a linear function with parameters $\\theta = W \\in \\mathbb{R}^{d \\times p}$:\n",
    "$$ z = f_\\theta(x) = W^\\top \\cdot x. $$\n",
    "The latent dimension $z$ is obtained from $x$ via a matrix $W$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "In this example, the data lives in a lower-dimensional 2D plane within a 3D space (image [credit](https://doc.plob.org/machine_learning/14_Dimensionality_Reduction.html)).\n",
    "\n",
    "<center><img width=50% src=\"img/pca_example_plane.png\"></center>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "Even linear dimensionality reduction is powerful. Here, in uncovers the geography of European countries from only DNA data\n",
    "\n",
    "<center><img width=50% src=\"img/dna_map.jpg\"></center>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Principal component analysis (PCA) assumes that \n",
    "* The projection subspace is spanned by a set of orthonormal vectors $w^{(1)}, w^{(2)}, \\ldots, w^{(p)}$\n",
    "* The data $x$ are approximated by a linear combination $\\tilde x$ of the $w^{(k)}$\n",
    "$$ x \\approx \\tilde x = \\sum_{k=1}^p w^{(k)} z_k = W z $$\n",
    "for some $z \\in \\mathcal{X}$ that are the coordinates of $\\tilde x$ in the basis $W$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "A natural objective is to minimize the reconstruction error\n",
    "$$J_1(W) = \\sum_{i=1}^n \\| x^{(i)} - \\tilde x^{(i)} \\|_2^2 =\\sum_{i=1}^n \\| x^{(i)} - W W^\\top x^{(i)} \\|_2^2$$\n",
    "between each input $x^{(i)}$ and its approximate reconstruction $$\\tilde x^{(i)} = W \\cdot z^{(i)} = W\\cdot W^\\top \\cdot x^{(i)}.$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "The variance objective is simply\n",
    "$$J_2(W) = \\hat{\\mathbb{E}}\\left[ \\| W^\\top x \\|^2 \\right] = \\frac{1}{n} \\sum_{i=1}^n \\| W^\\top x^{(i)}\\|_2^2.$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "The two are equivalent (figure credit: [Alex Williams](http://alexhwilliams.info/itsneuronalblog/2016/03/27/pca/))\n",
    "\n",
    "<center><img width=80% src=\"img/pca_two_views.png\"></center>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Recall that the positive semidefinite matrix $\\hat \\Sigma = \\frac{1}{n} X^\\top X$ has an *eigendecomposition*\n",
    "$$\\hat \\Sigma = Q \\Lambda Q^\\top = \\sum_{j=1}^d \\lambda_j q^{(j)} (q^{(j)})^\\top. $$\n",
    "* $Q$ is a matrix whose columns are orthonormal eigenvectors $q^{(j)}$ for $j = 1,2,\\ldots,d$.\n",
    "* $\\Lambda$ is a diagonal matrix of positive eigenvalues $\\lambda_1 \\geq \\lambda_2 \\geq \\ldots \\geq \\lambda_d$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "The variance objective can be written as follows.\n",
    "$$\\max_w J(w) = \\max_w \\sum_{j=1}^d \\lambda_j (w^\\top q^{(j)})^2$$\n",
    "Its optimum is attained by the top eigenvector $w = q^{(1)}$. The optimum is $J( q^{(1)}) = \\lambda_1$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "More generally when $p>1$, our objective is\n",
    "$$J(W) = \\sum_{k=1}^p \\sum_{j=1}^d \\lambda_j ((w^{(k)})^\\top q^{(j)})^2$$\n",
    "where $W$ is a matrix of orthonormal columns $w^{(1)}, w^{(2)}, \\ldots, w^{(p)}$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "<center><img width=40% src=\"img/algorithms13u.png\"></center>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "<left><img width=25% src=\"img/cornell_tech2.svg\"></left>\n",
    "# Part 3: Machine Learning in Practice\n",
    "\n",
    "We conclude with high-level considerations about how to apply machine learning."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Machine learning is an interactive process. At each iteration, the machine learning engineer needs to make a number of decisions."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "* Add more data?\n",
    "* Train the algorithm for longer?\n",
    "* Use a bigger model?\n",
    "* Add regularization?\n",
    "* Add new features?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "We prioritize these using a principled process (more on this in a few weeks)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Datasets for Model Development\n",
    "\n",
    "When developing machine learning models, it is customary to work with three datasets:\n",
    "* __Training set__: Data on which we train our algorithms.\n",
    "* __Development set__ (validation or holdout set): Data used for tuning algorithms.\n",
    "* __Test set__: Data used to evaluate the final performance of the model."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Model Development Workflow\n",
    "\n",
    "The typical way in which these datasets are used is:\n",
    "1. __Training:__ Try a new model and fit it on the training set."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "2. __Model Selection__: Estimate performance on the development set using metrics. Based on results, try a new model idea in step #1."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "3. __Evaluation__: Finally, estimate real-world performance on test set."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# How To Decide Which Algorithm to Use\n",
    "\n",
    "One factor is how much data you have. In the __small data__ (<10,000) regime, consider:\n",
    "* Linear models with hand-crafted features (LASSO, LR, NB, SVMs)\n",
    "* Kernel methods often work best (e.g., SVM + RBF kernel)\n",
    "* Non-parametric methods (kernels, nearest neighbors) are also powerful"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "In the following lectures, we will see algorithms for the __big data__ regime."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Some additional advice:\n",
    "* If interpretability matters, use decision trees or LASSO.\n",
    "* When uncertainty estimates are important use probabilistic methods.\n",
    "* If you know the data generating process, use generative models."
   ]
  }
 ],
 "metadata": {
  "accelerator": "GPU",
  "celltoolbar": "Slideshow",
  "colab": {
   "collapsed_sections": [],
   "name": "neural-ode.ipynb",
   "provenance": []
  },
  "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.0"
  },
  "rise": {
   "controlsTutorial": false,
   "height": 900,
   "help": false,
   "margin": 0,
   "maxScale": 2,
   "minScale": 0.2,
   "progress": true,
   "scroll": true,
   "theme": "simple",
   "width": 1200
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}