Lectures/lecture_14_summary.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Lecture #14: Hamiltonian Monte Carlo\n",
    "## AM 207: Advanced Scientific Computing\n",
    "### Stochastic Methods for Data Analysis, Inference and Optimization\n",
    "### Fall, 2021"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "<img src=\"fig/logos.jpg\" style=\"height:150px;\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Lecture #14 Summary\n",
    "\n",
    "We use the ***Gibbs distribution*** to transform between probability density functions and energy functions\n",
    "$$\n",
    "U(q) = - \\log \\pi(q), \\quad \\pi(q) = \\frac{1}{Z}\\exp\\left\\{ \\frac{-U(q)}{T}\\right\\},\\; T=1\n",
    "$$\n",
    "This allows us to use gradient information when we sample from $\\pi(q)$.\n",
    "\n",
    "### Hamiltonian Monte Carlo: with Exact Integration\n",
    "\n",
    "Let $\\pi(q)$ be our target distribution with $q\\in \\mathbb{R}^D$. We turn $\\pi$ into an energy function $U(q)$ by $U(q) = - \\log(\\pi(q))$. We choose a kinetic energy function $K(p)$. We fix $t = T$.\n",
    "\n",
    "0. start with a random $q^{(0)}\\in\\mathbb{R}^D$\n",
    "1. repeat:\n",
    "  1. (**kick-off**) sample a random momentum from the Gibbs distribution of $K(p)$, i.e.\n",
    "    $p^{(current)} \\sim \\frac{1}{Z} \\exp\\{-K(p)\\}$.\n",
    "    \n",
    "  2. (**simulate movment**) simulate Hamiltonian motion for a time interval of length $T$, i.e. start at $(q^{(current)}, p^{(current)})$ and find $(q^{(\\text{time }T)}, p^{(\\text{time }T)})$. You need to integrate $\\frac{dq}{dt} = \\frac{\\partial H}{\\partial p}, \\frac{dp}{dt} = -\\frac{\\partial H}{\\partial q}$!\n",
    "  \n",
    "  3. accept $q^{(current)}$ as a new sample: $q^{(current)} \\leftarrow q^{(\\text{time }T)}$.<br><br>\n",
    "  \n",
    "<img src=\"fig/skate.jpg\" style=\"height:500px;\">\n",
    "\n",
    "\n",
    "<br><br>\n",
    "### Hamiltonian Monte Carlo: with the Leap Frog Integrator\n",
    "\n",
    "Let $\\pi(q)$ be our target distribution with $q\\in \\mathbb{R}^D$. We turn $\\pi$ into an energy function $U(q)$ by $U(q) = - \\log(\\pi(q))$. We choose a kinetic energy function $K(p)$.\n",
    "\n",
    "**0.** start with a random $q^{(0)}\\in\\mathbb{R}^D$<br>\n",
    "\n",
    "**1.** repeat:<br>\n",
    "$\\quad$**A.** (**kick-off**) sample a random momentum from the Gibbs distribution of $K(p)$, i.e.\n",
    "    $p^{(current)} \\sim \\frac{1}{Z} \\exp(-K(p))$.<br>\n",
    "    \n",
    "$\\quad$**B.** (**simulate movement**) simulate Hamiltonian motion for $L$ steps each with time interval $\\epsilon$, using the leap-frog integrator. <br>\n",
    "$\\quad$$\\quad$**a.** Repeat for T - 1 times, for $p^{(\\text{step } 0)} = p^{(current)}$, $q^{(\\text{step } 0)} = q^{(current)}$:<br>\n",
    "$\\quad$$\\quad$$\\quad$**i.** (**half-step update for momentum**) $p^{(\\text{step } t + 1/2)} \\leftarrow  p^{(\\text{step } t)} - \\epsilon/2 \\frac{\\partial U}{\\partial q}(q^{(\\text{step } t)})$<br>\n",
    "$\\quad$$\\quad$$\\quad$**ii.*** (**full-step update for position**) $q^{(\\text{step } t + 1)} \\leftarrow q^{(\\text{step } t)} + \\epsilon \\frac{\\partial K}{\\partial p}(p^{(\\text{step } t + 1/2}))$<br>\n",
    "$\\quad$$\\quad$$\\quad$**iii.** (**half-step update for momentum**)$p^{(\\text{step } t + 1)} \\leftarrow  p^{(\\text{step } t + 1/2)} - \\epsilon/2 \\frac{\\partial U}{\\partial q}(q^{(\\text{step } t + 1)})$<br>\n",
    "$\\quad$$\\quad$**b.** (**reverse momentum**) $p^{(\\text{step } T)} \\leftarrow -p^{(\\text{step } T)}$<br>\n",
    "\n",
    "$\\quad$**C.** (**correction for simulation error**) implement Metropolis-Hasting accept mechanism: <br>\n",
    "$\\quad$$\\quad$**a.** compute $\\alpha = \\min\\left(1, \\exp\\left\\{H(q^{(current)}, p^{(current)}) - H(q^{(\\text{step } T)}, p^{(\\text{step } T)})\\right\\} \\right)$<br>\n",
    "$\\quad$$\\quad$**b.** sample $U\\sim U(0, 1)$, if $U \\leq \\alpha$ then accept, else keep old sample<br><br>\n",
    "\n",
    "<img src=\"fig/hamiltonian.jpg\" style=\"height:300px;\">"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "celltoolbar": "Slideshow",
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}