lectures/07_pose_tracking.ipynb

{
 "cells": [
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Markerless Pose Tracking\n",
    "\n",
    "We've talked a lot about extracting neural signals from ephys and ophys measurements, but what comes next? One of the fundamental goals of systems neuroscience is to understand the relationship between neural activity and behavioral outputs. To that end, we need to quantify  and model not only neural activity but animals' behavior as well. \n",
    "\n",
    "The study of natural behavior is called **ethology** {cite}`Tinbergen1963-fg`. With advances in deep learning, the field is undergoing a computational revival {cite}`Anderson2014-lg`. The once tedious process of labeling video data to track body posture and annotate behaviors has been largely automated by machine learning methods, which are well-suited to the task {cite}`Branson2009-mz, Machado2015-yi, Mathis2018-uz, Graving2019-fq, pereira2022sleap, Wu2020-fi, Bala2020-ch`. At the intersection of computational neuroscience and computational ethology is the (predictably named) emerging field of **computational neuroethology** {cite}`Datta2019-ji`.\n",
    "\n",
    "This chapter develops the basic methods underlying most markerless pose tracking systems.  There are a few key ideas. First, we can cast pose tracking &mdash; i.e. the task of labeling keypoints of interest in video frames, like an animals paws, snout, and tail, and tracking them over time &mdash; as a **supervised learning** problem. Given a few labeled frames, we can train a classifier to predict the keypoint locations in new frames.  Second, such classifiers require surprisingly few training examples, particularly when they are given features from networks that have been pre-trained on similar tasks. For example, deep neural networks trained for image classification on very large datasets like ImageNet may not immediately solve pose tracking problems, but the features they've learned for classifying cats and dogs may still be useful for tracking paws and snouts.  Using a neural network trained on one task to jump-start a model for a similar task is called **transfer learning**. \n",
    "\n",
    "We'll start with a simple logistic regression model for pose tracking and show how it can be implemented with a convolutional neural network (CNN). Then we'll show how the same ideas can be generalized to pre-trained CNNs for image classification, like very deep residual networks {cite}`He2016-rp`. "
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Supervised learning\n",
    "\n",
    "In supervised learning problems, our data consists of a set of tuples $\\{(\\mathbf{x}_n, y_n)\\}_{n=1}^N$ where $\\mathbf{x}_n$ are the **inputs** and $y_n$ are the **outputs**. Contrast this with the spike sorting and calcium deconvolution problems, which we framed as unsupervised matrix factorization problems.\n",
    "\n",
    "A simple way to frame the pose tracking problem as a supervised learning problem is to chop each image frame into **patches** and then assign each patch a binary **label** to indicate whether or not the patch contains a specific key point (e.g. \"left paw\"). Then we can train a binary **classifier** to predict the labels given the image patches. Of course, we usually want to track multiple key points at once, and we could simply train separate classifiers for each one. Ideally, the trained classifiers will **generalize** to new image patches from new video frames, giving us predictions about which patches contain which keypoints. Then, in post-processing, we can determine the most likely configuration of key points in future frames, using the classifiers' predictions.\n",
    "\n",
    ":::{figure} images/07_pose_tracking/paw_classifier.png\n",
    "---\n",
    "name: paw_track\n",
    "---\n",
    "A simple way to do pose tracking is to carve each video frame into small patches and label them as positive or negative examples of a key point. Then train a classifier to predict the label for new patches.\n",
    "This figure was adapted from fig. 1B of {cite:t}`Machado2015-yi`\n",
    ":::"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Logistic regression\n",
    "\n",
    "Let $\\mathbf{x}_n \\in \\mathbb{R}^P$ be the $n$-th image patch, flattened into a vector of $P$ pixels. Let $y_n \\in \\{0,1\\}$ be a binary label specifying whether the key point of interest is present in that frame. \n",
    "\n",
    "In **logistic regression**, we model the conditional distribution of the label given the image as,\n",
    "\n",
    "$$\n",
    "p(y_n \\mid \\mathbf{x}_n) = \\mathrm{Bern}\\big( \\sigma(\\mathbf{w}^\\top \\mathbf{x}_n) \\big)\n",
    "$$\n",
    "\n",
    "where $\\mathbf{w} \\in \\mathbb{R}^P$ are the **weights** for each pixel, and $\\sigma: \\mathbb{R} \\mapsto (0, 1)$ is the **logistic function**\n"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```{admonition} The Bernoulli distribution\n",
    "The [**Bernoulli distribution**](https://en.wikipedia.org/wiki/Bernoulli_distribution) is a distribution over binary variables $y \\in \\{0,1\\}$ with probability $p \\in [0,1]$. Its pmf can be written as,\n",
    "\n",
    "$$\n",
    "\\mathrm{Bern}(y; p) = p^{y} \\, (1-p)^{(1-y)}\n",
    "$$\n",
    "```"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```{admonition} The logistic (sigmoid) function\n",
    "\n",
    "The [**logistic (aka sigmoid) function**](https://en.wikipedia.org/wiki/Sigmoid_function) is a map from the reals to the interval $(0,1)$. It is defined as,\n",
    "\n",
    "$$\n",
    "\\sigma(a) = \\frac{1}{1 + e^{-a}} = \\frac{e^a}{1 + e^a}\n",
    "$$\n",
    "\n",
    "It asymptotes at $\\lim_{a \\to -\\infty} \\sigma(a) = 0$ and $\\lim_{a \\to \\infty} \\sigma(a) = 1$.  It is plotted below.\n",
    "\n",
    "Interestingly, the logistic function is symmetric in that,\n",
    "\n",
    "$$\n",
    "1 - \\sigma(a) = \\frac{1}{1+e^a} = \\sigma(-a).\n",
    "$$\n",
    "\n",
    "Its derivative is,\n",
    "\n",
    "$$\n",
    "\\begin{align*}\n",
    "\\sigma'(a) = \\frac{e^a}{(1+e^a)^2} = \\left( \\frac{e^a}{1+e^a}\\right) \\left( \\frac{1}{1+e^a}\\right) = \\sigma(a) \\sigma(-a).\n",
    "\\end{align*}\n",
    "$$\n",
    "\n",
    "The derivative is positive (i.e., the logistic function is monotonically increasing) and attains its maximum at $\\sigma'(0) = \\tfrac{1}{4}$. \n",
    "\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "tags": [
     "hide-input"
    ]
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEGCAYAAACKB4k+AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3deXhV5b328e8v80wgk0CYZRAFBAIoHjTOWqnWoVapdS5Fa/V46tT2tMW21lZbW8fyUmtF8eCA1ipinWq0FVEI80yIEELCEAiZ5zzvH4kh0IAJJHvtnX1/ritXsrNWwp2Hndz7WaM55xAREQEI8TqAiIj4D5WCiIi0UCmIiEgLlYKIiLRQKYiISIswrwMci+TkZDdw4ECvY1BRUUFsbKzXMfyCxuKA7Oxsxo8f73UMv6DnxQH+MBbZ2dlFzrmUtpZZIB+SmpGR4ZYuXep1DLKyssjMzPQ6hl/QWBxgZgTy71dn0vPiAH8YCzPLds5ltLVMm49ERKSFSkFERFqoFEREpIVKQUREWvikFMzsGTPbbWZrDrPczOwxM8sxs1VmNs4XuURE5GC+mik8C1xwhOUXAkOb36YDf/JBJhEROYRPSsE59zGw7wirXAI855osBhLNrLcvsomIyAGhM2fO9Mk/dP/99ycC02bOnPlUG8tmAAtnzpyZ1/z4UuCzmTNnFhy6rplNv//++//f/fffP720tLTPpk2bSExMZMqUKaSkpLB48WKmTZvGqFGjuOeee/jwww+JjY0lIyODwYMH8+677zJjxgxGjBjBbbfdRnZ2NpGRkYwcOZKxY8fyyiuvcO+99zJkyBBuuukmNm/eTEhICIMGDSIzM5Onn36aBx54gAEDBjBt2jR27dpFZGQkEyZM4LTTTuPxxx/n8ccfp0+fPnzjG9+gurqa/Px8Tj/9dCZOnMhDDz3EnDlzSElJ4bzzziMiIoL169czdepUxo0bx8yZM3n99dc9/5nKysoYM2ZMh3+mm2++mcjIyG71Mx3t/9O6devIzMzsVj/T0f4/TZo0iVtuuaVb/UxH+/80f/58nnrqqS75mQYNP4kf3n0PC975J8WN0Zx52kSqEvoz9/V/cMdtt1AYezz/c8cP2J2XUzhz5szZbf2t9tnJa2Y2EFjgnDupjWVvAQ865/7d/PgD4B7nXPaRvqdOXvM/GosDdPLaAXpeHHCksWhodOyvrGVfRdNbcWUdxZW1lFTVsb+yjpKqOkqr6iit/vJ9PWXV9ZRV11FT39juDNt+O/WwJ6/5y2Uu8oF+rR6nA/8xSxARCUTVdQ3sLKlmZ2k1iwvq2fxxLrtKqykqr2FPeQ17ymrYW15LcWUtjYd5HREeavSIjiAhOowe0eEkxkTQr1cM8VHhJESFERsZRlxkGHFRTe9jIkKb3zd9HBMRSnREKFHhoUT89vBZ/aUU3gBuM7MXgUlAiXOu0ONMIiLtUl3XwPZ9leQ1v+UXV5FfXMmO/VUU7q9mb0XtwV+waj1R4SGkxkeREh/JoORYMgb2Ijk2gl6xEfSKi6RXTAQ9Y8PpGRNBYkw40eGhmFmX/yw+KQUzmwdkAslmlg/8HAgHcM7NAhYCXwNygErgBl/kEhFpL+ccu8tq2LyrnM27y8jdU0FuUTlf7KmgoKT6oHWjwkNI7xlD38RoRvVNpHePqOa3aPI2rmLqOVOIjwzzyR/5jvJJKTjnrv6K5Q74vi+yiIh8leq6BtYXlrK+sIz1haVs2FnKxp1llFbXt6wTHxXG4JQ4Jg1OYkBSDAOTYumfFEP/XjEkxUYc9g9+1o4QEqLCffWjdJi/bD4SEfFEfUMjm3aVs3x7MSvy9rN6Rwmbd5fT0LxxPz4yjBG947n45D4MS4vn+NQ4jk+NIyUu0i9f6R8rlYKIBJWq2gaWby/m8y/2sWTrPpbn7aeytgGAXrERjOrbg3NOSOOkvj04sU8C6T2ju+Uf/8NRKYhIt1bf0MjK/BI+ySnik5wiluUVU9fgMIMTjkvgm+PTGTegJ2P79aRfr+AqgLaoFESk29lbXkPWxj18uHE3H2/aQ2l1PWZwYp8EbjxtEKcMTmLcgJ70iPbfbfteUSmISLewfV8l76zdybtrd7F02z4aHaTER3L+iceROTyVyUOS6Bkb4XVMv6dSEJGAtbOkmgWrCnhzVSErt+8HYMRx8fzgrKGcOzKNkb0TCAkJ7s1BHaVSEJGAUlXbwNtrCpmfnc+nuXtxDk7qm8B9F47gwpOOY0BSrNcRA5pKQUQCwpodJbzw2TbeXFlIeU09/XpFc/tZQ7n45D4MSYnzOl63oVIQEb9VXdfAglWFPL94Gyu37ycqPISLRvXhmxnpTBzYS5uGuoBKQUT8zt7yGuYuzuP5xVspKq9lSEosP//6SC4bl64jhrqYSkFE/EZ+cSWzPtrCK0vzqalv5MzhKdw8ZTCThyQF/fkDvqJSEBHPbS2q4KmsHF5btgMzuGxsOjdPGcTQtHivowUdlYKIeKawpIpH39/MK9n5hIUY15wygO+dMZjePaK9jha0VAoi4nP7K2t58sMc5ny6DRx855QB3HrmEFLjo7yOFvRUCiLiM3UNjcxdvI0/vr+Zsuo6LhuXzn+fM5T0njFeR5NmKgUR8Ymsjbv5xYJ15O6p4L+OT+anU0cy/DjtM/A3KgUR6VK7Sqt5Ynk1S/+xhEHJsTx9bQZnn5Cqo4n8lEpBRLpEQ6Pj+U+38rt3N1FT18Bd5w3ju6cPJjIs1OtocgQqBRHpdLl7yrln/iqWbitmytBkvt67givPGup1LGkHlYKIdJqGRsdfP/mCh9/ZSGRYCI9cOYZLx/blo48+8jqatJNKQUQ6RcH+Ku58aQWffbGPc05I5deXjiI1QYeYBhqVgogcs7dXF3Lfa6upb2jk4StGc8X4dO1IDlAqBRE5atV1Ddz/5lrmfb6dMek9ePSqsQxM1v0MAplKQUSOyra9FdwydxnrCkuZccYQ/ufcYUSEhXgdS46RSkFEOuz9dbu48+UVhJjxzPUZnDUizetI0klUCiLSbo2Njkc/2MyjH2zmpL4J/Onb4+nXS5eo6E5UCiLSLpW19dz1ykoWrt7J5ePSeeDSk4gK14lo3Y1KQUS+UsH+Kr773FLWFZbyk6+dwM1TBunoom5KpSAiR7S2oIQb/rqEqtoGnrluAmeOSPU6knQhlYKIHNbHm/Zwy9xsekSHM/+WybqqaRBQKYhIm+Zn53Pfq6s4PjWOZ2+YyHE9dHZyMPDZQcVmdoGZbTSzHDO7r43lPczsTTNbaWZrzewGX2UTkYP9+eNc7nplJZMG9+LlGaeqEIKIT2YKZhYKPAmcC+QDS8zsDefcularfR9Y55z7upmlABvN7AXnXK0vMooIOOf4w/ubeeyDzVw0qjd/+NbJOiEtyPjqf3sikOOcy23+I/8icMkh6zgg3poOaYgD9gH1PsonEvQaGx2/WLCOxz7YzJUZ6Tx29VgVQhDy1T6FvsD2Vo/zgUmHrPME8AZQAMQD33LONR76jcxsOjAdIC0tjaysrK7I2yHl5eV+kcMfaCwOFihj0egcz62tJSu/nvMGhHFB0j7+9XHnXe5az4sD/H0sfFUKbR3Q7A55fD6wAjgLGAK8Z2b/cs6VHvRFzs0GZgNkZGS4zMzMzk/bQVlZWfhDDn+gsThYIIxFY6PjJ6+vISs/j1szh3D3+cM7/RwEPS8O8Pex8NXcMB/o1+pxOk0zgtZuAF5zTXKAL4ARPsonEpS+LIR5n+fx/TO7phAksPiqFJYAQ81skJlFAFfRtKmotTzgbAAzSwOGA7k+yicSdJw7uBDuOk+FID7afOScqzez24B3gFDgGefcWjOb0bx8FvBL4FkzW03T5qZ7nXNFvsgnEmycc/xywXrmfd60yUiFIF/y2clrzrmFwMJDPjer1ccFwHm+yiMSzP74/mae+eQLrp88UJuM5CA63kwkyPz541we/WAz3xyfzs+mjlQhyEFUCiJB5OWl23lg4XouGtWb31w+mpAQFYIcTKUgEiT+uWEXP3ptNVOGJvOHb51MqApB2qBSEAkCy/OKufWFZZzQO54/XTNeZyrLYemZIdLN5e4p58Znl5AaH8Vfr59IXKQujiyHp1IQ6cb2ltdw/V+XEGLGczdOJCU+0utI4uf0kkGkm6qua+B7z2ezq7SaedNPYWByrNeRJACoFES6Iecc9766iqXbinli2ljG9e/pdSQJENp8JNIN/fH9zfx9RQF3nz+cqaP7eB1HAohKQaSbWbCqgEc/2Mzl49K5NXOI13EkwKgURLqRNTtKuOuVlYwf0JNfX3aSzlaWDlMpiHQTe8pqmP7cUnrGRDDrmvFEhoV6HUkCkHY0i3QDtfWN3DI3m70VtcyfMVmHnspRUymIdAO/WLCWpduKeezqsYxK7+F1HAlg2nwkEuDmZ+czd3Ee008fzMVjdKSRHBuVgkgAW1tQwk/+tppTBydxz/nDvY4j3YBKQSRAlVTWMWNuNj1jInh82ljCQvXrLMdO+xREAlBjo+POl1ews6Sal753Kslx2rEsnUMvLUQC0KyPt/DPDbv56dSRuoSFdCqVgkiA+Sx3L79/dxNTR/fmO6cM8DqOdDMqBZEAUlReww/mLad/rxgevGyUzliWTqdSEAkQDY2OO19aQUlVHU9OG0d8VLjXkaQb0o5mkQAx66Mt/GtzEb+5bBQj+yR4HUe6Kc0URAJA9rZ9PPLeJr4+pg/fmtDP6zjSjakURPxcSWUdt89bQd/EaB64VFc+la6lzUcifsw5x32vrWJXaTXzb5lMgvYjSBfTTEHEj/3f53m8vWYnd58/nJP7JXodR4KASkHET+XsLuOXC9YxZWgy350y2Os4EiRUCiJ+qKa+gdvnrSAmIozff3MMISHajyC+oX0KIn7o9+9uYl1hKX++NoPUhCiv40gQ0UxBxM/8e3MRsz/O5duT+nPuyDSv40iQ8VkpmNkFZrbRzHLM7L7DrJNpZivMbK2ZfeSrbCL+Yn9lLT98ZQVDUmL534tGeh1HgpBPNh+ZWSjwJHAukA8sMbM3nHPrWq2TCDwFXOCcyzOzVF9kE/EXzjl+8voa9pbX8pfrJhAdEep1JAlCvpopTARynHO5zrla4EXgkkPWmQa85pzLA3DO7fZRNhG/8MbKAt5aVcid5w7jpL66z7J4w1c7mvsC21s9zgcmHbLOMCDczLKAeOBR59xzh34jM5sOTAdIS0sjKyurK/J2SHl5uV/k8Acai4O1dyz2VjXyv59UcXxiCCewnays/K4N5mN6Xhzg72Phq1Jo63g6d8jjMGA8cDYQDXxqZoudc5sO+iLnZgOzATIyMlxmZmbnp+2grKws/CGHP9BYHKw9Y9HY6Pj2058RElLLM989nf5JMV0fzMf0vDjA38fCV6WQD7S+ilc6UNDGOkXOuQqgwsw+BsYAmxDpxv66aCuf5u7lN5eN6paFIIHFV/sUlgBDzWyQmUUAVwFvHLLO34EpZhZmZjE0bV5a76N8Ip7I2V3OQ//YwNkjUnX1U/ELPpkpOOfqzew24B0gFHjGObfWzGY0L5/lnFtvZv8AVgGNwNPOuTW+yCfihfqGRn748gpiIkJ58HLdRU38g8/OaHbOLQQWHvK5WYc8fhh42FeZRLz0VNYWVuaX8OS0caTG66xl8Q86o1nEA6vzS3jsg81ccnIfLhrd2+s4Ii1UCiI+VlPfwA9fWUFSXAS/uPgkr+OIHEQXxBPxsT+8t5lNu8r56w0T6BGjm+aIf9FMQcSHluUVM/vjLXwrox9nDteVXMT/qBREfKSqtoG7Xl5J7x7R/O/UE7yOI9ImbT4S8ZGH39lIblEFL9w8iXjda1n8VIdnCmYW23zVUxFpp8+/2MdfF33BNaf057Tjk72OI3JYX1kKZhZiZtPM7C0z2w1sAAqb73nwsJkN7fqYIoGrsraeu+evJL1nND+6UJuNxL+1Z6bwITAE+BFwnHOun3MuFZgCLAZ+Y2bXdGFGkYD20D82sm1vJQ9fMYbYSG2xFf/WnmfoOc65ukM/6ZzbB7wKvGpm2kAq0oZPt+zl2UVbuX7yQE4ZnOR1HJGv9JWl8GUhmFkScCVQDawFVjvnqlqvIyIHu+fVlQxIiuGeC4Z7HUWkXTqyo/lvQArwa5quT1RiZhu6JJVIN5FfXMXDV4whJkKbjSQwdOSZGu+c+4WZXeacO8PMLgeO76pgIoFsUU4RADdMHsTEQb08TiPSfh2ZKVQ3v68xs2jn3KvA17ogk0hAq6ip555XVwFw9/nabCSBpSMzhd+ZWS/gJeAZM1tE072XRaSVB99ez479VQBER+iUHgks7Z4pOOdedc7tc849QtN9EfoBl3RZMpEA9ElOEXMX53HTaYO8jiJyVL5ypmBm5pxzrT/nnHv+q9YRCTblNfXcM38Vg5Njuev84fzU60AiR6FdJ6+Z2Q/MrH/rT5pZhJmdZWZzgOu6Jp5I4Pj1wvUUlFTx8DdHExWuzUYSmNqzT+EC4EZgnpkNBoqBaJoK5V3gD865FV0XUcT//WvzHv7vszymnz6Y8QN0tJEErvacvFYNPGVmKcCDQBJQ5Zzb39XhRAJBWXUd985fxZCUWP7n3GFexxE5Jh05+uhnQAzQC1hmZvNUDCLwwFvr2Vlazau3TNZmIwl4Hb10djXwDk1HHn1qZid3fiSRwPHRpj28uGQ7008fwtj+Pb2OI3LMOjJT2OCc+3nzx/PN7FlgFnBWp6cSCQAlVU2bjYamxvHf5+gK8tI9dGSmUGRm47984JzbRNO1kESC0v1vrmVPeQ2PXHmyNhtJt9GRmcLtwItmlg2sBkYDX3RJKhE/9966Xby2bAe3n3U8o9J7eB1HpNN05IzmlcDJwLzmT30IXN0VoUT8WXFFLT96bTUn9E7gtrO02Ui6lw5dz9c5VwO81fwmEpR++vc1lFTV8tyNE4kI6/BtzkX8mp7RIh3w5soCFqwq5I6zhzKyT4LXcUQ6nUpBpJ12lVbz07+v4eR+icw4Y4jXcUS6hEpBpB2cc9z76iqq6xr4/ZVjCAvVr450T3pmi7TDi0u2k7VxD/deMIIhKXFexxHpMj4rBTO7wMw2mlmOmd13hPUmmFmDmV3hq2wiR5K3t5JfLVjHqYOTuO7UgV7HEelSPikFMwsFngQuBEYCV5vZyMOs91uaLqUh4rmGRsedL68gxIzfXTmGkBDzOpJIl/LVTGEikOOcy3XO1QIv0vZd234AvArs9lEukSOa9dEWsrcV84tvnEjfxGiv44h0uQ6dp3AM+gLbWz3OBya1XsHM+gKX0nQtpQmH+0ZmNh2YDpCWlkZWVlZnZ+2w8vJyv8jhD7rTWGwrbeCRT6uZcFwoifs3k5WV0+Hv0V3G4lh1p+fFsfL3sfBVKbQ15z709p1/BO51zjWYHX6K7pybDcwGyMjIcJmZmZ2V8ahlZWXhDzn8QXcZi+q6Bh54/N8kxUXy5+mnkxgTcVTfpzuMRWfoLs+LzuDvY+GrUsin6XLbX0oHCg5ZJ4OmaysBJANfM7N659zrvokocsBv3t7A5t3lzLlx4lEXgkgg8lUpLAGGmtkgYAdwFTCt9QrOuUFfftx8We4FKgTxwocbd/Psoq3ccNpAzhimCwFLcPFJKTjn6s3sNpqOKgoFnnHOrTWzGc3LZ/kih8hXKSqv4e5XVjI8LZ57LxjhdRwRn/PVTAHn3EJg4SGfa7MMnHPX+yKTSGvOOe6Zv4rS6nrm3jxJ90iQoKQzmkWaPb94G//csJsfXTiCEcfpYncSnFQKIsD6wlJ+9dZ6zhiWwvWTB3odR8QzKgUJepW19fxg3nJ6RIfz+yvHcKRDokW6O5/tUxDxV/e/sY4te8p5/sZJJMdFeh1HxFOaKUhQe3NlAS8t3c4tZwzhv4Ymex1HxHMqBQla2/ZW8OPXVjOufyJ3njvM6zgifkGlIEGpuq6B7//fMszg0avGEq6b5ogA2qcgQepXb61jzY5S/nxtBv16xXgdR8Rv6OWRBJ03VhYwd3Ee008fzLkj07yOI+JXVAoSVHL3lPOjV1cxfkBP7j5/uNdxRPyOSkGCRkVNPTPmZhMRFsLjV2s/gkhbtE9BgoJzjntfXUXO7nKeu3ESfXQXNZE26aWSBIVnPtnKglWF/PC84TofQeQIVArS7X3+xT5+vXA9541M49bMIV7HEfFrKgXp1gpLqrj1hWUM6BXD73RdI5GvpH0K0m1V1zUw/blsqusamPfdSSREhXsdScTvqRSkW/ryhjlrCkr483cyGJoW73UkkYCgzUfSLc36KJc3VhZw13nDOUcnqIm0m0pBup331+3ioXc28PUxfbRjWaSDVArSrazZUcLtLy5nVN8ePHT5aO1YFukglYJ0GztLqrlpzhISo8N5+toMoiNCvY4kEnC0o1m6hYqaem6as4SKmgZemXEqqQlRXkcSCUiaKUjAq29o5AfzlrO+sJTHp43lhN4JXkcSCViaKUhAc87x47+t5p8bdvOrb5zEmcNTvY4kEtA0U5CA9of3NvHy0nxuP+t4rjllgNdxRAKeSkEC1vOLt/HYP3P4VkY/3WNZpJOoFCQgvbmygJ/9fQ1nj0jlgUtP0qGnIp1EpSAB54P1u7jzpRVMGNCLJ6aNI0w3yxHpNPptkoCyaEsRt7ywjJF9EvjL9ToXQaSzqRQkYGRvK+bmOUsZmBTDnBsmEq+rnop0OpWCBIRlecVc98znpMZHMvemSfSMjfA6kki35LNSMLMLzGyjmeWY2X1tLP+2ma1qfltkZmN8lU3827K8Yq79y+ckx0Uwb/opOltZpAv5pBTMLBR4ErgQGAlcbWYjD1ntC+AM59xo4JfAbF9kE/92aCH07hHtdSSRbs1XM4WJQI5zLtc5Vwu8CFzSegXn3CLnXHHzw8VAuo+yiZ9atKWI7zz9mQpBxId8dZmLvsD2Vo/zgUlHWP8m4O22FpjZdGA6QFpaGllZWZ0U8eiVl5f7RQ5/0FljsWJ3PU+sqCEtxrhjlGPj8s/YeOzxfE7Piyb6HTnA38fCV6XQ1plFrs0Vzc6kqRT+q63lzrnZNG9aysjIcJmZmZ0U8ehlZWXhDzn8QWeMxRsrC3ji3RWM7NODOTdMDOidynpeNNHvyAH+Pha+KoV8oF+rx+lAwaErmdlo4GngQufcXh9lEz8yZ9FWZr65lgkDe/GX6zJ02KmIj/mqFJYAQ81sELADuAqY1noFM+sPvAZ8xzm3yUe5xE80Njp++48N/L+PczlvZBqPXT2WqHCdmCbiaz4pBedcvZndBrwDhALPOOfWmtmM5uWzgJ8BScBTzdexqXfOZfgin3irpr6Bu15ZxZsrC7j21AH8/OsnEhqiaxmJeMFn91Nwzi0EFh7yuVmtPr4ZuNlXecQ/FJXXcMvcbJZsLea+C0fwvdMH6+J2Ih7STXbEM+sLS7l5zlKKymt4/OqxfH1MH68jiQQ9lYJ44p21O7nzpRXER4XxyoxTGZ2e6HUkEUGlID5W39DII+9t4qmsLYxJ78HsazNI02UrRPyGSkF8pqi8htvnLWfRlr1cNaEfMy8+UUcYifgZlYL4xGe5e7njxRUUV9by0BWjuTKj31d/kYj4nEpBulR9QyOPfrCZJz/MoX+vGF67dTIn9unhdSwROQyVgnSZ7fsquePF5SzL2883x6cz8+ITiY3UU07En+k3VDpdY6Pjhc/zeHDhekJDjMeuHsvFOtxUJCCoFKRT7als5Jq/fMaiLXuZMjSZ31w+mr6JuuS1SKBQKUinqG9o5NlFW3n4kyrCw+p48LJRXDWhn85OFgkwKgU5Zsvzivnx39awvrCU0SmhPHXjFNJ7xngdS0SOgkpBjtqeshp+/+5GXlq6ndT4SP707XFEFW1QIYgEMJWCdFhNfQPPfrKVx/+ZQ3VdAzdMHsSd5w4lPiqcrKxAvD+aiHxJpSDt1tDoeGPlDh55bxPb91Vx9ohUfnzRCQxJifM6moh0EpWCfCXnHO+v383v3tnIxl1ljOydwHM3juL0YSleRxORTqZSkMNqbHS8u24XT3y4mTU7ShmUHMvjV4/lolG9CdFNcES6JZWC/Ie6hkYWrCpgVlYuG3eVMTAphoeuGM2lY/sSHhridTwR6UIqBWlRUlnHvCV5PPvJVnaWVjM0NY5HrzqZi0b1JkxlIBIUVArCmh0lvPDZNl5fXkBVXQOThyTx4GWjOGNYijYTiQQZlUKQKquu461Vhby4ZDsrtu8nKjyEi8f04brJA3UVU5EgplIIIvUNjXyau5fXlu3g7TWFVNc1MiQllp9NHcnl49LpERPudUQR8ZhKoZtrbHRk5xXz5soCFq4upKi8lvioMC4fl84V49M5uV+irk8kIi1UCt1QTX0Dn27Zyztrd/Heul0UldcQGRbC2SekcvGYPmQOT9VtMEWkTSqFbiK/uJKPNu3hww17WLSliMraBmIjQskckcp5I9M4+4Q04nSDGxH5CvorEaCKymv4/It9/DuniEU5RWzdWwlAes9orhifzpnDUzl1SJJmBCLSISqFANDY6MgtqmB5XjHZ24r5fOs+cvdUABAXGcYpg3tx7akDOX1YMkNS4rSPQESOmkrBzzjn2LG/ijU7Sli9o4TVO0pZuX0/JVV1ACREhTFhYC+uzOjHxEG9GN23h04sE5FOo1LwUEllHZt3l7F5dzkbd5axvrCUDTvLWgogNMQYmhrH10Ydx9h+PRnbP5EhKXE6oUxEuoxKoYtV1taTt6+SvL2VbN1bwRdFFeTuqSC3qII9ZTUt68VGhDL8uHimju7NiN4JjOrbgxHHxWufgIj4lErhGDQ2OvZV1rK1pIF31+6ksKSagpIqdhRXkV9cxY79VQf94QdIio1gUHIsmcNSGJoWx9DUeI5PjaNvYrRmACLiOZXCIeoaGimurKW4oo59FbXsrahhX0UtReW17Cmroai8hj1lNewurWZ3WQ31ja7pCz/NBiA81OiTGE16z2jOHJ5C/14x9E+KZUCvGAYkxZAYE+HhTycicmQ+KwUzuwB4FAgFnnbO/eaQ5da8/GtAJXC9c25ZR/6NuoZGKmsaqKitp7K2nrLqeipqGiivqaesuo7ymnrKq+spra6jtKrpfUlVq7fKOspq6pAJQToAAAVISURBVA+Tv+lVfnJcJCnxkQxJSSYtIZK0hCj25OVw7mkZ9E6MIjk2Uq/4RSRg+aQUzCwUeBI4F8gHlpjZG865da1WuxAY2vw2CfhT8/vD2rSrjFN+/QGVtfVU1zVS29DYrjzR4aEkRIeREBVOQnQ4xyVEMTwtnoTocHrFRtAzJpzEmAiSYiNIioskKS6CxOjwwx7lk1W7lTH9Etv1b4uI+DNfzRQmAjnOuVwAM3sRuARoXQqXAM855xyw2MwSzay3c67wcN80KjyU04clExMRRnREKLERocREhBETEUpMZBhxkaHERoQRG9lUAPFRTR9HhOkQThGRtviqFPoC21s9zuc/ZwFtrdMXOKgUzGw6MP3Lxw9/8+RODSrSmXQioQQaX5VCW78Z7ijWwTk3G5gNkJGR4ZYuXXrs6Y5RVlYWmZmZXsfwCxqLA8yMpomv6HlxgD+MxZFerPhqO0o+0K/V43Sg4CjWERGRLuSrUlgCDDWzQWYWAVwFvHHIOm8A11qTU4CSI+1PEBGRzueTzUfOuXozuw14h6ZDUp9xzq01sxnNy2cBC2k6HDWHpkNSb/BFNhEROcBn5yk45xbS9Ie/9edmtfrYAd/3VR4REflPOjZTRERaqBRERKSFSkFERFpYIB9HbWZ7gG1e5wCSgSKvQ/gJjcUBGosDNBYH+MNYDHDOpbS1IKBLwV+Y2VLnXIbXOfyBxuIAjcUBGosD/H0stPlIRERaqBRERKSFSqFzzPY6gB/RWBygsThAY3GAX4+F9imIiEgLzRRERKSFSkFERFqoFDqRmd1lZs7Mkr3O4hUze9jMNpjZKjP7m5kF3X1KzewCM9toZjlmdp/XebxiZv3M7EMzW29ma83sDq8zec3MQs1suZkt8DrL4agUOomZ9aPpHtR5Xmfx2HvASc650cAm4Ece5/GpVvcjvxAYCVxtZiO9TeWZeuCHzrkTgFOA7wfxWHzpDmC91yGORKXQef4A3EMbd4sLJs65d51z9c0PF9N0s6Rg0nI/cudcLfDl/ciDjnOu0Dm3rPnjMpr+GPb1NpV3zCwduAh42ussR6JS6ARmdjGwwzm30ussfuZG4G2vQ/jY4e41HtTMbCAwFvjM2ySe+iNNLxwbvQ5yJD67n0KgM7P3gePaWPQT4MfAeb5N5J0jjYVz7u/N6/yEps0HL/gymx9o173Gg4mZxQGvAv/tnCv1Oo8XzGwqsNs5l21mmV7nORKVQjs5585p6/NmNgoYBKxsvhl2OrDMzCY653b6MKLPHG4svmRm1wFTgbNd8J0Io3uNt2Jm4TQVwgvOude8zuOh04CLzexrQBSQYGZznXPXeJzrP+jktU5mZluBDOec11dB9ISZXQA8ApzhnNvjdR5fM7Mwmnawnw3soOn+5NOcc2s9DeYBa3qVNAfY55z7b6/z+IvmmcJdzrmpXmdpi/YpSGd7AogH3jOzFWY266u+oDtp3sn+5f3I1wMvB2MhNDsN+A5wVvNzYUXzK2XxY5opiIhIC80URESkhUpBRERaqBRERKSFSkFERFqoFEREpIVKQUREWqgURESkhUpBpJOZ2RVmttjMVprZv80sxetMIu2lk9dEOpmZJTnn9jZ//HOgyDn3pMexRNpFMwWRzne9mX1uZiuBW4FqrwOJtJeukirSiczsWpputHOWc67czD4GgvXaRxKANFMQ6VyjgEXNhXA5MBlY7XEmkXZTKYh0rjnA7Wb2L2AYkOucq/A4k0i7aUeziIi00ExBRERaqBRERKSFSkFERFqoFEREpIVKQUREWqgURESkhUpBRERa/H8qV5y1tMU7jgAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import torch\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "aa = torch.linspace(-5, 5, 100)\n",
    "plt.plot(aa, torch.sigmoid(aa))\n",
    "plt.axhline(0, color='k', lw=1)\n",
    "plt.axhline(1, color='k', ls=':', lw=1)\n",
    "plt.axvline(0, color='k', lw=1)\n",
    "plt.grid(\"True\")\n",
    "plt.xlim(-5, 5)\n",
    "plt.xlabel(r\"$a$\")\n",
    "_ = plt.ylabel(r\"$\\sigma(a)$\")"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Optimization\n",
    "\n",
    "Our goal is to estimate the weights $\\mathbf{w}_{\\mathsf{MLE}}$ that maximize the log likelihood of the training data, or equivalently minimize the negative log likelihood. Unlike most of the problems we've encountered thus far, we won't have a closed form solution for the weights, even if we try to do coordinate ascent. Instead we'll have to turn to more general optimization strategies like gradient descent and Newton's method. This is a good opportunity to introduce a few of these tools and some of the key concepts."
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Setup\n",
    "\n",
    "Let $\\mathbf{X} \\in \\mathbb{R}^{N \\times P}$ denote the matrix of inputs (with rows $\\mathbf{x}_n^\\top$) and $\\mathbf{y} = (y_1, \\ldots, y_N) \\in \\{0,1\\}^N$ denote the vector of output labels. The **negative** log likelihood is,\n",
    "\n",
    "$$\n",
    "\\begin{align*}\n",
    "\\mathcal{L}(\\mathbf{w}) &= -\\log p(\\mathbf{y} \\mid \\mathbf{w}, \\mathbf{X}) \\\\\n",
    "&= -\\sum_{n=1}^N \\log p(y_n \\mid \\mathbf{w}, \\mathbf{x}_n) \\\\\n",
    "&= -\\sum_{n=1}^N \\log \\mathrm{Bern}(y_n; \\sigma(\\mathbf{w}^\\top \\mathbf{x}_n)) \\\\\n",
    "&= -\\sum_{n=1}^N \\left(y_n \\mathbf{w}^\\top \\mathbf{x}_n - \\log(1 + e^{\\mathbf{w}^\\top \\mathbf{x}_n}) \\right)\n",
    "\\end{align*}\n",
    "$$"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### The gradient with respect to the weights\n",
    "The gradient with respect to $\\mathbf{w}$ is,\n",
    "\n",
    "$$\n",
    "\\begin{align*}\n",
    "\\nabla \\mathcal{L}(\\mathbf{w}) \n",
    "&= -\\sum_{n=1}^N \\left( y_n \\mathbf{x}_n - \\frac{e^{\\mathbf{w}^\\top \\mathbf{x}_n}}{1 + e^{\\mathbf{w}^\\top \\mathbf{x}_n}} \\mathbf{x}_n\\right) \\\\\n",
    "&= -\\sum_{n=1}^N \\big(y_n - \\sigma(\\mathbf{w}^\\top \\mathbf{x}_n)\\big) \\, \\mathbf{x}_n\n",
    "\\end{align*}\n",
    "$$\n",
    "\n",
    "Unfortunately, this is a nonlinear function of $\\mathbf{w}$ (due to the logistic function), and when we set to zero and try to solve for the weights, we find there is no closed-form solution.\n"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### The negative log likelihood is convex\n",
    "\n",
    "While there may not be a closed-form solution, the problem is not necessarily all that hard to solve. It turns out the negative log likelihood is a **convex function** of the weights &mdash; i.e., it looks like an upward facing bowl &mdash; so we can solve it with off-the-shelf optimization tools. \n",
    "\n",
    "To show that the objective function is convex, it suffices to show that it is twice-differentiable and its [**Hessian**](https://en.wikipedia.org/wiki/Hessian_matrix) (the matrix of second-order partial derivatives) is positive semi-definite (has eigenvalues $\\geq 0$).\n",
    "\n",
    "The Hessian of the negative log likelihood is,\n",
    "\n",
    "$$\n",
    "\\begin{align*}\n",
    "\\nabla^2 \\mathcal{L}(\\mathbf{w})\n",
    "&= \\sum_{n=1}^N \\sigma'(\\mathbf{w}^\\top \\mathbf{x}_n) \\mathbf{x}_n \\mathbf{x}_n^\\top \n",
    "\\end{align*}\n",
    "$$\n",
    "\n",
    "where $\\sigma'(\\mathbf{w}^\\top \\mathbf{x}_n)$ is the derivative of the logistic function (see above) evaluated at $\\mathbf{w}^\\top \\mathbf{x}_n$.\n",
    "\n",
    "Since this a sum of outer products ($\\mathbf{x}_n \\mathbf{x}_n^\\top$) with positive coefficients ($\\sigma(\\mathbf{w}^\\top \\mathbf{x}_n)$), the Hessian is positive semi-definite. \n",
    "\n",
    "```{admonition} Matrix derivatives\n",
    ":class: tip\n",
    "\n",
    "It takes some practice to become familiar with the rules of matrix calculus. I recommend the first chapters of [The Matrix Cookbook](https://www.math.uwaterloo.ca/~hwolkowi/matrixcookbook.pdf) {cite}`petersen2008matrix` for an introduction. \n",
    "```"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Gradient descent\n",
    "\n",
    "Since the the negative log likelihood is convex (equivalently, the log likelihood is concave), we have a host of tools at our disposal for maximum likelihood estimation. We don't need CVXPy to solve this problem (like we did for the previous chapter). Here, we can simply perform **gradient descent**. \n",
    "\n",
    "Let $\\mathbf{w}_0$ denote our initial setting of the weights. Gradient descent is an iterative algorithm that produces a sequence of weights $\\mathbf{w}_0, \\mathbf{w}_1, \\ldots$ that (under certain conditions) converges to a local optimum of the objective. Since the objective is convex, all local optima are global optima. The idea is straightforward, on each iteration we update the weights by taking a step in the direction of the gradient,\n",
    "\n",
    "$$\n",
    "\\begin{align*}\n",
    "\\mathbf{w}_{i+1} &= \\mathbf{w}_i - \\alpha_i \\nabla \\mathcal{L}(\\mathbf{w}_i)\n",
    "\\end{align*}\n",
    "$$\n",
    "\n",
    "where $\\alpha_i \\in \\mathbb{R}_+$ is the **learning rate** (aka step size) on iteration $i$, and $\\nabla \\mathcal{L}(\\mathbf{w}_i)$ is the gradient of the objective evaluated at the current weights $\\mathbf{w}_i$. \n",
    "\n",
    "```{admonition} Convergence of gradient descent\n",
    ":class: dropdown\n",
    "\n",
    "When $\\mathcal{L}$ is convex and $\\nabla \\mathcal{L}$ is [**Lipschitz continuous**](https://en.wikipedia.org/wiki/Lipschitz_continuity), and when the learning rates satisfy certain conditions (e.g., they are set by a [**backtracking line search**](https://en.wikipedia.org/wiki/Backtracking_line_search)), then gradient descent will converge to a local optimum. When the objective is convex, all local optima are global optima, so then gradient descent is guaranteed to find a global solution.\n",
    "```"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Newton's method\n",
    "\n",
    "Gradient descent uses first-order information (i.e., the gradient of the objective at the current weights) to determine the descent direction. We can obtain faster convergence rates using **second-order** information (i.e., the Hessian of the objective). \n",
    "\n",
    "The idea is to minimize a quadratic approximation of the objective given by a Taylor approximation around the current weights,\n",
    "\n",
    "$$\n",
    "\\begin{align*}\n",
    "\\mathcal{L}(\\mathbf{w}) \n",
    "&\\approx \\mathcal{L}(\\mathbf{w}_i) + (\\mathbf{w} - \\mathbf{w}_i)^\\top \\nabla \\mathcal{L}(\\mathbf{w}_i) \n",
    "+ \\frac{1}{2} (\\mathbf{w} - \\mathbf{w}_i)^\\top \\nabla^2 \\mathcal{L}(\\mathbf{w}_i) (\\mathbf{w} - \\mathbf{w}_i).\n",
    "\\end{align*}\n",
    "$$\n",
    "\n",
    "The minimum of this quadratic approximation has a closed form solution,\n",
    "\n",
    "$$\n",
    "\\mathbf{w}_{i+1} = \\mathbf{w}_i + \\nabla^2 \\mathcal{L}(\\mathbf{w}_i)^{-1} \\nabla \\mathcal{L}(\\mathbf{w}_i).\n",
    "$$\n",
    "\n",
    "Here, the descent direction is given by the inverse-Hessian times the gradient, $\\nabla^2 \\mathcal{L}(\\mathbf{w}_i)^{-1} \\nabla \\mathcal{L}(\\mathbf{w}_i)$. \n",
    "\n",
    "```{warning}\n",
    "Note that Newton's method assumes that the Hessian is invertible, which is almost surely the case for logistic regression with many data points. We can ensure invertibility by adding a multivariate normal prior on the weights, as we will introduce below.\n",
    "```\n",
    "\n",
    "Newton's method can be unstable in practice. A simple fix is to use the same descent direction, but that to vary the step size $\\alpha_i$,\n",
    "\n",
    "$$\n",
    "\\mathbf{w}_{i+1} = \\mathbf{w}_i + \\alpha_i \\nabla^2 \\mathcal{L}(\\mathbf{w}_i)^{-1} \\nabla \\mathcal{L}(\\mathbf{w}_i).\n",
    "$$\n",
    "\n",
    "For example, the step-size can be set to $\\alpha_i < 1$ to implement _damped_ Newton's method, or it can be determined by a backtracking line search. "
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Iteratively reweighted least squares (IRLS)\n",
    "\n",
    "The weight updates simplify nicely when we substitute in the form of the gradient and Hessian for logistic regression. Note that they can be written in matrix form as,\n",
    "\n",
    "$$\n",
    "\\begin{align*}\n",
    "\\nabla \\mathcal{L}(\\mathbf{w}) &= -\\mathbf{X}^\\top \\big(\\mathbf{y} - \\sigma(\\mathbf{X} \\mathbf{w}) \\big) \\\\\n",
    "\\nabla^2 \\mathcal{L}(\\mathbf{w}) &= \\mathbf{X}^\\top \\mathbf{S} \\mathbf{X}\n",
    "\\end{align*}\n",
    "$$\n",
    "\n",
    "where\n",
    "\n",
    "$$\n",
    "\\mathbf{S} = \\mathrm{diag}\\left([\\sigma'(\\mathbf{w}^\\top \\mathbf{x}_1), \\ldots, \\sigma'(\\mathbf{w}^\\top \\mathbf{x}_N)] \\right)\n",
    "$$\n",
    "\n",
    "is a diagonal **scaling (aka weighting) matrix**. Note that the scale factors are all positive since $\\sigma'(a) > 0$. \n",
    "\n",
    "Substituting these forms in and rearranging yield,\n",
    "\n",
    "$$\n",
    "\\begin{align*}\n",
    "\\mathbf{w}_{i+1} \n",
    "&= \\mathbf{w}_i + \\big(\\mathbf{X}^\\top \\mathbf{S} \\mathbf{X}\\big)^{-1} \\mathbf{X}^\\top \\big( \\mathbf{y} - \\sigma(\\mathbf{X} \\mathbf{w}) \\big) \\\\\n",
    "&= \\big(\\mathbf{X}^\\top \\mathbf{S} \\mathbf{X}\\big)^{-1} \\mathbf{X}^\\top \\mathbf{S} \\tilde{\\mathbf{y}},\n",
    "\\end{align*} \n",
    "$$\n",
    "\n",
    "where \n",
    "\n",
    "$$\n",
    "\\tilde{\\mathbf{y}} = \\mathbf{X} \\mathbf{w}_i + \\mathbf{S}^{-1} \\big(\\mathbf{y} - \\sigma(\\mathbf{X} \\mathbf{w}_i) \\big).\n",
    "$$\n",
    "\n",
    "In other words, the standard Newton method update can be viewed as the solution to a [**weighted least squares**](https://en.wikipedia.org/wiki/Weighted_least_squares) problem with weights $\\mathbf{S}$ and targets $\\tilde{\\mathbf{y}}$ that depend on the current weights $\\mathbf{w}_i$. Viewed this way, we see that Newton's method for logistic regression is equivalent to an algorithm called [**iteratively reweighted least squares (IRLS)**](https://en.wikipedia.org/wiki/Iteratively_reweighted_least_squares). "
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Computational complexity\n",
    "\n",
    "```{admonition} Exercise\n",
    ":class: tip\n",
    "\n",
    "Show that the computational complexity of gradient descent is $\\mathcal{O}(NP)$ whereas the complexity of Newton's method is $\\mathcal{O}(NP^2 + P^3)$. \n",
    "```"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Scaling up\n",
    "\n",
    "Though it converges faster, Newton's method quickly becomes intractable for large $N$ and $P$. For large $N$, even gradient descent becomes costly. There are a few ways to speed up computation:\n",
    "\n",
    "- **Quasi-Newton methods** like [**BFGS**](https://en.wikipedia.org/wiki/BFGS_method) replace the exact Hessian with an approximation and side-step the explicit matrix inversion.\n",
    "- [**Stochastic gradient descent (SGD)**](https://en.wikipedia.org/wiki/Stochastic_gradient_descent) uses subsets of data (a.k.a. _minibatches_) to approximate the gradient. Under fairly general conditions, it converges to a local optimum.\n",
    "- **Momentum** is often used in conjunction with SGD to keep the descent direction from changing too rapidly. This can address some of the limitations of regular gradient descent as well, e.g., where the updates overshoot in poorly conditioned problems. Related methods like Nesterov's accelerated gradient (see {cite:t}`sutskever2013importance`) can achieve second-order convergence rates using first-order information (under certain conditions).\n",
    "- Still, SGD (with momentum) requires **setting a learning rate**. Modern machine learning packages like [`torch.optim`](https://pytorch.org/docs/stable/optim.html) implement a number of optimizers that automatically tune the learning rates, like **AdaGrad** {cite}`duchi2011adaptive`, **RMSProp** {cite}`rmsprop`, and **Adam** {cite}`kingma2014adam`. "
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Pose tracking with convolutional neural networks (CNNs)\n",
    "\n",
    "Remember we started by treating each data point as patch $\\mathbf{x}_n$ and a binary label $y_n$ specifying whether a specific key point (e.g. \"left paw\") is present. Of course, in practice we want to classify all the patches in an image in parallel, and we want to predict more than one type of key point.\n",
    "\n",
    "Let $\\mathbf{X}_n \\in \\mathbb{R}^{P_H \\times P_W}$ denote the $n$-th image and $\\mathbf{Y}_{n,k} \\in \\{0,1\\}^{P_H \\times P_W}$ denote the binary mask of where in the $k$-th key point is in the $n$-th image. Both are $P_H$ pixels in height and $P_W$ pixels wide. Assume the patches are $P_h \\times P_w$ in size, with $P_h < P_H$ and $P_w < P_W$.  Finally, ket $\\mathbf{W}_k \\in \\mathbb{R}^{P_h \\times P_w}$ denote the weights for the $k$-th key point. \n",
    "\n",
    "We can think of each image as having $P_H \\cdot P_W$ patches and corresponding labels, one centered on each pixel. In our simple logistic regression model, each patch's label is modeled as a conditionally independent Bernoulli random variable,\n",
    "\n",
    "$$\n",
    "\\begin{align*}\n",
    "p(\\mathbf{Y} \\mid \\mathbf{X}, \\mathbf{W}) \n",
    "&= \\prod_{n=1}^N \\prod_{k=1}^K \\prod_{i=1}^{P_H} \\prod_{j=1}^{P_W} \\mathrm{Bern}(y_{n,k,i,j}; \\sigma(a_{n,k,i,j}))\n",
    "\\end{align*}\n",
    "$$\n",
    "\n",
    "where $a_{n,k,i,j} \\in \\mathbb{R}$ is the **activation** for that specific image, key point, and pixel. The activations are given by a **2D cross-correlation**,\n",
    "\n",
    "$$\n",
    "\\begin{align*}\n",
    "a_{n,k,i,j} \n",
    "&= \\sum_{d=1}^{P_h} \\sum_{d'=1}^{P_w} w_{k,d,d'} x_{n,i+d-\\frac{P_h}{2},j+d'-\\frac{P_w}{2}} \\\\\n",
    "&= [\\mathbf{X_n} \\star \\mathbf{W}_k]_{i,j}.\n",
    "\\end{align*}\n",
    "$$\n",
    "\n",
    "In fact, the activations for an entire batch of $N$ images and $K$ key points (i.e., output channels) can be computed in a single call to [`F.conv2d`](https://pytorch.org/docs/stable/generated/torch.nn.functional.conv2d.html), with the appropriate padding. In lab, we'll make use of the [`torch.nn.Conv2d`](https://pytorch.org/docs/stable/generated/torch.nn.Conv2d.html) class, which encapsulates the weights of a 2D convolution layer and makes it easy to train such models.\n"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Convolutional neural networks\n",
    "\n",
    "Framed this way, we can view the logistic regression model as a **one-layer convolutional neural network (CNN)**. This view also suggests an obvious direction for improvement. The activations of the logistic regression model are **linear functions** of the pixels. In practice, a good key point detector may need **nonlinear features** of the images. Moreover, the features necessary to predict one keypoint (e.g., the left paw) may be similar to those needed for another (e.g., the right). \n",
    "\n",
    "Convolutional neural networks allow both nonlinear feature learning and feature sharing between outputs. The idea is straightforward: stack multiple convolutional layers on top of each other, feeding the output of one as the input to the next. \n",
    "\n",
    "**Residual networks (ResNets)** enable very deep CNNs to be stably trained by adding **skip connections** whereby the input is fed straight to the output of a layer, thereby allowing the convolution to capture the difference (i.e., residual) between the input and output.\n",
    "\n",
    "These notes will not comprehensively cover CNNs and ResNets. Instead, please consult the many great online resources, like {cite:t}`Goodfellow-et-al-2016` and the PyTorch [tutorials](https://pytorch.org/tutorials/beginner/blitz/cifar10_tutorial.html)."
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Working with minimal labeled data.\n",
    "\n",
    "While deep neural networks like CNNs and ResNets have the capacity to learn nonlinear mappings from images to key point predictions, they also have many free parameters to train.\n",
    "\n",
    "For modern machine learning tasks and benchmarks, these models are trained on millions of labeled images, like the ImageNet dataset for image classification. \n",
    "\n",
    "Unfortunately, we don't (yet) have comparable datasets of labeled images of animals in the lab (however, see {cite:t}`marshall2021continuous`) Moreover, environments change from one lab to the next, making it challenging to collect one dataset that covers all the bases.\n",
    "\n",
    "Ideally, we would like to train a pose tracking model for a specific lab environment using minimal labeled data.\n",
    "\n",
    "Modern methods for markerless animal pose tracking take two approaches:\n",
    "1. Use a **lightweight architecture** with fewer free parameters, like the modular UNet in SLEAP {cite}`pereira2022sleap`\n",
    "2. Use **transfer learning** to adapt the weights of a pre-trained network to the problem of animal pose tracking {cite}`Mathis2018-uz`."
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Transfer learning\n",
    "\n",
    "**Transfer learning** is when you take a model trained for one task (e.g. a ResNet trained for image classification with ImageNet) and adapt it to another task (e.g. classifying key points in behavioral video). If the tasks are similar enough, then it may require new minimal training data to adapt the pre-trained network. \n",
    "\n",
    "For example, DeepLabCut {cite}`Mathis2018-uz` originally converted a ResNet-50 (a very deep CNN) trained on ImageNet by rerouting the output of a middle layer to make predictions about key point locations in animal videos. \n",
    "We can think of the activations at the middle layer as a highly nonlinear **feature bank** that represents the original image in a form that may be more amenable to key point prediction with simple models.\n",
    "\n",
    "Finally, using a small number (100s) of hand-labeled video frames, DLC can train the mapping from intermediate features to key point predictions. It can also **fine-tune** the weights of the ResNet to obtain very good accuracy in key point prediction.\n",
    "\n",
    "There are, of course, many details and design questions to answer. Which layer of the ResNet do you reroute? How complex is the mapping from features to predictions? How big should you make the prediction targets? We will work through these concerns in Lab 3."
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Structured prediction\n",
    "\n",
    "The output of these CNNs is a map of key point probabilities for each image, key point, and pixel. However, the key points are obviously not independent &mdash; knowing where the left elbow is tells you a lot about where the left shoulder is likely to be found. Ideally, we would like to make a unified prediction about the collection of key point locations. \n",
    "\n",
    "One way to tackle this problem is with **structured prediction** {cite}`felzenszwalb2005pictorial,insafutdinov2016deepercut`. Let $\\ell_{n,k} \\in [0,P_h] \\times [0,P_w]$ be a random variable denoting the location of key point $k$ in image $n$ (in pixel coordinates). We can represent a joint distribution over the collection of key points $\\boldsymbol{\\ell}_n = (\\ell_{n,1}, \\ldots, \\ell_{n,K})$ as a **factor graph**,\n",
    "\n",
    "$$\n",
    "\\begin{align*}\n",
    "p(\\boldsymbol{\\ell}_n) &\\propto \n",
    "\\prod_{k=1}^K e^{f_k(\\ell_{n,k}; \\mathbf{x}_n)} \n",
    "\\prod_{i=1}^K \\prod_{j=1}^K e^{g_{ij}(\\ell_{n,i}, \\ell_{n,j})},\n",
    "\\end{align*}\n",
    "$$\n",
    "\n",
    "where $f_k$ are called **unary potentials** and $g_{ij}$ are called **pairwise potentials**. The unary potentials are functions that specify how likely key point $k$ is to be found at a location $\\ell_{n,k}$ given the observed image $\\mathbf{x}_n$; the pairwise potentials specify how likely the joint configuration of $(\\ell_{n,i}, \\ell_{n,j})$ is. For example, the pairwise potentials could rule out configurations in which the left elbow is meters away from the left shoulder.\n",
    "\n",
    "Note that this joint distribution is **unnomrmalized**. That isn't a problem for maximuim a posteriori (MAP) estimation of the key point locations, but it does make learning the parameters of the potential functions a bit harder. We'll come back to this later in the course."
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3D Triangulation \n",
    "\n",
    "Finally, we have so far focused on finding key points in 2D images, but there is growing interest in tracking 3D pose from multiple camera views. \n",
    "\n",
    "One way to solve this problem is by training a CNN that synthesizes multiple 2D views into a 3D representation {cite}`dunn2021geometric`. Another is to use **triangulation** to infer a 3D key point location from multiple 2D views {cite}`nath2019using,karashchuk2021anipose,zhang2021animal`.\n",
    "\n",
    ":::{figure} images/07_pose_tracking/triangulation.png\n",
    "---\n",
    "name: triangulation\n",
    "---\n",
    "Projective geometry makes far away objects appear smaller. Image modified from [Wikipedia](https://en.wikipedia.org/wiki/Triangulation_(computer_vision)).\n",
    ":::\n",
    "\n",
    "GIMBAL {cite}`zhang2021animal` combines structured prediction and triangulation to infer 3D pose from multiple 2D key point estimates. Changing notation slightly, let $\\mathbf{y}_{c,k} \\in \\mathbb{R}^2$ denote the estimated location of key point $k$ in the coordinates of camera $c$, and let $\\mathbf{x}_k \\in \\mathbb{R}^3$ denote the 3D location in \"world coordinates.\" [**Projective geometry**](https://en.wikipedia.org/wiki/Projective_geometry) tells how to map 3D locations to 2D images,\n",
    "\n",
    "$$\n",
    "\\mathbf{y}_{c,k} \\approx \\frac{1}{w} \\begin{bmatrix} u \\\\ v \\end{bmatrix} \n",
    "$$\n",
    "\n",
    "where\n",
    "\n",
    "$$\n",
    "\\begin{bmatrix} u \\\\ v \\\\ w \\end{bmatrix} = \\mathbf{A}_c \\mathbf{x}_k + \\mathbf{b}_c\n",
    "$$\n",
    "\n",
    "are the 3D coordinates after affine transformation (parameterized by $\\mathbf{A}_c$ and $\\mathbf{b}_c$) into the frame of camera $c$. Note that the transformation is **nonlinear** due to the scaling by $w^{-1}$. \n",
    "\n",
    "GIMBAL combines this nonlinear camera model with a structured prior distribution on 3D key point locations $\\mathbf{x} = (\\mathbf{x}_1, \\ldots, \\mathbf{x}_K)$ obtained from a hierarchical von Mises-Fisher-Gaussian distribution. "
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Conclusion\n",
    "\n",
    "This chapter introduced the basics of modern animal pose tracking methods using supervised learning, deep convolutional neural networks, transfer learning, and structured prediction. Many of these techniques will reappear later in the course as we talk about encoding and decoding neural spike trains. In lab, you will practice implementing these techniques and building a simple pose tracker."
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Further reading\n",
    "\n",
    "As always, there is much more to say on each of the topics covered in this chapter. There is growing work on **semi-supervised** pose tracking to leverage the wealth of unlabeled image frames {cite}`wu2020deep,whiteway2021partitioning`. There is still work to be done on multi-animal pose tracking {cite}`lauer2022multi,pereira2022sleap`, and recent machine learning challenges like [MABe](https://sites.google.com/view/mabe22/home) have sought to make progress on these problems. Finally, experimentalists are starting to use these kinds of techniques for **real-time** pose tracking to guide closed-loop optogenetic stimulation of neural circuits {cite}`markowitz2023spontaneous`. While there is still room for methodological improvement, these methods have already transformed the landscape of systems neuroscience."
   ]
  }
 ],
 "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.7.11 (default, Jul 27 2021, 07:03:16) \n[Clang 10.0.0 ]"
  },
  "orig_nbformat": 4,
  "vscode": {
   "interpreter": {
    "hash": "5b28c5bd4ee93d765ebe901023d5522822fb8ad083dac3187c5545022f913719"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}