diff --git a/docs/overlap_add_math.ipynb b/docs/overlap_add_math.ipynb new file mode 100644 index 0000000..083953d --- /dev/null +++ b/docs/overlap_add_math.ipynb @@ -0,0 +1,742 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "8330f3d0-1e58-4b8f-9619-f0baa21a4f16", + "metadata": {}, + "source": [ + "# **Compute Sliding Dot Product using Convolution**" + ] + }, + { + "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 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 to compute convolution faster when $T$ is very long?** This is where we learn about 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", + "$$Q'[j] = Q[m-1-j], \\quad 0 \\le j \\le m-1,$$" + ] + }, + { + "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": [ + "$$ C[idx] = \\sum_{i=-\\infty}^{\\infty}{T[i] \\times Q'[idx-i]}$$" + ] + }, + { + "cell_type": "markdown", + "id": "8028ee30-3244-4c79-86ca-22e2b8876694", + "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.\n", + "\n", + "\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", + "\n", + "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", + "$$ C[idx] = \\sum_{i=0}^{n-1}{T[i] \\times Q'[idx-i]}$$" + ] + }, + { + "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", + "$$\\bigcup_{i=0}^{n-1} [\\,i,\\ i+m-1\\,]$$\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 **exist 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 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", + "$$ C[idx] = \\sum_{i=0}^{n-1}{T[i] \\times Q'[idx-i]}$$\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", + "> **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 in range [m-1,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", + "$$\\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", + "\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", + "$$\\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", + "\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", + "\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", + "$$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", + "\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", + "$$C_{N}[idx] = \\sum_{i=0}^{n-1} T[i] \\times \\tilde{Q'}_{N}[(idx-i)_{N}]$$\n", + "\n", + "\n", + "\n", + "So, instead of the original equations, we need to show the following ones are equivalent:\n", + "\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", + "\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", + "\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", + "$$max(j) = max(idx-i) = max(idx) + max(-i) = max(idx) - min(i) = (N-1) - 0 = N-1$$ \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", + "$$min(j) = min(idx-i) = min(idx) + min(-i) = min(idx) - max(i) = 0 - (n-1)$$\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", + "$$ -(n-1)+N \\le j+N \\le -1+N $$\n", + "\n", + "Recall that $N=n+m-1$. Let's plug that value into $-(n-1)+N$. The inequality becomes:\n", + "\n", + "$$(n+m-1)-(n-1) \\le j+N \\le N-1 \\implies m \\le j+N \\le N-1$$\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", + "$$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", + "\n", + "we can compute it by taking it to the frequency domain:\n", + "\n", + "$$C_{N} = IFFT\\left(\n", + "FFT(\\tilde{T}_{N})\n", + "\\times\n", + "FFT(\\tilde{Q'}_{N})\n", + "\\right)$$\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", + "$$ C[idx] = \\sum_{i=0}^{n-1}{T[i] \\times Q'[idx-i]}, \\quad 0 \\le idx \\le n+m-2$$\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", + "$$ C[idx] = \\sum_{i=0}^{n-1}{T[i] \\times Q'[idx-i]}, \\quad 0 \\le idx \\le n-1$$\n", + "\n", + "and\n", + "\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", + "\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", + "$$ max(j)=max(idx-i)=max(idx)+max(-i)=max(idx) - min(i)=(n-1)-0=n-1$$\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", + "$$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", + "$$j+n \\le max(j+n) \\le max(j)+n \\le -1+n$$\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 RFFT 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", + "$$ C[idx] = \\sum_{i=-\\infty}^{\\infty}{T[i] \\times Q'[idx-i]} $$\n", + "\n", + "$$ C[idx] = \\sum_{i=-\\infty}^{\\infty}{(T^{(1)}[i] \\times Q'[idx-i] + T^{(2)}[i] \\times Q'[idx-i])} $$\n", + "\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", + "\n", + "$$ C[idx] = C^{(1)}[idx] + C^{(2)}[idx] $$" + ] + }, + { + "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$ twice. Instead, we compute it for lengths \"$s + m - 1$\" and \"$(n-s) + m - 1$\". 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": "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.\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": "d7ba4836-7b72-4ad1-a4c9-dd92879843a9", + "metadata": {}, + "source": [ + "**Let's see this in action!**" + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "id": "c042721c-dc3c-4e8d-a10b-0f7fcd07be8d", + "metadata": {}, + "outputs": [], + "source": [ + "import numpy as np\n", + "import numba\n", + "\n", + "\n", + "def direct_convolve(Q, T):\n", + " # linear convolution in one shot\n", + " return np.convolve(Q, T, mode='full')\n", + " \n", + "\n", + "def overlap_add_convolve(Q, T, s=None):\n", + " if s is None:\n", + " s = 2 * len(Q)\n", + " \n", + " m = len(Q)\n", + " n = len(T)\n", + "\n", + " # break T into chunks with size s,\n", + " # compute linear convolution between each and Q\n", + " # and add it to its corresponding slice in output\n", + " out = np.zeros(n + m - 1, dtype=np.float64)\n", + " \n", + " n_chunks = int(np.ceil(n / s))\n", + " for i in range(n_chunks):\n", + " start_chunk = i * s\n", + " stop_chunk = min(start_chunk + s, n)\n", + " T_chunk = T[start_chunk:stop_chunk]\n", + "\n", + " l_conv = len(T_chunk) + m - 1\n", + " start_conv = start_chunk\n", + " stop_conv = start_conv + l_conv\n", + "\n", + " out[start_conv:stop_conv] += np.convolve(Q, T_chunk, mode='full')\n", + "\n", + " return out\n", + "\n", + "\n", + "# input\n", + "T = np.random.rand(1234)\n", + "Q = np.random.rand(123)\n", + "\n", + "ref = direct_convolve(Q, T)\n", + "\n", + "s = 2 * len(Q) # overlap-add\n", + "comp = overlap_add_convolve(Q, T, s=2*len(Q))\n", + "\n", + "np.testing.assert_allclose(ref, comp)" + ] + }, + { + "cell_type": "markdown", + "id": "26748568-b6d5-4852-9a5d-aa11102562be", + "metadata": {}, + "source": [ + "To improve the efficiency, different avenues can be explored:\n", + "* Use numba's prange to parallelize the computation\n", + "* Compute cost of operation and find an optimal value for the chunk size $s$" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "10b4039b-1f79-40b7-a30c-3fd67ae654b4", + "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 +}