diff --git a/docs/overlap_add_convolution.gif b/docs/overlap_add_convolution.gif new file mode 100644 index 0000000..995e023 Binary files /dev/null and b/docs/overlap_add_convolution.gif differ diff --git a/docs/overlap_add_math.ipynb b/docs/overlap_add_math.ipynb new file mode 100644 index 0000000..df595ed --- /dev/null +++ b/docs/overlap_add_math.ipynb @@ -0,0 +1,1200 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "8330f3d0-1e58-4b8f-9619-f0baa21a4f16", + "metadata": {}, + "source": [ + "# **Compute Sliding Dot Product using Convolution**" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "id": "9aa1b0cf-4684-462d-8e1a-b1e790b582d1", + "metadata": {}, + "outputs": [], + "source": [ + "import functools\n", + "import numpy as np\n", + "\n", + "from IPython.display import HTML\n", + "import base64" + ] + }, + { + "cell_type": "markdown", + "id": "5c9355b9-0bd4-4740-9721-7394a19f4588", + "metadata": {}, + "source": [ + "# 0. Objective" + ] + }, + { + "cell_type": "markdown", + "id": "199eb655-fde9-4ee1-82e4-55bad5892f22", + "metadata": {}, + "source": [ + "The paper [\"Matrix Profile I\"](https://www.cs.ucr.edu/~eamonn/PID4481997_extend_Matrix%20Profile_I.pdf) uses the [MASS algorithm](https://www.cs.unm.edu/~mueen/FastestSimilaritySearch.html) to compute the distance between a query $Q$ and every subsequence of length $len(Q)$ in $T$. As part of this algorithm, the sliding dot product (sdp) is calculated. The [MASS paper](https://link.springer.com/epdf/10.1007/s10618-024-01005-2?sharing_token=067pAxnaLDvz89q1n_GJt_e4RwlQNchNByi7wbcMAY6HNwOWuMxQSNE3HcKcuL8siHB8L8krJpchQVaGvicUgoegxxV7BWgaU4Y9evg1FU1LGtxlvM9A5UrxrtkqDHyvoOPk_ttVNps-6-LSRf1vuRgFWJI7qrkqWGitjsXfRsw%3D) proposes different variants of MASS for computing the sliding dot product. All of the proposed methods use some form of convolution at their core. This notebook is created to help us understand how convolution comes into the picture for computing the sliding dot product. In particular, this notebook answers the following questions:\n", + "\n", + "1. **How is the sliding dot product (sdp) related to linear convolution?** This can help us discover the fundamental relationship between sdp and linear convolution\n", + "2. **How can we compute the convolution faster?** Convolution is a well-studied area. The goal is to review the existing methods to eventually compute sdp faster\n", + "3. **How can we compute convolution faster when $T$ is very long?** This is where we learn about the overlap-add method, a divide-and-conquer approach for convolution" + ] + }, + { + "cell_type": "markdown", + "id": "03f9d2f8-1560-4a5d-93a2-ec82fb422f6f", + "metadata": {}, + "source": [ + "The content of this notebook is based on the book \"Discrete-Time Signal Processing\" by Alan V. Oppenheim and Ronald W. Schafer. The book is [publicly available in MIT OpenCourseWare](https://ocw.mit.edu/courses/res-6-dtsp-discrete-time-signal-processing/resources/mitres_6-dtsp_s26_thirdedition_pdf/). Most textbooks use $x$ and $h$ notation as the two signals in convolution. However, in this notebook, we would like to use $Q$ and $T$ as those are the names of arrays that are used in STUMPY. Also, we would like to introduce $Q'$, which is simply the reverse of $Q$. So, if $Q$ has length $m$ (indices $0,\\ldots,m-1$), then its reverse $Q'$ is an array with the same length, and it is defined as follows:\n", + "\n", + "\\begin{equation}\n", + "Q'[j] = Q[m-1-j], \\quad 0 \\le j \\le m-1\n", + "\\end{equation}" + ] + }, + { + "cell_type": "markdown", + "id": "3dda604a-3a3e-42f4-90c7-7b01ada16337", + "metadata": {}, + "source": [ + "# 1. How is the sliding dot product(sdp) related to linear convolution?\n", + "\n", + "We first start with the definition of linear convolution for finite sequences. Then, we will discover its relationship with sdp." + ] + }, + { + "cell_type": "markdown", + "id": "877af78c-0026-44c0-8295-0f75db075d97", + "metadata": {}, + "source": [ + "## 1.1 The Linear Convolution of Two Finite-Length Sequences $Q'$ and $T$" + ] + }, + { + "cell_type": "markdown", + "id": "9ee474bd-1229-4749-80fe-f90c4a135931", + "metadata": {}, + "source": [ + "Consider two finite-length sequences, $Q'$ and $T$, of lengths $m$ and $n$, respectively, where $m < n$. Their linear convolution, $C$, can be computed as follows:" + ] + }, + { + "cell_type": "markdown", + "id": "ad4a63b5-2d91-4875-b4c2-bc7ce11b6af9", + "metadata": {}, + "source": [ + "\\begin{equation}\n", + "C[idx] = \\sum_{i=-\\infty}^{\\infty}{T[i] \\times Q'[idx-i]}\n", + "\\end{equation}" + ] + }, + { + "cell_type": "markdown", + "id": "7b2a4283-f6de-4ef8-9bac-2f23e62aca02", + "metadata": {}, + "source": [ + "where $T[.]$ and $Q'[.]$ are zeros at any index that is outside of their range. The linear convolution can be seen as a flip and slide operation. So, $Q'$ is flipped (reversed) and slid across $T$. The output at index $idx$ is the sum of element-wise product of overlapping elements. So, as long as there is at least one overlapping element, there can be an output." + ] + }, + { + "cell_type": "markdown", + "id": "c7134e06-57df-41a7-aa82-c0e42ecb5ade", + "metadata": {}, + "source": [ + "
\n", + " Note: Linear convolution is NOT the same as the sliding dot product (sdp). In sdp, an overlap needs to cover the full query. However, in linear convolution, an overlap of one element can still contribute to the output.\n", + "
" + ] + }, + { + "cell_type": "markdown", + "id": "8028ee30-3244-4c79-86ca-22e2b8876694", + "metadata": {}, + "source": [ + "As stated earlier, $T[i]$ is zero if $i$ is outside of the range $0 \\le i \\le n-1$. Hence, the boundaries of the sum operation in the previous equation can be changed to $[0, n-1]$, i.e.\n", + "\n", + "\\begin{equation}\n", + "C[idx] = \\sum_{i=0}^{n-1}{T[i] \\times Q'[idx-i]}\n", + "\\end{equation}" + ] + }, + { + "cell_type": "markdown", + "id": "593a8276-49cb-43d2-9300-8d29e67d3104", + "metadata": {}, + "source": [ + "Note that $Q'[idx-i]$ can be a non-zero value only when $0 \\le idx-i \\le m-1$, or equivalently $i \\le idx \\le i+m-1$. So, as the index $i$ changes from $0$ to $n-1$, the interval $i \\le idx \\le i+m-1$ changes accordingly, and the range of $idx$ becomes the union of these intervals, i.e.\n", + "\n", + "\\begin{equation}\n", + "\\bigcup_{i=0}^{n-1} [\\,i,\\ i+m-1\\,]\n", + "\\end{equation}\n", + "\n", + "So,\n", + "* min(idx) is coming from the lower bound of interval $[\\,i,\\ i+m-1\\,]$ when $i=0$. This gives: $min(idx)=0$\n", + "* max (idx) is coming from the upper bound of interval $[\\,i,\\ i+m-1\\,]$ when $i=n-1$. This gives: $max(idx) = (n-1)+m-1 = n+m-2$\n", + "\n", + "When $idx$ changes from $0$ to $n+m-2$, there **exists at least one $i$** such that both $T[i]$ and $Q'[idx-i]$ can be non-zero values. This mathematically shows that the linear convolution $C$ can have values at indices $\\set{0, 1, ..., n + m - 2}$. The length of the linear convolution is $n + m - 1$." + ] + }, + { + "cell_type": "markdown", + "id": "ecb14f7d-cde6-42db-aece-2fb40536d68c", + "metadata": {}, + "source": [ + "## 1.2 Sliding dot product (of $Q$ and $T$) as a slice of linear convolution (of $Q'$ and $T$)" + ] + }, + { + "cell_type": "markdown", + "id": "e638d3d7-cf2a-4355-aff3-a9192400defb", + "metadata": {}, + "source": [ + "Let's check out the values of linear convolution at different indices. For convenience, the linear convolution equation is repeated below:\n", + "\n", + "\\begin{equation}\n", + "C[idx] = \\sum_{i=0}^{n-1}{T[i] \\times Q'[idx-i]}\n", + "\\end{equation}\n", + "\n", + "Recall that the linear convolution $C$ can have values at indices $\\set{0, 1, ..., n + m - 2}$.\n", + "\n", + "> **NOTE** Let $T_{i}$ denote the sequence $\\{T[i], T[i+1], ..., T[i+m-1]\\}$\n", + "\n", + "* $idx=0 \\implies C[0]=T[0]Q'[0]$\n", + "* $idx=1 \\implies C[1]=T[0]Q'[1] + T[1]Q'[0]$\n", + "* ...\n", + "* $\\textcolor{red}{idx=m-1 \\implies C[m-1]=T[0]Q'[m-1] + T[1]Q'[m-2] + ... + T[m-1]Q'[0] = T[0]Q[0] + ... + T[m-1]Q[m-1] = T_{0}.Q}$\n", + "* $\\textcolor{red}{idx=m \\implies C[m]=T[1]Q'[m-1] + T[2]Q'[m-2] + ... + T[m]Q'[0] = T[1]Q[0] + ... + T[m]Q[m-1] = T_{1}.Q}$\n", + "* ...\n", + "* $\\textcolor{red}{idx=n-1 \\implies C[n-1]=T[n-m]Q'[(n-1)-(n-m)] + T[n-m+1]Q'[(n-1)-(n-m+1)] + ... + T[n-1]Q'[0] = ... = T_{n-m}.Q}$\n", + "* ...\n", + "* $idx=n+m-2 \\implies C[n+m-2]=T[n-1]Q'[(n+m-2)-(n-1)] = T[n-1]Q'[m-1]$\n", + "\n", + "When $m-1 \\le idx \\le n-1$ (the cases shown in red), $C[idx]$ actually becomes the dot product between $T_{idx-(m-1)}$ and $Q$ (the reverse of $Q'$). Therefore, the linear convolution of $Q'$ and $T$ contains the sliding dot product between $Q$ and $T$ in the slice in range $[m-1,n)$. \n", + "
\n", + "
\n", + "\n", + "#### **Let's summarize what we have learned so far:**\n", + "* The linear convolution $C$, between $Q'$ and $T$, has values at indices $\\set{0, 1, ..., n + m - 2}$. The length of convolution is $n + m - 1$.\n", + "* The sliding dot product between $Q$ and $T$ is a slice of $C$, the linear convolution between $Q'$ (the reverse of $Q$) and $T$. The slice is for the range $[m-1, n)$.\n", + "\n", + "\n", + "So, once we calculate the linear convolution between $Q'$ and $T$, we have the sliding dot product between $Q$ and $T$ for free!! Note that we haven't talked about the efficiency of computation. This is discussed in the next section.\n", + "\n", + "
\n", + " Note: To compute the sdp between Q and T, we can compute the linear convolution between Q' (the reverse of Q) and T, and select the slice(m-1,n)\n", + "
" + ] + }, + { + "cell_type": "markdown", + "id": "8475be6a-9515-43ae-a8d7-083b8db3af51", + "metadata": {}, + "source": [ + "# 2. How can we compute the convolution faster?\n", + "\n", + "In this section, we make three attempts, and each attempt is to help us eventually improve the efficiency of computing sdp." + ] + }, + { + "cell_type": "markdown", + "id": "44c2f68d-3424-4417-bf4d-474b31373688", + "metadata": {}, + "source": [ + "## 2.1 Attempt I: Circular Convolution via Fast Fourier Transform (FFT) \n", + "\n", + "To understand this part, we first need to understand what \"Circular Convolution\" is, and how it is related to the linear convolution. Then, we will see how FFT can speed it up." + ] + }, + { + "cell_type": "markdown", + "id": "52bcf566-3ac0-483d-81a3-57fa3b5766c3", + "metadata": {}, + "source": [ + "### 2.1.1 First Step: Convert linear convolution to circular convolution" + ] + }, + { + "cell_type": "markdown", + "id": "2c4176e2-60f3-410f-ad2d-6c1d29cca6cf", + "metadata": {}, + "source": [ + "**The periodic convolution** between two sequences that are both periodic and their periods are the same, say $N$, is a periodic sequence with same period $N$. We can recover this periodic sequence by computing one period of it. The following equation is called **\"Circular Convolution\"** and it shows how the values of one period, for range $[0, N)$, are computed.\n", + "\n", + "\\begin{equation}\n", + "\\tilde{C}_{N}[idx] = \\sum_{i=0}^{N-1} \\tilde{T}_{N}[i] \\times \\tilde{Q'}_{N}[(idx-i)_{N}], \\quad 0 \\le idx \\le N-1\n", + "\\end{equation}\n", + "\n", + "where, the subscript $N$ represents one period. So, $\\tilde{T}_{N}$ and $\\tilde{Q'}_{N}$ represent one period of the periodic sequences $\\tilde{T}$ and $\\tilde{Q'}$, respectively. Similarly $\\tilde{C}_{N}$ represents one period of the periodic output. The index $(idx-i)_{N}$ means \"$(idx-i)$ modulo $N$\", and it is always between $0$ and $N-1$." + ] + }, + { + "cell_type": "markdown", + "id": "e10c5304-a1a1-4cbd-9b26-266e229bdc57", + "metadata": {}, + "source": [ + "We now show that the linear convolution between $Q'$ and $T$ is equivalent to the circular convolution $C_{N}$ between $\\tilde{Q'}$ and $\\tilde{T}$ if the three following conditions are satisfied.\n", + "\n", + "* **Condition 1:** $N = n + m - 1$\n", + "* **Condition 2:** $\\tilde{T}_{N}$ is $T$, padded with $N-n$ zeros\n", + "* **Condition 3:** $\\tilde{Q'}_{N}$ is $Q'$, padded with $N-m$ zeros" + ] + }, + { + "cell_type": "markdown", + "id": "a7e66371-4cd4-4673-957c-b7c548dcbfde", + "metadata": {}, + "source": [ + "#### **Proof:** \n", + "\n", + "Let's write down the equations for the linear convolution and the circular convolution below:\n", + "\n", + "\\begin{equation}\n", + "\\text{Linear Convolution:} \\quad\\quad C[idx] = \\sum_{i=0}^{n-1}{T[i] \\times Q'[idx-i]}, \\quad 0 \\le idx < n+m-1\n", + "\\end{equation}\n", + "\n", + "\\begin{equation}\n", + "\\text{Circular Convolution:} \\quad\\quad C_{N}[idx] = \\sum_{i=0}^{N-1} \\tilde{T}_{N}[i] \\times \\tilde{Q'}_{N}[(idx-i)_{N}], \\quad 0 \\le idx \\le N-1\n", + "\\end{equation}\n", + "\n", + "We want to show that, at any index $idx$, the two values on the left-hand side of equations above are equal when: \n", + "* **Condition 1:** $N = n + m - 1$\n", + "* **Condition 2:** $\\tilde{T}_{N}$ is $T$, padded with $N-n$ zeros\n", + "* **Condition 3:** $\\tilde{Q'}_{N}$ is $Q'$, padded with $N-m$ zeros\n", + " \n", + "**According to the \"condition 1\": $N=n+m-1$**
\n", + "This means $N \\ge n$ since $m$ is at least one. Therefore, we can the rewrite the circular convolution equation by separating the cases where the iterator $i$ is smaller than $n$ from the cases where it is not.\n", + "\n", + "\\begin{equation}\n", + "C_{N}[idx] = \\sum_{i=0}^{n-1} \\tilde{T}_{N}[i] \\times \\tilde{Q'}_{N}[(idx-i)_{N}] + \\sum_{i=n}^{N-1} \\tilde{T}_{N}[i] \\times \\tilde{Q'}_{N}[(idx-i)_{N}]\n", + "\\end{equation}\n", + "\n", + "\n", + "**According to the \"condition 2\": $\\tilde{T}_{N}$ is $T$, padded with $N-n$ zeros.** Therefore, $\\tilde{T}_{N}[i]$ is zero if $n \\le i \\le N-1$. Hence, the second sum is always 0. In addition, in the first summation, $\\tilde{T}_{N}$ can be replaced with $T$ as $i$ falls inside the range $[0, n-1]$. This results in simplifying the circular convolution equation to:\n", + "\n", + "\\begin{equation}\n", + "C_{N}[idx] = \\sum_{i=0}^{n-1} T[i] \\times \\tilde{Q'}_{N}[(idx-i)_{N}]\n", + "\\end{equation}\n", + "\n", + "\n", + "So, instead of the original equations, we need to show the following ones are equivalent:\n", + "\n", + "\\begin{equation}\n", + "\\text{Linear Convolution:} \\quad\\quad C[idx] = \\sum_{i=0}^{n-1}{T[i] \\times Q'[idx-i]}, \\quad 0 \\le idx < n+m-1\n", + "\\end{equation}\n", + "\n", + "\\begin{equation}\n", + "\\text{Circular Convolution:} \\quad\\quad C_{N}[idx] = \\sum_{i=0}^{n-1} T[i] \\times \\tilde{Q'}_{N}[(idx-i)_{N}], \\quad 0 \\le idx < N \\quad \\text{, where N=n+m-1}\n", + "\\end{equation}\n", + "\n", + "As we can see, these two equations are identical except for the term $Q'$ and $\\tilde{Q'}$. So, we only need to show that $\\tilde{Q'}_{N}[(idx-i)_{N}]$ and $Q'[idx-i]$ give the same value. Let $j$ denote $idx-i$. So, all we need to do is to prove that $Q'[j]$ and $\\tilde{Q'}_{N}[(j)_{N}]$ give the same value for different values of $j$.\n", + "\n", + "**Case I: $j \\ge 0$**
\n", + "\n", + "Let's compute the upper bound for $j$. Recall that: \n", + "* $j=idx-i$\n", + "* $0 \\le idx \\le N-1$ (according to range of indices of one period in circular convolution)\n", + "* $0 \\le i \\le n-1$ (according to range of indices of elements in $T$)\n", + "\n", + "Therefore,\n", + "\n", + "\\begin{equation}\n", + "max(j) = max(idx-i) = max(idx) + max(-i) = max(idx) - min(i) = (N-1) - 0 = N-1\n", + "\\end{equation}\n", + "\n", + "So, $j$ changes from 0 to $N-1$. Therefore, $j$ module $N$ becomes $j$, meaning: $\\tilde{Q'}_{N}[(j)_{N}]=Q'_{N}[j]$. \n", + "\n", + "* If $0 \\le j < m$: $Q'_{N}[j] = Q'[j]$ since, according to **the condition 3**, $Q'_{N}$ is simply $Q'$ but padded with zeros.\n", + "* If $m \\le j < N$: $Q'[j]$ is considered to be zero in the linear convolution formula since the index $j$ is not in the range of indices of $Q'$. According to **the condition 3**, $Q'_{N}[j]$ is also zero for $m \\le j < N$ as it falls into the zero-padding zone. So, in this scenario, both $Q'_{N}[j]$ and $Q'[j]$ are zero.\n", + "\n", + "Proof is now complete for this case.\n", + "\n", + "**Case II: $j < 0$**
\n", + "\n", + "In this case, $Q'[j]$ is zero in the linear convolution as it is outside of the range of indices of $Q'$. We now need to prove that $\\tilde{Q'}_{N}[(j)_{N}]$ is zero as well. \n", + "\n", + "Let's compute the lower bound for $j$. Recall that: \n", + "* $j=idx-i$\n", + "* $0 \\le idx \\le N-1$\n", + "* $0 \\le i \\le n-1$\n", + "\n", + "Therefore,\n", + "\n", + "\\begin{equation}\n", + "min(j) = min(idx-i) = min(idx) + min(-i) = min(idx) - max(i) = 0 - (n-1)\n", + "\\end{equation}\n", + "\n", + "So: $-(n-1) \\le j \\le -1$. \n", + "\n", + "Recall that $\\tilde{Q'}$ is a periodic sequence. Therefore, $\\tilde{Q'}_{N}[(j)_{N}] == \\tilde{Q'}_{N}[(j+N)_{N}]$, we can compute the latter instead of the former. The boundaries for $j$ is shown above. Accordingly, the boundaries of $j+N$ are as follows:\n", + "\n", + "\\begin{equation}\n", + "-(n-1)+N \\le j+N \\le -1+N\n", + "\\end{equation}\n", + "\n", + "Recall that $N=n+m-1$. Let's plug that value into $-(n-1)+N$. The inequality becomes:\n", + "\n", + "\\begin{equation}\n", + "(n+m-1)-(n-1) \\le j+N \\le N-1 \\implies m \\le j+N \\le N-1\n", + "\\end{equation}\n", + "\n", + "Since $j+N$ is between $m$ and $N-1$, its value modulo $N$ becomes $j+N$ again, meaning $\\tilde{Q'}_{N}[(j+N)_{N}] = \\tilde{Q'}_{N}[j+N]$. Since the index $j+N$ is $\\ge m$, it falls into the zero-padding zone, and it is zero. \n", + "\n", + "Proof is now complete for this case.\n", + "\n", + "---" + ] + }, + { + "cell_type": "markdown", + "id": "66a707c1-8004-4fe6-b954-ffe039eda441", + "metadata": {}, + "source": [ + "We just showed that the linear convolution can be computed in the form of circular convolution when the arrays are zero-padded properly. We now need to understand how this allows us to compute the linear convolution (hence the sdp) faster!" + ] + }, + { + "cell_type": "markdown", + "id": "8648ff47-4d67-4bec-a91e-4ec6d799cc0e", + "metadata": {}, + "source": [ + "### 2.1.2 Second Step: Compute Circular Convolution using Fast Fourier Transform (FFT)" + ] + }, + { + "cell_type": "markdown", + "id": "a6b6ff9c-d1e2-4422-ab1b-e7f4915ea0f9", + "metadata": {}, + "source": [ + "The previous part showed that linear convolution between $Q'$ and $T$ can be computed via circular convolution when both arrays are zero-padded till they reach the length $N=n+m-1$, where $m$ and $n$ are the length of $Q'$ and $T$, respectively. However, the computation via the circular convolution formula provided in the previous section is not faster. In fact, it has the time complexity of $O(N^2)$! So, why did we go through all that trouble to show that linear convolution can be computed via circular convolution? This is because the circular convolution can be computed more efficiently when it is computed with the help of Fast Fourier Transform (FFT). So, instead of computing the circular convolution in time domain, i.e.\n", + "\n", + "\\begin{equation}\n", + "C_{N}[idx] = \\sum_{i=0}^{N-1} \\tilde{T}_{N}[i] \\times \\tilde{Q'}_{N}[(idx-i)_{N}], \\quad 0 \\le idx \\le N-1\n", + "\\end{equation}\n", + "\n", + "we can compute it by taking it to the frequency domain:\n", + "\n", + "\\begin{equation}\n", + "C_{N} = IFFT\\left(\n", + "FFT(\\tilde{T}_{N})\n", + "\\times\n", + "FFT(\\tilde{Q'}_{N})\n", + "\\right)\n", + "\\end{equation}\n", + "\n", + "and **the time complexity becomes $O(NlogN)$**.\n", + "\n", + "Since I can now compute the circular convolution faster, it means that I can calculate the linear convolution faster! and therefore I can obtain the sliding dot product faster!!" + ] + }, + { + "cell_type": "markdown", + "id": "181cc549-9b2a-48ac-9077-71174cb49f1c", + "metadata": {}, + "source": [ + "## 2.2 Attempt II: Avoid Unnecessary Computation in Circular Convolution\n", + "\n", + "The time complexity of FFT-based circular convolution is $O(NlogN)$, and $N$ is $n+m-1$, the size of output in linear convolution. Can we reduce $N$ itself? \n", + "\n", + "Let's revisit the linear convolution:\n", + "\n", + "\\begin{equation}\n", + "C[idx] = \\sum_{i=0}^{n-1}{T[i] \\times Q'[idx-i]}, \\quad 0 \\le idx \\le n+m-2\n", + "\\end{equation}\n", + "\n", + "\n", + "\n", + "The length of output is $n + m - 1$. However, the sliding dot product (sdp) between $Q$ and $T$ is in range $[m-1, n)$. Therefore, if all we care about is the sliding dot product, the computation can stop at index $n$. In other words, $idx$ (of linear convolution) can be from $0$ to $n-1$. How about the Circular convolution? If, in the Circular Convolution, we accordingly set the length of period $N$ to $n$, does the slice in range $[m-1,n)$ still reflects the values of linear convolution in that same slice? We need to show that\n", + "\n", + "\\begin{equation}\n", + "C[idx] = \\sum_{i=0}^{n-1}{T[i] \\times Q'[idx-i]}, \\quad 0 \\le idx \\le n-1\n", + "\\end{equation}\n", + "\n", + "and\n", + "\n", + "\\begin{equation}\n", + "C_{N=n}[idx] = \\sum_{i=0}^{n-1}T[i] \\times \\tilde{Q'}_{n}[(idx-i)_{n}], \\quad 0 \\le idx \\le n-1\n", + "\\end{equation}\n", + "\n", + "give the same result for the slice $idx \\in \\{m-1,...,n-1\\}$. To prove this, we just need to show that $Q'[idx-i]$ and $\\tilde{Q'}_{n}[(idx-i)_{n}]$ have the same value when $m-1 \\le idx \\le n-1$.
\n", + "\n", + "\n", + "**Proof:** Let $j$ denote $idx-i$. We need to check two different cases:\n", + "\n", + "**Case I: $j \\ge 0$**
\n", + "\n", + "Let's compute the upper bound for $j$. Recall that: \n", + "\n", + "* $j=idx-i$\n", + "* $m-1 \\le idx \\le n-1$\n", + "* $0 \\le i \\le n-1$\n", + "\n", + "\\begin{equation}\n", + "max(j)=max(idx-i)=max(idx)+max(-i)=max(idx) - min(i)=(n-1)-0=n-1\n", + "\\end{equation}\n", + "\n", + "Since $0 \\le j \\le n-1$, $\\tilde{Q'}_{n}[(j)_{n}] = \\tilde{Q'}_{n}[j]$, and this is the same as $Q'[j]$ for the given range of $j$. \n", + "\n", + "Proof is now complete for this case.\n", + "\n", + "**Case II: $j < 0$**
\n", + "Note that $Q'[j]=0$ for $j < 0$. So, we need to show $\\tilde{Q'}_{n}[(j)_{n}]$ becomes zero as well in this case.\n", + "\n", + "Let's compute the lower bound for $j$. Recall that:\n", + "\n", + "* $j=idx-i$\n", + "* $m-1 \\le idx \\le n-1$\n", + "* $0 \\le i \\le n-1$\n", + "\n", + " \n", + "Recall that $\\tilde{Q'}$ is a periodic sequence. Therefore, $\\tilde{Q'}_{n}[(j)_{n}] = \\tilde{Q'}_{n}[(j+n)_{n}]$. Note that: \n", + "\n", + "\\begin{equation}\n", + "j+n \\ge min(j+n) \\ge min(j)+n \\ge min(idx-i)+n \\ge min(idx)+min(-i) + n \\ge min(idx) - max(i) + n \\ge ((m-1) - (n-1)) + n \\ge m\n", + "\\end{equation}\n", + "\n", + "and\n", + "\n", + "\\begin{equation}\n", + "j+n \\le max(j+n) \\le max(j)+n \\le -1+n\n", + "\\end{equation}\n", + "\n", + "So, $j+n$ is in range $[m, n-1]$. Therefore, $j+n$ module $n$ becomes $j+n$. Subsequently, $\\tilde{Q'}_{n}[(j+n)_{n}] = \\tilde{Q'}_{n}[j+n]$. As shown above, The index $j+n$ is at least $m$, meaning it falls into the zero-padding zone of $\\tilde{Q'}_{n}$, and its value is zero. \n", + "\n", + "Proof is now complete for this case. \n", + "\n", + "---" + ] + }, + { + "cell_type": "markdown", + "id": "03cbd90d-a7e4-466d-af15-81d9a1c4ce6c", + "metadata": {}, + "source": [ + "We just showed that the sliding dot product can still be achieved by computing a circular convolution for length $n$. This can be computed in $O(nlogn)$." + ] + }, + { + "cell_type": "markdown", + "id": "125eb500-0452-4b7b-974e-cb17c855b08e", + "metadata": {}, + "source": [ + "## 2.3 Attempt III: Use RFFT/IRFFT instead of FFT/IFFT\n", + "\n", + "In matrix profile / STUMPY, the assumption is that the time series data are all real-valued arrays. As stated in [Scipy's document](https://docs.scipy.org/doc/scipy/reference/generated/scipy.fft.rfft.html):\n", + "\n", + "> When the DFT is computed for purely real input, the output is Hermitian-symmetric, i.e., the negative frequency terms are just the complex conjugates of the corresponding positive-frequency terms.\n", + "\n", + "RFFT/IRFFT is a certain type of FFT/IFFT that leverages this property to perform faster fourier transform. It turns out that the FFT and IFFT in Eq. (18) can be replaced with RFFT and IRFFT when inputs are real-valued arrays. This should improve the performance of computing the circular convolution when input sequences have only real values." + ] + }, + { + "cell_type": "markdown", + "id": "97fcc328-35f9-4b66-9598-1b5d8f1f97ac", + "metadata": {}, + "source": [ + "# 3. How to compute convolution faster when $T$ is long?\n", + "\n", + "In this section, we learn about a divide-and-conquer approach in convolution. We first need to learn about the \"linearity\" property in linear convolution. Then, we learn how it allows us to divide the problem into similar problems. We then compute circular convolution on those smaller problems and combine their outputs." + ] + }, + { + "cell_type": "markdown", + "id": "1e45ccf7-1564-4c41-b32e-4de9b1b201a3", + "metadata": {}, + "source": [ + "## 3.1 The \"Linearity\" property in Linear Convolution" + ] + }, + { + "cell_type": "markdown", + "id": "a084578e-cb20-437c-a70c-3e33f8da8ac3", + "metadata": {}, + "source": [ + "Suppose the array $T$ can be written as $T^{(1)} + T^{(2)}$. In other words: $T[i] = T^{(1)}[i] + T^{(2)}[i]$, then:\n", + "\n", + "\\begin{equation}\n", + "C[idx] = \\sum_{i=-\\infty}^{\\infty}{T[i] \\times Q'[idx-i]}\n", + "\\end{equation}\n", + "\n", + "\\begin{equation}\n", + "C[idx] = \\sum_{i=-\\infty}^{\\infty}{(T^{(1)}[i] \\times Q'[idx-i] + T^{(2)}[i] \\times Q'[idx-i])}\n", + "\\end{equation}\n", + "\n", + "\\begin{equation}\n", + "C[idx] = \\sum_{i=-\\infty}^{\\infty}{T^{(1)}[i] \\times Q'[idx-i]} + \\sum_{i=-\\infty}^{\\infty}{T^{(2)}[i] \\times Q'[idx-i]}\n", + "\\end{equation}\n", + "\n", + "\\begin{equation}\n", + "C[idx] = C^{(1)}[idx] + C^{(2)}[idx]\n", + "\\end{equation}" + ] + }, + { + "cell_type": "markdown", + "id": "1f5be13b-82a7-46d9-9cb1-229c878e67f8", + "metadata": {}, + "source": [ + "This shows that a linear convolution has the property of \"linearity\". Now, we use this property to show that we can compute the linear convolution between $Q'$ and $T$ by breaking $T$ into smaller parts." + ] + }, + { + "cell_type": "markdown", + "id": "fecee5fc-3c10-4c2e-8c66-bce17c3650bd", + "metadata": {}, + "source": [ + "## 3.2 The overlap-add method \n", + "\n", + "To understand this method, we first learn about its core logic. Then, we will discuss its efficiency." + ] + }, + { + "cell_type": "markdown", + "id": "f5711673-578d-4a21-a183-7abf510f0ca2", + "metadata": {}, + "source": [ + "### 3.2.1 The overlap-add method: Core logic\n", + "\n", + "Suppose $T$ is $\\{t_{0}, t_{1}, ..., t_{n-1}\\}$. I can create two sequences $T^{(1)}$ and $T^{(2)}$ such that they both have the same length as $T$, and they have the following values:\n", + "\n", + "* $T^{(1)}$: $\\{t_{0}, t_{1}, ..., t_{s-1}, 0, ..., 0\\}$\n", + "* $T^{(2)}$: $\\{0, 0, ..., 0, t_{s}, t_{s+1}, ..., t_{n-1}\\}$\n", + "\n", + "Note that, at any index $i$, $T[i] = T^{(1)}[i] + T^{(2)}[i]$." + ] + }, + { + "cell_type": "markdown", + "id": "d033e311-81c1-4613-b5ba-be43c3c0c24c", + "metadata": {}, + "source": [ + "We can compute the linear convolution between each new sequence and $Q'$:\n", + "\n", + "* $C^{(1)}$: Linear convolution between $Q'$ and $T^{(1)}$. $C^{(1)}$ is an array with length $n+m-1$, and $C^{(1)} = \\{c_{0}^{(1)}, c_{1}^{(1)}, ..., c_{s+m-2}^{(1)}, 0, ..., 0\\}$\n", + "* $C^{(2)}$: Linear convolution between $Q'$ and $T^{(2)}$. $C^{(2)}$ is an array with length $n+m-1$, and $C^{(2)} = \\{0, ..., 0, c_{s}^{(2)}, ..., c_{n+m-2}^{(2)}\\}$" + ] + }, + { + "cell_type": "markdown", + "id": "a014c77d-fdb1-4f42-ab13-cd6af26ae451", + "metadata": {}, + "source": [ + "The elements of $C$, the convolution between $T$ and $Q'$, can be computed by summing $C^{(1)}$ and $C^{(2)}$:\n", + "\n", + "\n", + "* $c_{0} = c^{(1)}_{0} + 0$\n", + "* $c_{1} = c^{(1)}_{1} + 0$\n", + "* ...\n", + "* $c_{s-1} = c^{(1)}_{s-1} + 0$\n", + "* $\\textcolor{red}{c_{s} = c^{(1)}_{s} + c^{(2)}_{s}}$\n", + "* ...\n", + "* $\\textcolor{red}{c_{s+m-2} = c^{(1)}_{s+m-2} + c^{(2)}_{s+m-2}}$\n", + "* $c_{s+m-1} = 0 + c^{(2)}_{s+m-1}$\n", + "* $c_{s+m} = 0 + c^{(2)}_{s+m}$\n", + "* ...\n", + "* $c_{n+m-2} = 0 + c^{(2)}_{n+m-2}$" + ] + }, + { + "cell_type": "markdown", + "id": "f5e2d3b1-b653-441c-a9a2-49f4677464f9", + "metadata": {}, + "source": [ + "As shown above, each of the elements in $\\{\\textcolor{red}{c_{s}, ..., c_{s+m-2}}\\}$ requires elements from both $C^{(1)}$ and $C^{(2)}$. Although the (non-zero) values of $T^{(1)}$ stop at index $s$, its corresponding linear convolution $C^{(1)}$ has values till $s+m-2$. Therefore, it will have some overlaps with $C^{(2)}$, whose (non-zero) values start right from index $s$. This is why the method is called \"overlap-add\"." + ] + }, + { + "cell_type": "markdown", + "id": "44608755-9b19-418a-b690-37da32ec4365", + "metadata": {}, + "source": [ + "Note that we have not improved the efficiency yet. In fact, we are computing the linear convolution twice without reducing the size of problem. Note that some of elements of $C^{(1)}$ and $C^{(2)}$ are zero. If we can avoid computing those elements, then we might have a chance to improve the efficiency! \n", + "\n", + "Let's revisit the input arrays:\n", + "\n", + "* $T^{(1)}$: $\\{t_{0}, t_{1}, ..., t_{s-1}, 0, ..., 0\\}$\n", + "* $T^{(2)}$: $\\{0, 0, ..., 0, t_{s}, t_{s+1}, ..., t_{n-1}\\}$" + ] + }, + { + "cell_type": "markdown", + "id": "a6ac56c0-4516-497f-961b-d62acf74f87e", + "metadata": {}, + "source": [ + "For each case, we know that the length of its corresponding linear convolution with $Q'$ is $n+m-1$. We also know that some elements of convolution become zero. **Therefore, we can just compute the convolution for non-zero elements and include its contribution to its corresponding slice in $C$**. \n", + "\n", + "So, we can initialize an empty array $C$ of length $n + m - 1$, and fill it with zeros. Then:\n", + "\n", + "1. Compute the linear convolution between $Q'$ and $\\{t_{0}, t_{1}, ..., t_{s-1}\\}$. Its output will have length $L=s + m - 1$. We can use it to update $C[0:0+L]$.\n", + "2. Compute the linear convolution between $Q'$ and $\\{t_{s}, t_{1}, ..., t_{n-1}\\}$. Its output will have length $L=(n-s) + m - 1$. We can use it to update $C[s:s+L]$. Note that the index is shifted by $s$ as the second chunk's contribution starts from index $s$." + ] + }, + { + "cell_type": "markdown", + "id": "0a43e23b-d913-4751-9631-3b318bfe16ed", + "metadata": {}, + "source": [ + "This allows us to avoid computing circular convolution for length $n + m - 1$ in both cases. As observed, we compute it for lengths \"$s + m - 1$\" and \"$(n-s) + m - 1$\" instead. We are improving the performance of calculation for each chunk as we are reducing the length of circular convolution. Note that we have not discussed the overall efficiency yet. " + ] + }, + { + "cell_type": "markdown", + "id": "a2f37dba-0dd5-45d9-aa78-243f5b5909c8", + "metadata": {}, + "source": [ + "**Where does the name \"overlap-add\" come from?**\n", + "\n", + "Let's revisit the example above. The first chunk of `T` contributes to the `slice(0, L)` of the output and the second chunk of `T` contributes to the `slice(s, s+L)` of the output, where `s` is the point of split in `T` and `L=s+m-1`. So, the first slice stops at `L` and the second slice starts at `s < L`. Hence, there is an overlap. So, the term \"overlap\" refers to the overlap between the slices of outputs that correspond to the adjacent, non-overlapping chunks of input. The term \"add\" refers to the fact that the overlapped parts are summed together." + ] + }, + { + "cell_type": "markdown", + "id": "081d6fa7-fc57-41e8-89d1-299bfabb1e92", + "metadata": {}, + "source": [ + "**Example:**
\n", + "The following animation shows the computation of a full linear convolution between `T=[1, 2, ..., 12]` and `Qr=[3, 2, 1]` via overlap-add approach. In this example, `T` is split into two non-overlapping chunks `T1=[1,2,3,4]` and `T2=[5, ..., 12]`." + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "id": "b1abff79-8a87-4013-a769-923a2a430ecf", + "metadata": {}, + "outputs": [ + { + "data": { + "text/html": [ + "" + ], + "text/plain": [ + "" + ] + }, + "execution_count": 3, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "with open(\"overlap_add_convolution.gif\", \"rb\") as f:\n", + " data = base64.b64encode(f.read()).decode()\n", + "\n", + "HTML(f'')" + ] + }, + { + "cell_type": "markdown", + "id": "2c410e5b-17f4-41b0-8790-9b61273e003f", + "metadata": {}, + "source": [ + "### 3.2.2 The overlap-add method: Efficiency" + ] + }, + { + "cell_type": "markdown", + "id": "c34a99c1-fc92-4fb9-a5b1-3e1776e915bd", + "metadata": {}, + "source": [ + "In the previous section, we showed that the linear convolution between $Q'$ and $T$ can be computed by breaking $T$ into two chunks and performing linear convolution on each. We also learned that we do not have to compute convolution for the full length $n+m-1$ each time, and instead, we can do it on smaller size and just update its corresponding slice in the convolution $C$. Recall that the linear convolution can be computed faster with FFT-based circular convolution. \n", + "\n", + "So, to compute the linear convolution in the previous example via overlap-add method, we need to:\n", + "\n", + "1. Compute linear convolution between $Q'$ and $\\{t_{0}, t_{1}, ..., t_{s-1}\\}$: This can be computed via FFT-based circular convolution with length $N = s + m - 1$. This will contribute to the slice $C[0:0+N]$\n", + "2. Compute linear convolution between $Q'$ and $\\{t_{s}, t_{1}, ..., t_{n-1}\\}$: This can be computed via FFT-based circular convolution with length $N = (n-s) + m - 1$. This will contribute to the slice $C[s:s+N]$. Note that the chunk of $T$ starts from index $s$ in this case. Therefore, its contribution to convolution also starts from $s$." + ] + }, + { + "cell_type": "markdown", + "id": "520fd6db-9c41-4b2c-8f14-3b4718d7bd46", + "metadata": {}, + "source": [ + "The heavy part of computation is coming from the Fourier Transforms:\n", + "\n", + "1. Circular convolution that corresponds to first chunk of $T$\n", + " 1. FFT of $Q'$, padded with 0s till it reaches the length $s + m - 1$\n", + " 2. FFT of $T[:s]$, padded with 0s till it reaches the length $s + m - 1$\n", + " 3. element-wise product of two FFTs, and compute its IFFT\n", + "2. Circular convolution that corresponds to second chunk of $T$\n", + " 1. FFT of $Q'$, padded with 0s till it reaches the length $(n-s) + m - 1$\n", + " 2. FFT of $T[s:]$, padded with 0s till it reaches the length $(n-s) + m - 1$\n", + " 3. element-wise product of two FFTs, and compute its IFFT" + ] + }, + { + "cell_type": "markdown", + "id": "9dbcff71-8c1e-4538-a9a5-7939302ee06c", + "metadata": {}, + "source": [ + "What happens if $n = 2 \\times s$? In this case, the length of circular convoluton in the second case, i.e. $(n-s) + m - 1$, becomes equal to $s + m - 1$, the length of the first convolution. This allows us to:\n", + "* Avoid computing the FFT of zero-padded $Q'$ for the second case\n", + "* Once the coefficients of FFT are computed, they can be reused to compute FFT for other arrays that have the same length of $s + m - 1$\n", + "* Once the coefficients of an IFFT are computed, they can be reused to compute IFFT for other arrays that have the same length of $s + m - 1$" + ] + }, + { + "cell_type": "markdown", + "id": "4874492e-1057-412b-8219-58cbbc4de078", + "metadata": {}, + "source": [ + "Now, let's see how it looks like when $n$ has multiple chunks of size $s$. Suppose $k = n // s$, and $r = n - ks$. Let's write down the chunks and their contribution to the convolution $C$, whose length is $n + m - 1$." + ] + }, + { + "cell_type": "markdown", + "id": "de128f43-f2a6-42f5-b52b-a6bdfdd9db20", + "metadata": {}, + "source": [ + "| T
chunk start | T
chunk stop | length of linear convolution
between Q' and the chunk of T
( = input size for circular convolution) | Contribution to C
slice start | Contribution to C
slice stop |\n", + "| ----------- | ---------- | ---------------------------------------- | ------------------------ | ----------------------- |\n", + "| 0 | s | s + m -1 | 0 | s + m - 1 |\n", + "| s | 2s | s + m - 1 | s | s + (s + m - 1) |\n", + "| 2s | 3s | s + m - 1 | 2s | 2s + (s + m - 1) |\n", + "| ... | | | | |\n", + "| ... | | | | |\n", + "| ... | | | | |\n", + "| (k-1)s | ks | s + m - 1 | (k-1)s | (k-1)s + (s + m - 1) |\n", + "| ks | n | r + m - 1 | ks | ks + (r + m - 1) = n + m - 1 |" + ] + }, + { + "cell_type": "markdown", + "id": "e27d0462-0a1a-45f9-9942-bca08350d13d", + "metadata": {}, + "source": [ + "A few things to notice here:\n", + "\n", + "* Last chunk might have a size that is less than $s$, and hence the length of input for circular convolution becomes different. However, we can still zero-pad the inputs with zeros till they reach the size $s + m - 1$, which is the same as other chunks. In such case, the values up to index $r + m - 1$ (exclusive) will be considered from the output, where $r = n \\% s$\n", + "* The contribution of last slice stops at index $n+m-1$, which is aligned with the length of $C$, the linear convolution between $Q'$ and full $T$." + ] + }, + { + "cell_type": "markdown", + "id": "b6950aa2-21f3-4797-bc42-85814b41b8b9", + "metadata": {}, + "source": [ + "# 4. Conclusion" + ] + }, + { + "cell_type": "markdown", + "id": "099a313a-31dc-4677-8a25-4cd46ea2aa72", + "metadata": {}, + "source": [ + "We started with the sliding dot product (sdp) between a query $Q$ and a time series $T$, and saw that it is a slice of the linear convolution between $Q'$ (the reverse of $Q$) and $T$. Two key ideas make computing this convolution fast, even when $T$ is very long.\n", + "\n", + "## 4.1 Linear convolution can be computed via circular convolution\n", + "\n", + "Let $Q'$ and $T$ have lengths $m$ and $n$, respectively ($m < n$). Their full linear convolution $C$ has $n + m - 1$ values, and we can obtain it from a circular convolution as long as the period is large enough to hold all of them. So, we zero-pad both $Q'$ and $T$ until their length becomes\n", + "\n", + "\\begin{equation}\n", + "N = n + m - 1\n", + "\\end{equation}\n", + "\n", + "and compute their circular convolution. The circular convolution is identical to the full linear convolution. This is useful because a circular convolution can be computed with FFT in $O(N\\log N)$, instead of $O(N^2)$:\n", + "\n", + "\\begin{equation}\n", + "C = IFFT\\left(FFT(\\tilde{T}_{N}) \\times FFT(\\tilde{Q'}_{N})\\right)\n", + "\\end{equation}\n", + "\n", + "where $\\tilde{T}_{N}$ and $\\tilde{Q'}_{N}$ are the zero-padded versions of $T$ and $Q'$.\n", + "\n", + "If we only care about the **valid** portion of the convolution, which is exactly the sdp, we can do even better and use a shorter period:\n", + "\n", + "\\begin{equation}\n", + "N = \\max(m, n) = n\n", + "\\end{equation}\n", + "\n", + "The slice $[m-1, n)$ is still exactly the valid convolution, which is equivalent to sdp, the sliding dot product between $Q$ and $T$.\n", + "\n", + "## 4.2 Overlap-add: divide and conquer for long $T$\n", + "\n", + "When $T$ is long, we do not need to convolve it with $Q'$ in one shot. Because convolution is linear, we can:\n", + "\n", + "1. **Divide:** split $T$ into non-overlapping chunks. The chunks can have any sizes.\n", + "2. **Conquer:** compute the full linear convolution between $Q'$ and each chunk separately. Each of these can be computed with FFT-based circular convolution of length $N = len(chunk) + m - 1$.\n", + "3. **Combine:** add the output of each chunk into the output $C$. A chunk that starts at index $start$ contributes to the slice of $C$ that starts at $start$. The slices of the output that correspond to the neighboring chunks overlap by $m-1$ elements, and the overlapped elements are summed. This is where the name \"overlap-add\" comes from. So, again, there is no overlap in the chunks of $T$ but the overlap is between the slices of the output that correspond to those chunks.\n", + "\n", + "The final result is identical to the full linear convolution between $Q'$ and the whole $T$. When all chunks have the same size $s$, every circular convolution has the same length $N = s + m - 1$. In that case, the FFT of the zero-padded $Q'$ can be computed only once and reused for all chunks. Furthermore, the twiddle factors for forward fourier transform can be computed only once since it depends on the size of input only. Same is true for backward fourier transform.\n", + "\n", + "#### **In a nutshell:**\n", + "* A full/valid linear convolution between $T$ (with length of $n$) and $Q'$ (with length of $m$) can be computed via **FFT-based circular convolution**. For full (which contains the valid slice as well), set $N=n+m-1$. For valid only, you can optimize further by setting $N=n$.\n", + "* A full linear convolution can be computed via **Overlap-add**. A long $T$ is broken into chunks. The full linear convolution between each chunk and $Q'$ is computed, and the (overlapping) outputs are added up. A valid linear convolution can then be obtained by selecting the slice(m-1, n) from the final output." + ] + }, + { + "cell_type": "markdown", + "id": "d7ba4836-7b72-4ad1-a4c9-dd92879843a9", + "metadata": {}, + "source": [ + "# Bonus: Code" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "id": "d28a2641-98c9-46f4-9052-ad429089ddfb", + "metadata": {}, + "outputs": [], + "source": [ + "# imports\n", + "import functools\n", + "import numpy as np" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "id": "f7886382-d2a8-4fee-bc6c-3d4a09a4e375", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "seed: 155504\n", + "m: 250\n", + "n: 461\n" + ] + } + ], + "source": [ + "# inputs\n", + "seed = np.random.randint(0, 1000000)\n", + "np.random.seed(seed)\n", + "m, n = np.sort(np.random.choice(np.arange(1, 1000), size=2))\n", + "m = int(m)\n", + "n = int(n)\n", + "\n", + "Q = np.random.rand(m)\n", + "T = np.random.rand(n)\n", + "\n", + "ref_full = np.convolve(Q, T, mode='full')\n", + "ref_valid = np.convolve(Q, T, mode='valid')\n", + "\n", + "valid_slice = slice(m-1, n) # Section 1.2\n", + "np.testing.assert_allclose(ref_valid, ref_full[valid_slice])\n", + "\n", + "print('seed: ', seed)\n", + "print('m: ', m)\n", + "print('n: ', n)" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "id": "0dd21bfc-902c-4cfc-b9dd-ba69f6e44b31", + "metadata": {}, + "outputs": [], + "source": [ + "# Section 2.1\n", + "# linear convolution can be computed via fft-based circular convolution\n", + "def fft_circular_convolution(Q, T, N=None):\n", + " \"\"\"\n", + " Returns circular convolution between padded Q and padded T, where each\n", + " is padded with 0s till they reach the length of `N`. When `N` is None, \n", + " it is default to `n + m - 1` and the output is equivalent to\n", + " full linear convolution.\n", + " \"\"\"\n", + " m = len(Q)\n", + " n = len(T)\n", + " if N is None: \n", + " # default to computing \n", + " # full linear convolution\n", + " # between Q and T\n", + " N = n + m - 1\n", + "\n", + " # pad Q and T till they reach length N\n", + " padded_Q = np.zeros(N, dtype=np.float64)\n", + " padded_T = np.zeros(N, dtype=np.float64)\n", + "\n", + " padded_Q[:m] = Q\n", + " padded_T[:n] = T\n", + "\n", + " # N-circular convolution between padded_Q and padded_T\n", + " out = np.fft.ifft(\n", + " np.fft.fft(padded_Q) * np.fft.fft(padded_T)\n", + " ).real\n", + "\n", + " return out\n", + "\n", + "# full linear convolution\n", + "comp = fft_circular_convolution(Q, T, N=n+m-1)\n", + "np.testing.assert_allclose(comp, ref_full)" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "id": "b50f83f7-838a-4abc-88cb-9379da270445", + "metadata": {}, + "outputs": [], + "source": [ + "# Section 2.2. \n", + "# Avoid unnecessary computation When we care about 'valid' convolution\n", + "\n", + "# valid linear convolution\n", + "comp = fft_circular_convolution(Q, T, N=n)[valid_slice] # Note that setting N=n is ok for `valid`\n", + "np.testing.assert_allclose(comp, ref_valid)" + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "id": "62b8fb3a-3b4a-467c-b9d9-cdff2c2b5dd6", + "metadata": {}, + "outputs": [], + "source": [ + "# Section 2.3 \n", + "# Use RFFT/IRFFT instead of FFT/IFFT\n", + "def rfft_circular_convolution(Q, T, N=None):\n", + " \"\"\"\n", + " Returns circular convolution between padded Q and padded T, where each\n", + " is padded with 0s till they reach the length of `N`. When `N` is None, \n", + " it is default to `n + m - 1` and the output is equivalent to\n", + " full linear convolution.\n", + " \"\"\"\n", + " m = len(Q)\n", + " n = len(T)\n", + " if N is None:\n", + " N = m + n - 1\n", + " \n", + " padded_Q = np.empty(N, dtype=np.float64)\n", + " padded_Q[:m] = Q\n", + " padded_Q[m:] = 0.0\n", + "\n", + " padded_T = np.empty(N, dtype=np.float64)\n", + " padded_T[:n] = T\n", + " padded_T[n:] = 0.0\n", + "\n", + " out = np.fft.irfft(\n", + " np.fft.rfft(padded_Q) * np.fft.rfft(padded_T), n=N\n", + " ).real\n", + "\n", + " return out\n", + "\n", + "# test1: full linear convolution \n", + "comp = rfft_circular_convolution(Q, T, N=n+m-1)\n", + "np.testing.assert_allclose(comp, ref_full)\n", + "\n", + "# test2: valid linear convolution\n", + "comp = rfft_circular_convolution(Q, T, N=n)[valid_slice] # note that we are fine with N=n\n", + "np.testing.assert_allclose(comp, ref_valid)" + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "id": "540113b0-a8b4-41a5-aaa9-6eedb18ee31a", + "metadata": {}, + "outputs": [], + "source": [ + "# Section 3.2.1\n", + "# overlap-add method\n", + "def _create_random_chunks(n):\n", + " \"\"\"\n", + " Returns a sorted list of tuples such that each tuple (start, stop) is\n", + " the start and stop of a slice. No two slices have overlap, and the union\n", + " of all slices cover the range(0, n). \n", + " \"\"\"\n", + " n_cuts = np.random.randint(low=0, high=n, size=1)\n", + " IDX = [0] + sorted(np.random.choice(np.arange(1, n), size=n_cuts, replace=False)) + [n]\n", + " T_chunks_range = []\n", + " for i in range(len(IDX)-1):\n", + " T_chunks_range.append((IDX[i], IDX[i+1]))\n", + "\n", + " return T_chunks_range\n", + "\n", + "\n", + "def overlap_add_convolve(Q, T, T_chunks_range, full_conv_func=None):\n", + " \"\"\"\n", + " Returns the linear convolution between Q and T.\n", + " It gives the same result as np.convolve(Q, T, mode='full')\n", + "\n", + " Parameters\n", + " ----------\n", + " Q : numpy.ndarray\n", + " A 1D array representing a query data\n", + "\n", + " T : numpy.ndarray\n", + " A 1D array representing a time series data\n", + "\n", + " T_chunks_range : list\n", + " A list of tuples, such that each tuple shows\n", + " start and stop index of a slice of T. The ranges\n", + " provided in the tuples have no overlap with each\n", + " other, and the ranges cover the full T\n", + "\n", + " full_conv_func : callable\n", + " A callable object that computes full linear convolution between\n", + " two arrays. When None (default), it is automatically set to\n", + " `np.convolve(..., mode='full')`\n", + "\n", + " Returns\n", + " -------\n", + " The full linear convolution between Q and T\n", + " \"\"\"\n", + " if full_conv_func is None:\n", + " full_conv_func = functools.partial(np.convolve, mode='full')\n", + " \n", + " m = len(Q)\n", + " n = len(T)\n", + "\n", + " # compute full linear convolution\n", + " # by leveraging the \"linearity property\" (Section 3.1)\n", + " out = np.zeros(n + m - 1, dtype=np.float64)\n", + " for (start, stop) in T_chunks_range: \n", + " # full linear conv between Q and T[start:stop] contributes\n", + " # to a certain slice in output:\n", + " output_slice_start = start\n", + " output_slice_stop = output_slice_start + ((stop - start) + m - 1)\n", + " \n", + " # add contribution\n", + " out[output_slice_start: output_slice_stop] += full_conv_func(Q, T[start:stop]) \n", + "\n", + " return out\n", + "\n", + "\n", + "# Note: \n", + "# In all of the following tests, full_conv_func gives full linear convolution\n", + "\n", + "T_chunks_range = _create_random_chunks(n)\n", + "\n", + "# test1\n", + "# directly use linear convolution to compute full linear convolution between each chunk and Q\n", + "comp = overlap_add_convolve(Q, T, T_chunks_range, full_conv_func=np.convolve) \n", + "np.testing.assert_allclose(comp, ref_full)\n", + "\n", + "# test2\n", + "# use fft-based circular convolution to compute full linear convolution between each chunk and Q\n", + "comp = overlap_add_convolve(Q, T, T_chunks_range, full_conv_func=fft_circular_convolution)\n", + "np.testing.assert_allclose(comp, ref_full)\n", + "\n", + "# test3\n", + "# use rfft-based circular convolution to compute full linear convolution between each chunk and Q\n", + "comp = overlap_add_convolve(Q, T, T_chunks_range, full_conv_func=rfft_circular_convolution)\n", + "np.testing.assert_allclose(comp, ref_full)\n", + "\n", + "# test4\n", + "# Similar to test3 but for mode='valid'\n", + "comp = overlap_add_convolve(Q, T, T_chunks_range, full_conv_func=rfft_circular_convolution)[valid_slice]\n", + "np.testing.assert_allclose(comp, ref_valid)" + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "id": "5f9c97b8-a76b-4430-81f8-ef12bf12f3d3", + "metadata": {}, + "outputs": [], + "source": [ + "# Section 3.2.2\n", + "# Choose fix-sized chunks to\n", + "# (i) avoid recomputing fft for Q \n", + "# (ii) compute twiddle factors only once (hopefully) for forward transform and once for backward transform\n", + "\n", + "def overlap_add_rfftconvolve(Q, T, s):\n", + " \"\"\"\n", + " Returns the full linear convolution between Q and T, computed via the\n", + " overlap-add method with fixed-size chunks of T and rfft-based circular\n", + " convolution. It gives the same result as np.convolve(Q, T, mode='full')\n", + "\n", + " Parameters\n", + " ----------\n", + " Q : numpy.ndarray\n", + " A 1D array representing a query data\n", + "\n", + " T : numpy.ndarray\n", + " A 1D array representing a time series data\n", + "\n", + " s : int\n", + " The size of each chunk of T. T is split into non-overlapping chunks\n", + " of length `s`. If `len(T)` is not a multiple of `s`, the last chunk\n", + " is zero-padded to length `s`. Each chunk is convolved with Q via an\n", + " rfft-based circular convolution of length `N = s + m - 1`, where\n", + " `m = len(Q)`. Since all chunks have the same size, the rfft of the\n", + " zero-padded Q is computed only once and reused for all chunks.\n", + "\n", + " Returns\n", + " -------\n", + " numpy.ndarray\n", + " The full linear convolution between Q and T, with length\n", + " `len(T) + len(Q) - 1`\n", + " \"\"\" \n", + " m = len(Q)\n", + " n = len(T)\n", + " n_chunks_T = int(np.ceil(n/s))\n", + "\n", + " # Ensure T can be splitted into chunks of size `s` \n", + " n_new = s * n_chunks_T\n", + " T_new = np.empty(n_new, dtype=np.float64)\n", + " T_new[:n] = T\n", + " T_new[n:] = 0.0\n", + " \n", + " # Compute circular convolution between Q and each chunk of T\n", + " N = s + m - 1\n", + " padded_T_chunks = np.pad(\n", + " T_new.reshape(-1, s), \n", + " ((0,0), (0, N-s))\n", + " )\n", + " padded_Q = np.pad(Q, (0, N-m))\n", + " \n", + " RFFT_padded_T_chunks = np.fft.rfft(padded_T_chunks, axis=1) # fft on 2D array\n", + " RFFT_padded_Q = np.fft.rfft(padded_Q) # only once\n", + "\n", + " # Since `T` was extended initially, its new length \n", + " # is considered when computing length of output\n", + " out = np.zeros(n_new + m - 1, dtype=np.float64)\n", + " for i in range(n_chunks_T):\n", + " output_slice_start = i * s # same as start of chunk\n", + " output_slice_stop = output_slice_start + N\n", + " out[output_slice_start: output_slice_stop] += np.fft.irfft(RFFT_padded_T_chunks[i] * RFFT_padded_Q, n=N).real\n", + "\n", + " # Note1: We can also take IRFFT on 2D array once and then use each of its row in the for-loop\n", + " # Note2: We can try to vectorize the summation part as well\n", + " # Note3: Optiomal length of `s` can be computed empirically or mathematically (but approximately)\n", + "\n", + " return out[:(n+m-1)]\n", + "\n", + "\n", + "# test full\n", + "s_random = np.random.randint(1, n)\n", + "for s in [1, 2, m//2, m-1, m, m+1, n//2, n-1, s_random]:\n", + " comp = overlap_add_rfftconvolve(\n", + " Q, T, s=s\n", + " )\n", + " np.testing.assert_allclose(comp, ref_full)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "bbf7851f-b734-4446-987a-aa7b91fb2ec5", + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3 (ipykernel)", + "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.14.4" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +}