Added a new tutorial describing the matrix profile

Author: seanlaw

docs/Tutorial_0.ipynb

@@ -0,0 +1,305 @@
+{
+ "cells": [
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "# Tutorial 0 - The Matrix Profile\n",
+ "\n",
+ "At its core, the STUMPY library efficiently computes something called a <i><b>matrix profile</b>, a vector that stores the z-normalized Euclidean distance between any subsequence within a time series and its nearest neigbor</i>.\n",
+ "\n",
+ "To fully understand what this means, let's take a step back and start with a simple illustrative example along with a few basic definitions:"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Time Series with Length n = 13"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 1,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "time_series = [0, 1, 3, 2, 9, 1, 14, 15, 1, 2, 2, 10, 7]\n",
+ "n = len(time_series)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "To analyze this time series with length `n = 13`, we could visualize the data or calculate global summary statistics (i.e., mean, median, mode, min, max). If you had a much longer time series, then you may even feel compelled to build an ARIMA model, perform anomaly detection, or attempt a forecasting model but these methods can be complicated and may often have false positives or no interpretable insights.\n",
+ "\n",
+ "<img src=\"images/time_series_viz.jpeg\" width=800>\n",
+ "\n",
+ "However, if we were to apply <i>Occam's Razor</i>, then what is the most <b>simple and intuitive</b> approach that we could take analyze to this time series?\n",
+ "\n",
+ "To answer this question, let's start with our first defintion:"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Subsequence /ˈsəbsəkwəns/ noun\n",
+ "\n",
+ "### <i>a part or section of the full time series</i>\n",
+ "\n",
+ "So, the following are all considered subsequences of our `time_series` since they can all be found in the time series above.\n",
+ "\n",
+ "<div style='margin-left: auto;margin-right: auto; text-align: center;'><img src=\"images/subsequence_a.jpeg\" width=400><img src=\"images/subsequence_b.jpeg\" width=400><img src=\"images/subsequence_c.jpeg\" width=400></div>"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 2,
+ "metadata": {},
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "[0, 1]\n",
+ "[9, 1, 14]\n",
+ "[3, 2, 9, 1, 14, 15, 1, 2]\n"
+ ]
+ }
+ ],
+ "source": [
+ "print(time_series[0:2])\n",
+ "print(time_series[4:7])\n",
+ "print(time_series[2:10])"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "We can see that each subsequence can have a different sequence length that we'll call `m`. So, for example, if we choose `m = 4`, then we can think about how we might compare any two subsequences of the same length."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 3,
+ "metadata": {},
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "[0, 1, 3, 2] [1, 2, 2, 10]\n"
+ ]
+ }
+ ],
+ "source": [
+ "m = 4\n",
+ "i = 0 # starting index for the first subsequence\n",
+ "j = 8 # starting index for the second subsequence\n",
+ "\n",
+ "subseq_1 = time_series[i:i+m]\n",
+ "subseq_2 = time_series[j:j+m]\n",
+ "\n",
+ "print(subseq_1, subseq_2)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "<div style='margin-left: auto;margin-right: auto; text-align: center;'><img src=\"images/subsequences.jpeg\" width=800></div>"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "One way to compare any two subsequences is to calculate what is called the Euclidean distance.\n",
+ "\n",
+ "## Euclidean Distance /yo͞oˈklidēən/ /ˈdistəns/ noun\n",
+ "\n",
+ "### <i>the straight-line distance between two points</i>\n",
+ "\n",
+ "<div style='margin-left: auto;margin-right: auto; text-align: center;'><img src=\"images/euclidean_distance.jpeg\" width=800></div>"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 4,
+ "metadata": {},
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "The square root of 67 = 8.18535277187245\n"
+ ]
+ }
+ ],
+ "source": [
+ "import math\n",
+ "\n",
+ "D = 0\n",
+ "for k in range(m):\n",
+ " D += (time_series[i+k] - time_series[j+k])**2\n",
+ "print(f\"The square root of {D} = {math.sqrt(D)}\")"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Distance Profile - Pairwise Euclidean Distances\n",
+ "\n",
+ "Now, we can take this a step further where he keep one subsequence the same (reference subsequence), change the second subsequence in a sliding window manner, and compute the Euclidean distance for each window. The resulting vector of pairwise Euclidean distances is also known as a <i><b>distance profile</b></i>.\n",
+ "\n",
+ "<div style='margin-left: auto;margin-right: auto; text-align: center;'><img src=\"images/pairwise_euclidean_distance.gif\" width=800></div>\n",
+ "\n",
+ "Of course, not all of these distances are useful. Specifically, the distance for the self match (or trivial match) isn't informative since the distance will be always be zero when you are comparing a subsequence with itself. So, we'll ignore it and, instead, take note of the next smallest distance from the distance profile and choose that as our best match:\n",
+ "\n",
+ "<div style='margin-left: auto;margin-right: auto; text-align: center;'><img src=\"images/trivial_match.jpeg\" width=800></div>\n",
+ "\n",
+ "Next, we can shift our reference subsequence over one element at a time and repeat the same sliding window process to compute the distance profile for each new reference subsequence.\n",
+ "\n",
+ "<div style='margin-left: auto;margin-right: auto; text-align: center;'><img src=\"images/distance_matrix.gif\" width=800></div>\n",
+ "\n",
+ "## Distance Matrix\n",
+ "\n",
+ "If we take all of the distance profiles that were computed for each reference subsequence and stack them one on top of each other then we get something called a <i><b>distance matrix</b></i>\n",
+ "\n",
+ "<div style='margin-left: auto;margin-right: auto; text-align: center;'><img src=\"images/distance_matrix.jpeg\" width=800></div>\n",
+ "\n",
+ "## The Real Problem - The Brute Force Approach\n",
+ "\n",
+ "Now, it might seem pretty straightforward at this point but what we need to do is consider how to compute this full distance matrix efficiently. Let's start with the brute force approach: "
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 5,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "for i in range(n-m+1):\n",
+ " for j in range(n-m+1):\n",
+ " D = 0\n",
+ " for k in range(m):\n",
+ " D += (time_series[i+k] - time_series[j+k])**2\n",
+ " D = math.sqrt(D)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "At first glance, this may not look too bad but if we start considering both the computational complexity as well as the spatial complexity then we begin to understand the real problem. It turns out that, for longer time series (i.e., <i>n >> 10,000</i>) the computational complexity is <i>O(n<sup>2</sup>m)</i> (as evidenced by the three for loops in the code above) and the spatial complexity for storing the full distance matrix is <i>O(n<sup>2</sup>)</i>. \n",
+ "\n",
+ "To put this into perspective, imagine if you had a single sensor that collected data 20 times/min over the course of 5 years. This would result:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 6,
+ "metadata": {},
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "There would be n = 52416000 data points\n"
+ ]
+ }
+ ],
+ "source": [
+ "n = 20 * 60 * 24 * 364 * 5 # 20 times/min x 60 mins/hour x 24 hours/day x 365 days/year x 5 years\n",
+ "print(f\"There would be n = {n} data points\")"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Assuming that each calculation in the inner loop takes 0.0000001 seconds then this would take:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 7,
+ "metadata": {},
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "It would take 137371850.1792 seconds to compute\n"
+ ]
+ }
+ ],
+ "source": [
+ "time = 0.0000001 * (n * n - n)/2\n",
+ "print(f\"It would take {time} seconds to compute\")"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "<b>Which is equivalent to 1,598.7 days (or 4.4 years) and 11.1 PB of memory to compute!</b> So, it is clearly not feasible to compute the distance matrix using our naive brute force method. And this is where the <i><b>matrix profile</b></i> comes into play. Recall,\n",
+ "\n",
+ "## Matrix Profile /ˈmātriks/ /ˈprōˌfīl/ noun \n",
+ "\n",
+ "### a vector that stores the (z-normalized) Euclidean distance between any subsequence within a time series and its nearest neigbor\n",
+ "\n",
+ "Practically, what this means is that the matrix profile is only interested in storing the smallest non-trivial distances from each distance profile, which significantly reduces the spatial complexity to O(n):\n",
+ "\n",
+ "<div style='margin-left: auto;margin-right: auto; text-align: center;'><img src=\"images/matrix_profile.gif\" width=800></div>\n",
+ "\n",
+ "But we still haven't figured out how to reduce the computational complexity and this is precisely the problem that is being addressed by [STUMPY](https://github.com/TDAmeritrade/stumpy)."
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## STUMPY\n",
+ "\n",
+ "In the fall of 2016, researchers from the [University of California, Riverside](https://www.cs.ucr.edu/~eamonn) and the [University of New Mexico](https://www.cs.unm.edu/~mueen/) published a beautiful set of [back-to-back papers](https://www.cs.ucr.edu/~eamonn/MatrixProfile.html) that described an <u>exact method</u> called <i><b>STOMP</b></i> for computing the matrix profile for any time series with a computational complexity of O(n<sup>2</sup>)! They also further demonstrated this using GPUs and this method approach is called <i><b>GPU-STOMP</b></i>.\n",
+ "\n",
+ "With the academics, data scientists, and developers in mind, we have taken some of these concepts and have open sourced STUMPY, a powerful and scalable library that efficiently computes the matrix profile. And, thanks to other open source software such as [Numba](http://numba.pydata.org/) and [Dask](https://dask.org/), our implementation can be both highly parallelized (for a single server) and highly distributed (across multiple servers). We've tested STUMPY on as many as 256 CPU cores that are spread across 32 servers and have achieved similiar [performance](https://github.com/TDAmeritrade/stumpy#performance) to the published GPU-STOMP work."
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Conclusion\n",
+ "\n",
+ "We are excited to share STUMPY with you and please reach out and tell us know how it has enabled your time series analysis work. We'd love to hear from you!"
+ ]
+ }
+ ],
+ "metadata": {
+ "kernelspec": {
+ "display_name": "Python 3",
+ "language": "python",
+ "name": "python3"
+ },
+ "language_info": {
+ "codemirror_mode": {
+ "name": "ipython",
+ "version": 3
+ },
+ "file_extension": ".py",
+ "mimetype": "text/x-python",
+ "name": "python",
+ "nbconvert_exporter": "python",
+ "pygments_lexer": "ipython3",
+ "version": "3.6.8"
+ }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 2
+}

docs/Tutorial_2.ipynb

@@ -4,7 +4,7 @@
  "cell_type": "markdown",
  "metadata": {},
  "source": [
- "# Analyzing Taxi Time Series Data with Time Series Chains\n",
+ "# Tutorial 2 - Analyzing Taxi Time Series Data with Time Series Chains\n",
  "\n",
  "This example utilizes the main takeways from the research papers: [Matrix Profile VII](https://www.cs.ucr.edu/~eamonn/chains_ICDM.pdf).\n",
  "\n",
@@ -260,7 +260,7 @@
  "cell_type": "markdown",
  "metadata": {},
  "source": [
- "# Resources\n",
+ "## Resources\n",
  "\n",
  "[Matrix Profile VII](https://www.cs.ucr.edu/~eamonn/chains_ICDM.pdf)"
  ]