scientific_details_of_algorithms/streaming_median/streamingMedian.py.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Streaming Algorithms in Machine Learning\n",
    "\n",
    "In this notebook, we will use an extremely simple \"machine learning\" task to learn about streaming algorithms. We will try to find the median of some numbers in batch mode, random order streams, and arbitrary order streams.\n",
    "The idea is to observe first hand the advantages of the streaming model as well as to appreciate some of the complexities involved in using it.\n",
    "\n",
    "The task at hand will be to approximate the median (model) of a long sequence of numbers (the data). This might seem to have little to do with machine learning. We are used to thinking of a median, $m$, of number $x_1,\\ldots,x_n$ in the context of statistics as the number, $m$, which is smaller than at most half the values $x_i$ and larger than at most half the values $x_i$.\n",
    "\n",
    "Finding the median, however, also solves a proper machine learning optimization problem (albeit a simple one). The median minimizes the following clustering-like objective function\n",
    "$$m = \\min_x \\frac1n\\sum_i|x - x_i|.$$\n",
    "In fact, the median is the solution to the well studied k-median clustering problem in one dimension and $k=1$. Moreover, the extension to finding all quantiles is common in feature transformations and an important ingedient in speeding up decission tree training."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Batch Algorithms\n",
    "Let's first import a few libraries and create some random data. \n",
    "Our data will simple by $100,000$ equally spaced points between $0$ and $1$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "n = 100000\n",
    "data = np.linspace(0, 1, n)\n",
    "np.random.shuffle(data)\n",
    "\n",
    "\n",
    "def f(x, data):\n",
    "    return sum(abs(x - datum) for datum in data) / len(data)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's look at the data to make sure everything is correct."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "# Plotting every 10th point\n",
    "plt.scatter(range(0, n, 100), data[0:n:100], vmin=0, vmax=1.0)\n",
    "plt.ylim((0.0, 1.0))\n",
    "plt.xlim((0, n))\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Computing the median brute force is trivial. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from math import floor\n",
    "\n",
    "\n",
    "def batchMedian(data):\n",
    "    n = len(data)\n",
    "    median = sorted(data)[int(floor(n / 2))]\n",
    "    return median\n",
    "\n",
    "\n",
    "median = batchMedian(data)\n",
    "print(\"The median found if {}\".format(median))\n",
    "print(\"The objective value is {}\".format(f(median, data)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The result is, of course, correct ($0.5$).\n",
    "To get the median we sorted the data in $O(n\\log n)$ time even though QuickSelect would have been faster ($O(n)$). The algorithm speed is not the main issue here though. The main drawback of this algorithm is that it must store the entire dataset in memory. For either sorting or quickSelect the algorithm must also duplicate the array. Binary search is also a possible solution which doesn't require data duplication but does require $O(\\log(n))$ passes over the data.\n",
    "When the data is large this is either very expensive or simply impossible."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Streaming Algorithms (Random Order, SGD)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the streaming model, we assume only an iterator over the data is given. That is, we can only make one pass over the data. Moreover, the algorithm is limited in its memory footprint and the limit is much lower than the data size. Otherwise, we could \"cheat\" by storing all the data in memory and executing the batch mode algorithm.\n",
    "\n",
    "Gradient Descent (GD) type solutions are extremely common in this setting and are, de facto, the only mechanism for optimizing neural networks. In gradient descent, a step is taken in the direction opposite of the gradient. In one dimension, this simply means going left if the derivative is positive or right if the derivative is negative. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "xs = list(np.linspace(-1.0, 2.0, 50))\n",
    "ys = [f(x, data) for x in xs]\n",
    "plt.plot(xs, ys)\n",
    "plt.ylim((0.0, 2.0))\n",
    "plt.xlim((-1.0, 2.0))\n",
    "ax = plt.axes()\n",
    "ax.arrow(-0.5, 1.1, 0.3, -0.3, head_width=0.05, head_length=0.1, fc=\"k\", ec=\"k\")\n",
    "ax.arrow(1.5, 1.1, -0.3, -0.3, head_width=0.05, head_length=0.1, fc=\"k\", ec=\"k\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In **Stochastic Gradience Descent**, one only has a stochastic (random) unbiased estimator of the gradient. So, instead of computing the gradient of $\\frac1n\\sum_i|x - x_i|$ we can compute the gradient of $|x - x_i|$ where $x_i$ is chosen **uniformly at random** from the data. Note that a) the derivative of $|x - x_i|$ is simply $1$ if $x > x_i$ and $-1$ otherwise and b) the *expectation* of the derivative is exactly equal to the derivative of the overall objective function. \n",
    "\n",
    "Comment: the authors of the paper below suggest essentially this algorithm but do not mention the connection to SGD for some reason.\n",
    "\n",
    "Frugal Streaming for Estimating Quantiles: One (or two) memory suffices: Qiang Ma, S. Muthukrishnan, Mark Sandler"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from math import sqrt\n",
    "\n",
    "\n",
    "def sgdMedian(data, learningRate=0.1, initMedianValue=0):\n",
    "    median = initMedianValue\n",
    "    for (t, x) in enumerate(data):\n",
    "        gradient = 1.0 if x < median else -1.0\n",
    "        median = median - learningRate * gradient / sqrt(t + 1)\n",
    "    return median\n",
    "\n",
    "\n",
    "median = sgdMedian(data, learningRate=0.1, initMedianValue=0)\n",
    "\n",
    "print(\"The median found if {}\".format(median))\n",
    "print(\"The objective value is {}\".format(f(median, data)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The result isn't exactly $0.5$ but it is pretty close. If this was a real machine learning problem, matching the objective up to the 5th digit of the true global minimum would have been very good.\n",
    "\n",
    "Why does this work? Let's plot our objective function to investigate further."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It should not come as a big surprise to you that the objective function is convex. After all, it is the sum of convex functions (absolute values). It is a piece-wise linear curve that approximates a parabole in the range $(0,1)$ and is linear outside that range. Therefore, gradient descent is guaranteed to converge."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "SGD significantly more efficient than sorting or even QuickSelect. More importantly, its memory footprint is tiny, a handful of doubles, *regardless of the size of the data*!!!\n",
    "\n",
    "This is a huge advantage when operating with large datasets or with limited hardware.\n",
    "Alas, SGD has some subtleties that make it a little tricky to use sometimes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# SGD needs to be initialized carfully\n",
    "median = sgdMedian(data, learningRate=0.1, initMedianValue=100.0)\n",
    "\n",
    "print(\"The median found if {}\".format(median))\n",
    "print(\"The objective value is {}\".format(f(median, data)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# SGD needs to set step sizes corectly (controled via the learing rate)\n",
    "median = sgdMedian(data, learningRate=0.001, initMedianValue=0.0)\n",
    "\n",
    "print(\"The median found if {}\".format(median))\n",
    "print(\"The objective value is {}\".format(f(median, data)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "These issues are usually alleviated by adaptive versions of SGD. Enhancements to SGD such as second order (based) methods, adaptive learning rate, and momentum methods may help in these situations but still require tuning in many cases. A common approach is to use many epochs."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "median = 0.0\n",
    "numEpochs = 100\n",
    "for i in range(numEpochs):\n",
    "    median = sgdMedian(data, learningRate=0.001, initMedianValue=median)\n",
    "\n",
    "print(\"The median found if {}\".format(median))\n",
    "print(\"The objective value is {}\".format(f(median, data)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "While clearly much less efficient than a single pass, increasing the number of epochs seemed to have solved the problem. Machine learning practitioners can relate to this result. That is, SGD is a great algorithm IF one finds good parameters for initialization, learning rate, number of epochs etc.\n",
    "\n",
    "One of the main challenges in designing fundamentally better SGD-based streaming algorithms is in adaptively controlling these parameters during the run of the algorithm.\n",
    "\n",
    "It is important to mention that there are also fundamentally better algorithms than SGD for this problem. See for example:\n",
    "\n",
    "_Sudipto Guha, Andrew McGregor_ <br>\n",
    "Stream Order and Order Statistics: Quantile Estimation in Random-Order Streams. <br>\n",
    "_SIAM J. Comput. 38(5): 2044-2059 (2009)_\n",
    "\n",
    "Unfortunately, we don't have time to dive into that..."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Trending data poses a challenge...\n",
    "SGD and the other above algorithm have a fundamental drawback. They inherently rely on the fact that the data is random. For SGD, the gradient of the loss on a single point (or minibatch) must be an estimator of the global gradient. This is not true if trends in data make its statistics change (even slightly) over time. Let's simulate this with our data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "\n",
    "# SGD also  depends on the data being reandomly suffeled\n",
    "n, k = len(data), 10\n",
    "minibatches = [data[i : i + k] for i in range(0, n, k)]\n",
    "minibatches.sort(key=sum)\n",
    "trendyData = np.array(minibatches).reshape(n)\n",
    "\n",
    "# Plotting every 10th point in the trending dataset\n",
    "plt.scatter(range(0, n, 100), trendyData[0:n:100], vmin=0, vmax=1.0)\n",
    "plt.ylim((0.0, 1.0))\n",
    "plt.xlim((0, n))\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "median = sgdMedian(trendyData, learningRate=0.1, initMedianValue=0.0)\n",
    "\n",
    "print(\"The median found if {}\".format(median))\n",
    "print(\"The objective value is {}\".format(f(median, data)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Streaming Algorithms (single pass, arbitrary order)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "One way not to be fooled by trends in data and to sample from it. \n",
    "The algorithm uses Reservoir Sampling to obtain $k$ (in this case $k=1000$) uniformly chosen samples from the stream. Then, compute the batch median of the sample. \n",
    "\n",
    "The main drawback of sampling is that we now use more memory. Roughly the sample size $k$ ($k=1000$ here). This much more than $O(1)$ needed for SGD. Yet, it has some very appealing properties. Sampling very efficient ($O(1)$ per update), it is very simple to implement, it doesn't have any numeric sensitivities or tunable input parameters, and it is provably correct. \n",
    "\n",
    "_(For the sake of simplicity below we use python's builtin sample function rather than recode reservoir sampling)_"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from random import sample\n",
    "\n",
    "\n",
    "def sampleMedian(data):\n",
    "    k = 300\n",
    "    samples = sample(list(data), k)\n",
    "    return batchMedian(samples)\n",
    "\n",
    "\n",
    "median = sampleMedian(trendyData)\n",
    "print(\"The median found if {}\".format(median))\n",
    "print(\"The objective value is {}\".format(f(median, data)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As you can see, sampling provides relatively good results.\n",
    "\n",
    "Nevertheless, there is something deeply dissatisfying about it. The algorithm was given $100,000$ points and used on $1,000$ of them. I other words, it would have been just as accurate had we collected only $1\\%$ of the data.\n",
    "\n",
    "Can we do better? Can an algorithm simultaneously take advantage of all the data, have a fixed memory footprint, and not be sensitive to the order in which the data is consumed? The answer is _yes!_. These are known in the academic literature as Sketching (or simply streaming) algorithms.\n",
    "\n",
    "Specifically for approximating the median (or any other quantile), there is a very recent result that shows how best to achieve that:\n",
    "\n",
    "_Zohar S. Karnin, Kevin J. Lang, Edo Liberty_ <br>\n",
    "Optimal Quantile Approximation in Streams. <br>\n",
    "FOCS 2016: 71-78\n",
    "\n",
    "The following code is a hacky version of the algorithm described in the paper above. Warning: this function will not work for streams much longer than $100,000$!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from kll300 import KLL300\n",
    "from bisect import bisect\n",
    "\n",
    "\n",
    "def sketchMedian(data):\n",
    "    sketch = KLL300()\n",
    "    for x in data:\n",
    "        sketch.update(x)\n",
    "        assert sketch.size <= 300  # making sure there is no cheating involved...\n",
    "    items, cdf = sketch.cdf()\n",
    "    i = bisect(cdf, 0.5)\n",
    "    median = items[i]\n",
    "    return median\n",
    "\n",
    "\n",
    "median = sketchMedian(trendyData)\n",
    "print(\"The median found if {}\".format(median))\n",
    "print(\"The objective value is {}\".format(f(median, data)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that sketchMedian and sampleMedian both retain at most 300 items from the stream.\n",
    "Still, the sketching solution is significantly more accurate.\n",
    "Note that both sampling and sketching are randomized algorithms. \n",
    "It could be that sampling happens to be more accurate than sketching for any single run. But, as a whole, you should expect the sketching algorithm to be much more accurate. \n",
    "\n",
    "If you are curious about what sketchMedian actually does, you should look here:\n",
    "* Academic paper - https://arxiv.org/abs/1603.05346\n",
    "* JAVA code as part of datasketches - https://github.com/DataSketches/sketches-core/tree/master/src/main/java/com/yahoo/sketches/kll\n",
    "* Scala code by Zohar Karnin - https://github.com/zkarnin/quantiles-scala-kll\n",
    "* Python experiments by Nikita Ivkin - https://github.com/nikitaivkin/quantilesExperiments\n",
    "\n",
    "The point is, getting accurate and stable streaming algorithms is complex. This is true even for very simple problems (like the one above). But, if one can do that, the benefits are well worth it."
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "conda_python3",
   "language": "python",
   "name": "conda_python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.2"
  },
  "notice": "Copyright 2017 Amazon.com, Inc. or its affiliates. All Rights Reserved.  Licensed under the Apache License, Version 2.0 (the \"License\"). You may not use this file except in compliance with the License. A copy of the License is located at http://aws.amazon.com/apache2.0/ or in the \"license\" file accompanying this file. This file is distributed on an \"AS IS\" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. See the License for the specific language governing permissions and limitations under the License."
 },
 "nbformat": 4,
 "nbformat_minor": 2
}