From 65bdf79a8fd5d039b28a6490ed365288c5bb8e00 Mon Sep 17 00:00:00 2001 From: NimaSarajpoor Date: Tue, 1 Sep 2026 23:51:54 -0400 Subject: [PATCH 1/8] added notebook to explain math ffor convolution --- docs/overlap_add_math.ipynb | 222 ++++++++++++++++++++++++++++++++++++ 1 file changed, 222 insertions(+) create mode 100644 docs/overlap_add_math.ipynb diff --git a/docs/overlap_add_math.ipynb b/docs/overlap_add_math.ipynb new file mode 100644 index 0000000..346a29e --- /dev/null +++ b/docs/overlap_add_math.ipynb @@ -0,0 +1,222 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "e984e902-bf76-43ae-a5d5-55e59357b6bf", + "metadata": {}, + "source": [ + "# From Linear Convolution to Overlap-Add Method!" + ] + }, + { + "cell_type": "markdown", + "id": "e9725572-a39b-4fb5-992c-5f3d22a80797", + "metadata": {}, + "source": [ + "This notebook explains the math behind the overlap-add method in convolution. The content 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/)" + ] + }, + { + "cell_type": "markdown", + "id": "03f9d2f8-1560-4a5d-93a2-ec82fb422f6f", + "metadata": {}, + "source": [ + "Most textbooks use $x$ and $h$ notation in the convolution. However, in this notebook, we would like to use $Q$ and $T$ as those are the name of arrays in the realm of STUMPY. We would also like to introduce $Q'$, which is simply the reverse of $Q$." + ] + }, + { + "cell_type": "markdown", + "id": "3dda604a-3a3e-42f4-90c7-7b01ada16337", + "metadata": {}, + "source": [ + "# Linear Convolution of Two Finite-Length Sequences" + ] + }, + { + "cell_type": "markdown", + "id": "9ee474bd-1229-4749-80fe-f90c4a135931", + "metadata": {}, + "source": [ + "Consider two finte-length sequences $Q'$ (of length $m$) and $T$ (of length $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 for any index that is outside of their range. So, $T[i]$ is zero if $i$ is outside of the range $0 \\le i \\le n-1$. Hence, the equation above is equivalent to:\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]$ is not zero when $0 \\le idx-i \\le m-1$, or equivalently $i \\le idx \\le i+m-1$. \n", + "\n", + "Therefore:\n", + "\n", + "* $idx \\ge i, \\quad i \\ge 0 \\implies idx \\ge 0$\n", + "* $idx \\le i+m-1, \\quad i \\le n-1 \\implies idx \\le n+m-2$" + ] + }, + { + "cell_type": "markdown", + "id": "e638d3d7-cf2a-4355-aff3-a9192400defb", + "metadata": {}, + "source": [ + "This shows that the linear convolution $C$ has values at indices $\\set{0, 1, ..., n + m - 2}$, and it is 0 otherwise. Therefore, the length of output in linear convolution is $n+m-1$. \n", + "\n", + "Furthermore, let's check out the values for different indices:\n", + "* $idx=0 \\implies C[0]=T[0]Q'[0]$\n", + "* $idx=1 \\implies C[0]=T[0]Q'[1] + T[1]Q'[0]$\n", + "* ...\n", + "* $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$\n", + "* $idx=m \\implies C[m]=T[1]Q'[m-1] + T[2]Q'[m-2] + ... + T[m]Q'[0] = T_{1}.Q$\n", + "* ...\n", + "* $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+1}.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", + "As observed, when $m-1 \\le idx \\le n-1$, $C[idx]$ becomes the dot product between a subsequnce of $T$ and $Q$, which is the reverse of $Q'$. Therefore, the `range(m-1,n)` gives the sliding dot product between $Q$ and $T$." + ] + }, + { + "cell_type": "markdown", + "id": "48c7f9d8-1b7c-4abc-886b-27a015b78ed6", + "metadata": {}, + "source": [ + "**How can we leverage this to speed up the computation of sliding dot product?**" + ] + }, + { + "cell_type": "markdown", + "id": "52bcf566-3ac0-483d-81a3-57fa3b5766c3", + "metadata": {}, + "source": [ + "# Option I: Convert to Circular Convolution and use FFT-IFFT" + ] + }, + { + "cell_type": "markdown", + "id": "2c4176e2-60f3-410f-ad2d-6c1d29cca6cf", + "metadata": {}, + "source": [ + "Circular convolution is defined between two sequences that are both periodic and their period are the same, say $N$. Their circular convolution is also a periodic sequence, with period $N$, and it can be computed as follows:\n", + "\n", + "$$C_{N}[idx] = \\sum_{i=0}^{N-1} \\tilde{Q'}[i] \\times \\tilde{T}[(idx-i)_{N}], \\quad 0 \\le idx \\le N-1$$\n", + "\n", + "where, $\\tilde{Q'}$ and $\\tilde{T}$ are both periodic sequence with period $N$. $C_{N}$ represents one period of N-Circular convolution. This can also be computed via FFT, i.e. $$C_{N} = IFFT( FFT(\\tilde{Q'}_{N}) \\times FFT(\\tilde{T}_{N}) ) $$" + ] + }, + { + "cell_type": "markdown", + "id": "e10c5304-a1a1-4cbd-9b26-266e229bdc57", + "metadata": {}, + "source": [ + "If there is a way to compute the linear convolution via circular convolution, then we can take advantage of FFT by using eq (4). The good news is that there is a way! The linear convolution between $Q'$ and $T$ can be obtained by performing N-circular convolution between $\\tilde{Q'}$ and $\\tilde{T}$, where:\n", + "\n", + "* $N \\ge n + m - 2$\n", + "* $\\tilde{Q'}_{N}$ is $Q'$ but zero-padded with $N-m$ zeros\n", + "* $\\tilde{T}_{N}$ is $T$ but zero-padded with $N-n$ zeros" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "id": "b6529137-fcdc-4d65-aa93-da3746739f94", + "metadata": {}, + "outputs": [], + "source": [ + "# WIP" + ] + }, + { + "cell_type": "markdown", + "id": "97fcc328-35f9-4b66-9598-1b5d8f1f97ac", + "metadata": {}, + "source": [ + "# Option II: Overlap-add method\n", + "This is a divide and conquer algorithm. It applies divide on \"linear convolution\" and conquer each via circular convolution (See: Option I)" + ] + }, + { + "cell_type": "markdown", + "id": "1e45ccf7-1564-4c41-b32e-4de9b1b201a3", + "metadata": {}, + "source": [ + "## Linearity 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 long $T$ by breaking $T$ into smaller parts." + ] + }, + { + "cell_type": "markdown", + "id": "fecee5fc-3c10-4c2e-8c66-bce17c3650bd", + "metadata": {}, + "source": [] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "acf5516d-e52d-457c-a8c5-f4d1b7ae8f1c", + "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 +} From f10007723333fb1133840f0aa122a315360a635e Mon Sep 17 00:00:00 2001 From: NimaSarajpoor Date: Thu, 3 Sep 2026 01:30:28 -0400 Subject: [PATCH 2/8] complete Option 1 part and added bonus --- docs/overlap_add_math.ipynb | 96 +++++++++++++++++++++++++++++++++---- 1 file changed, 88 insertions(+), 8 deletions(-) diff --git a/docs/overlap_add_math.ipynb b/docs/overlap_add_math.ipynb index 346a29e..8d0b9e7 100644 --- a/docs/overlap_add_math.ipynb +++ b/docs/overlap_add_math.ipynb @@ -115,9 +115,9 @@ "source": [ "Circular convolution is defined between two sequences that are both periodic and their period are the same, say $N$. Their circular convolution is also a periodic sequence, with period $N$, and it can be computed as follows:\n", "\n", - "$$C_{N}[idx] = \\sum_{i=0}^{N-1} \\tilde{Q'}[i] \\times \\tilde{T}[(idx-i)_{N}], \\quad 0 \\le idx \\le N-1$$\n", + "$$C_{N}[idx] = \\sum_{i=0}^{N-1} \\tilde{T}[i] \\times \\tilde{Q'}[(idx-i)_{N}], \\quad 0 \\le idx \\le N-1$$\n", "\n", - "where, $\\tilde{Q'}$ and $\\tilde{T}$ are both periodic sequence with period $N$. $C_{N}$ represents one period of N-Circular convolution. This can also be computed via FFT, i.e. $$C_{N} = IFFT( FFT(\\tilde{Q'}_{N}) \\times FFT(\\tilde{T}_{N}) ) $$" + "where, $\\tilde{Q'}$ and $\\tilde{T}$ are both periodic sequence with period $N$. $C_{N}$ represents one period of N-Circular convolution. This can also be computed via FFT, i.e. $C_{N} = IFFT( FFT(\\tilde{T}_{N}) \\times FFT(\\tilde{Q'}_{N}) )$." ] }, { @@ -125,21 +125,101 @@ "id": "e10c5304-a1a1-4cbd-9b26-266e229bdc57", "metadata": {}, "source": [ - "If there is a way to compute the linear convolution via circular convolution, then we can take advantage of FFT by using eq (4). The good news is that there is a way! The linear convolution between $Q'$ and $T$ can be obtained by performing N-circular convolution between $\\tilde{Q'}$ and $\\tilde{T}$, where:\n", + "If there is a way to compute the linear convolution via circular convolution, then we can use FFT to compute the linear convolution. Actually, there is a way! The linear convolution between $Q'$ and $T$ is equivalent to the N-circular convolution between $\\tilde{Q'}$ and $\\tilde{T}$, where:\n", "\n", - "* $N \\ge n + m - 2$\n", + "* $N \\ge n + m - 1$, to cover the $n+m-1$ data points of linear convolution\n", "* $\\tilde{Q'}_{N}$ is $Q'$ but zero-padded with $N-m$ zeros\n", "* $\\tilde{T}_{N}$ is $T$ but zero-padded with $N-n$ zeros" ] }, { - "cell_type": "code", - "execution_count": 2, + "cell_type": "markdown", "id": "b6529137-fcdc-4d65-aa93-da3746739f94", "metadata": {}, - "outputs": [], "source": [ - "# WIP" + "So:\n", + "\n", + "$$\n", + "\\tilde{Q'}_{N}[i] = \\begin{cases} \n", + " Q'[i] & \\text{if } 0 \\le i \\le m-1 \\\\\n", + " 0 & \\text{if } m \\le i \\le N-1\n", + " \\end{cases}\n", + "$$\n", + "\n", + "$$\n", + "\\tilde{T}_{N}[i] = \\begin{cases} \n", + " T[i] & \\text{if } 0 \\le i \\le n-1 \\\\\n", + " 0 & \\text{if } n \\le i \\le N-1\n", + " \\end{cases}\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "id": "5ce64782-66d0-413b-a116-9ac3555966b6", + "metadata": {}, + "source": [ + "Let's break down the circular convolution equation:\n", + "\n", + "$$C_{N}[idx] = \\sum_{i=0}^{m-1} \\tilde{T}[i] \\times \\tilde{Q'}[(idx-i)_{N}] + \\sum_{i=m}^{n-1} \\tilde{T}[i] \\times \\tilde{Q'}[(idx-i)_{N}] + \\sum_{i=n}^{N-1} \\tilde{T}[i] \\times \\tilde{Q'}[(idx-i)_{N}]$$\n", + "\n", + "$$C_{N}[idx] = \\sum_{i=0}^{m-1} T[i] \\times \\tilde{Q'}[(idx-i)_{N}] + \\sum_{i=m}^{n-1} T[i] \\times \\tilde{Q'}[(idx-i)_{N}]$$\n", + "\n", + "$$C_{N}[idx] = \\sum_{i=0}^{n-1} T[i] \\times \\tilde{Q'}[(idx-i)_{N}]$$" + ] + }, + { + "cell_type": "markdown", + "id": "a7e66371-4cd4-4673-957c-b7c548dcbfde", + "metadata": {}, + "source": [ + "Let's compare the linear convolution in eq(2) and the N-circular convolution in eq(9). To show that they give the same value, we need to compare $Q'[idx-i]$ with $\\tilde{Q'}[(idx-i)_{N}]$.\n", + "\n", + "Let $j$ denote $idx-i$. So, all we need to do is to compare $Q'[j]$ with $\\tilde{Q'}[(j)_{N}]$.\n", + "\n", + "**Case A: $j \\ge 0$**
\n", + "\n", + "In this case, $max(j) = max(idx) - min(i) = (N-1) - 0 = N-1 \\implies \\tilde{Q'}[(j)_{N}]=Q'_{N}[j]$\n", + "\n", + "\n", + "**Case B: $j < 0$**
\n", + "\n", + "In this case, $Q'[j]=0$.\n", + "\n", + "In this case, $\\tilde{Q'}[(j)_{N}] = Q'_{N}[j+N]$. Note that: \n", + "$$j+N \\ge min(j)+N \\ge (min(idx)-max(i)) + N \\ge (0 - (n-1)) + N$$ \n", + "\n", + "Since $Q'_{N}$ is padded with zeros for indices that are at least $m$, therefore the value $\\tilde{Q'}[(j)_{N}] = Q'_{N}[j+N]$ can be 0 if the condition $-(n-1)+N \\ge m$ is satisfied. The condition can be written as $N \\ge n + m - 1$, which is satisfied!" + ] + }, + { + "cell_type": "markdown", + "id": "181cc549-9b2a-48ac-9077-71174cb49f1c", + "metadata": {}, + "source": [ + "# Option I: Bonus\n", + "\n", + "Let's revisit the linear convolution, shown in eq (2). Recall that the length of output is $n + m - 1$. The sliding dot product between $Q$ and $T$ is in $range(m-1, n)$ though. Therefore, if the goal is to calculate the sliding dot product, the computation can stop at index $n$. In other words, $idx$ can change from 0 to $n$ (exclusive). How about the N-Circular convolution in this case? If we set $N$ to $n$ instead of $n+m-1$, does the outcome have the correct answer for the slice in $range(m-1,n)$? We need to show that eq(2) and eq(8) give the same result for $m-1 \\le idx \\le n-1$ when $N=n$.\n", + "\n", + "**Case A: $j \\ge 0$**
\n", + "\n", + "Same as before!\n", + "\n", + "**Case B: $j < 0$**
\n", + "In this case, $Q'[j]=0$.\n", + "\n", + "In this case, $\\tilde{Q'}[(j)_{N}] = Q'_{N}[j+N]$. Note that: \n", + "$$j+N \\ge min(j)+N \\ge (min(idx)-max(i)) + N \\ge ((m-1) - (n-1)) + N$$ \n", + "\n", + "Since $Q'_{N}$ is padded with zeros for indices that are at least $m$, therefore the value $\\tilde{Q'}[(j)_{N}] = Q'_{N}[j+N]$ can be 0 if the condition $m-n + N \\ge m$ is satisfied. The condition can be written as $N \\ge n$, which is satisfied!" + ] + }, + { + "cell_type": "markdown", + "id": "03cbd90d-a7e4-466d-af15-81d9a1c4ce6c", + "metadata": {}, + "source": [ + "This shows that the sliding dot product of $Q$ and $T$ can be computed via $n$-Circular convolution rather than $(n+m-1)$-circular convolution" ] }, { From 520990368fc273da1ca82b8f4afdc1a50abd6e1b Mon Sep 17 00:00:00 2001 From: NimaSarajpoor Date: Thu, 3 Sep 2026 22:32:28 -0400 Subject: [PATCH 3/8] minor changes --- docs/overlap_add_math.ipynb | 62 ++++++++++++++++++++----------------- 1 file changed, 33 insertions(+), 29 deletions(-) diff --git a/docs/overlap_add_math.ipynb b/docs/overlap_add_math.ipynb index 8d0b9e7..fa43df6 100644 --- a/docs/overlap_add_math.ipynb +++ b/docs/overlap_add_math.ipynb @@ -5,7 +5,7 @@ "id": "e984e902-bf76-43ae-a5d5-55e59357b6bf", "metadata": {}, "source": [ - "# From Linear Convolution to Overlap-Add Method!" + "# From Linear Convolution to Overlap-Add Method" ] }, { @@ -13,7 +13,7 @@ "id": "e9725572-a39b-4fb5-992c-5f3d22a80797", "metadata": {}, "source": [ - "This notebook explains the math behind the overlap-add method in convolution. The content 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/)" + "This notebook explains the math behind linear convolution, circular convolution, and the overlap-add method in convolution. The content 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/)" ] }, { @@ -21,7 +21,7 @@ "id": "03f9d2f8-1560-4a5d-93a2-ec82fb422f6f", "metadata": {}, "source": [ - "Most textbooks use $x$ and $h$ notation in the convolution. However, in this notebook, we would like to use $Q$ and $T$ as those are the name of arrays in the realm of STUMPY. We would also like to introduce $Q'$, which is simply the reverse of $Q$." + "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 name of arrays in the realm of STUMPY. Also, we would like to introduce $Q'$, which is simply the reverse of $Q$." ] }, { @@ -37,7 +37,7 @@ "id": "9ee474bd-1229-4749-80fe-f90c4a135931", "metadata": {}, "source": [ - "Consider two finte-length sequences $Q'$ (of length $m$) and $T$ (of length $n$). Their linear convolution, $C$, can be computed as follows:" + "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:" ] }, { @@ -45,7 +45,7 @@ "id": "ad4a63b5-2d91-4875-b4c2-bc7ce11b6af9", "metadata": {}, "source": [ - "$$ C[idx] = \\sum_{i=-\\infty}^{\\infty}{T[i] \\times Q'[idx-i]} $$" + "$$ C[idx] = \\sum_{i=-\\infty}^{\\infty}{T[i] \\times Q'[idx-i]}, \\quad 0 \\le idx < n + m - 1$$" ] }, { @@ -55,7 +55,7 @@ "source": [ "where $T[.]$ and $Q'[.]$ are zeros for any index that is outside of their range. So, $T[i]$ is zero if $i$ is outside of the range $0 \\le i \\le n-1$. Hence, the equation above is equivalent to:\n", "\n", - "$$ C[idx] = \\sum_{i=0}^{n-1}{T[i] \\times Q'[idx-i]} $$" + "$$ C[idx] = \\sum_{i=0}^{n-1}{T[i] \\times Q'[idx-i]}, \\quad 0 \\le idx < n + m - 1$$" ] }, { @@ -63,7 +63,7 @@ "id": "593a8276-49cb-43d2-9300-8d29e67d3104", "metadata": {}, "source": [ - "Note that $Q'[idx-i]$ is not zero when $0 \\le idx-i \\le m-1$, or equivalently $i \\le idx \\le i+m-1$. \n", + "Note that $Q'[idx-i]$ can be a non-zero value when $0 \\le idx-i \\le m-1$, or equivalently $i \\le idx \\le i+m-1$. \n", "\n", "Therefore:\n", "\n", @@ -76,9 +76,7 @@ "id": "e638d3d7-cf2a-4355-aff3-a9192400defb", "metadata": {}, "source": [ - "This shows that the linear convolution $C$ has values at indices $\\set{0, 1, ..., n + m - 2}$, and it is 0 otherwise. Therefore, the length of output in linear convolution is $n+m-1$. \n", - "\n", - "Furthermore, let's check out the values for different indices:\n", + "This shows that the linear convolution $C$ has values at indices $\\set{0, 1, ..., n + m - 2}$, and it is 0 otherwise. Furthermore, let's check out the values at different indices:\n", "* $idx=0 \\implies C[0]=T[0]Q'[0]$\n", "* $idx=1 \\implies C[0]=T[0]Q'[1] + T[1]Q'[0]$\n", "* ...\n", @@ -89,15 +87,9 @@ "* ...\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", - "As observed, when $m-1 \\le idx \\le n-1$, $C[idx]$ becomes the dot product between a subsequnce of $T$ and $Q$, which is the reverse of $Q'$. Therefore, the `range(m-1,n)` gives the sliding dot product between $Q$ and $T$." - ] - }, - { - "cell_type": "markdown", - "id": "48c7f9d8-1b7c-4abc-886b-27a015b78ed6", - "metadata": {}, - "source": [ - "**How can we leverage this to speed up the computation of sliding dot product?**" + "As observed, when $m-1 \\le idx \\le n-1$, $C[idx]$ becomes the dot product between $Q$ (the reverse of $Q'$) and a subsequnce of $T$. Therefore, $range(m-1,n)$ gives the sliding dot product between $Q$ and $T$. \n", + "\n", + "The time complexity of calcularing the linear convolution is $O(nm)$. Can we do better? Yes!" ] }, { @@ -159,11 +151,11 @@ "id": "5ce64782-66d0-413b-a116-9ac3555966b6", "metadata": {}, "source": [ - "Let's break down the circular convolution equation:\n", + "Let's prove this! We start with breaking down the circular convolution equation:\n", "\n", - "$$C_{N}[idx] = \\sum_{i=0}^{m-1} \\tilde{T}[i] \\times \\tilde{Q'}[(idx-i)_{N}] + \\sum_{i=m}^{n-1} \\tilde{T}[i] \\times \\tilde{Q'}[(idx-i)_{N}] + \\sum_{i=n}^{N-1} \\tilde{T}[i] \\times \\tilde{Q'}[(idx-i)_{N}]$$\n", + "$$C_{N}[idx] = \\sum_{i=0}^{n-1} \\tilde{T}[i] \\times \\tilde{Q'}[(idx-i)_{N}] + \\sum_{i=n}^{N-1} \\tilde{T}[i] \\times \\tilde{Q'}[(idx-i)_{N}]$$\n", "\n", - "$$C_{N}[idx] = \\sum_{i=0}^{m-1} T[i] \\times \\tilde{Q'}[(idx-i)_{N}] + \\sum_{i=m}^{n-1} T[i] \\times \\tilde{Q'}[(idx-i)_{N}]$$\n", + "In the second term, $i \\ge n$. Therefore, $\\tilde{T}[i]$ is zero. Hence, the equation above becomes:\n", "\n", "$$C_{N}[idx] = \\sum_{i=0}^{n-1} T[i] \\times \\tilde{Q'}[(idx-i)_{N}]$$" ] @@ -173,9 +165,9 @@ "id": "a7e66371-4cd4-4673-957c-b7c548dcbfde", "metadata": {}, "source": [ - "Let's compare the linear convolution in eq(2) and the N-circular convolution in eq(9). To show that they give the same value, we need to compare $Q'[idx-i]$ with $\\tilde{Q'}[(idx-i)_{N}]$.\n", + "Let's compare the linear convolution in eq(2) and the N-circular convolution in eq(7). To show that they give the same value, we need to compare $Q'[idx-i]$ with $\\tilde{Q'}[(idx-i)_{N}]$.\n", "\n", - "Let $j$ denote $idx-i$. So, all we need to do is to compare $Q'[j]$ with $\\tilde{Q'}[(j)_{N}]$.\n", + "Let's define: $j = idx-i$. So, all we need to do is to compare $Q'[j]$ with $\\tilde{Q'}[(j)_{N}]$.\n", "\n", "**Case A: $j \\ge 0$**
\n", "\n", @@ -197,9 +189,9 @@ "id": "181cc549-9b2a-48ac-9077-71174cb49f1c", "metadata": {}, "source": [ - "# Option I: Bonus\n", + "# Option I: Enhanced\n", "\n", - "Let's revisit the linear convolution, shown in eq (2). Recall that the length of output is $n + m - 1$. The sliding dot product between $Q$ and $T$ is in $range(m-1, n)$ though. Therefore, if the goal is to calculate the sliding dot product, the computation can stop at index $n$. In other words, $idx$ can change from 0 to $n$ (exclusive). How about the N-Circular convolution in this case? If we set $N$ to $n$ instead of $n+m-1$, does the outcome have the correct answer for the slice in $range(m-1,n)$? We need to show that eq(2) and eq(8) give the same result for $m-1 \\le idx \\le n-1$ when $N=n$.\n", + "Let's revisit the linear convolution, shown in eq (2). Recall that the length of output is $n + m - 1$. The sliding dot product between $Q$ and $T$ is in $range(m-1, n)$ though. Therefore, if the goal is to calculate the sliding dot product, the computation can stop at index $n$. In other words, $idx$ can change from 0 to $n$ (exclusive). How about the N-Circular convolution in this case? If we set $N$ to $n$ instead of $n+m-1$, does the outcome have the correct answer for the slice in $range(m-1,n)$? We need to show that eq(2) and eq(7) give the same result for $m-1 \\le idx \\le n-1$ when $N=n$.\n", "\n", "**Case A: $j \\ge 0$**
\n", "\n", @@ -209,9 +201,13 @@ "In this case, $Q'[j]=0$.\n", "\n", "In this case, $\\tilde{Q'}[(j)_{N}] = Q'_{N}[j+N]$. Note that: \n", - "$$j+N \\ge min(j)+N \\ge (min(idx)-max(i)) + N \\ge ((m-1) - (n-1)) + N$$ \n", + "$$j+N \\ge min(j)+N \\ge (min(idx)-max(i)) + N \\ge ((m-1) - (n-1)) + N \\ge m - n + N$$ \n", "\n", - "Since $Q'_{N}$ is padded with zeros for indices that are at least $m$, therefore the value $\\tilde{Q'}[(j)_{N}] = Q'_{N}[j+N]$ can be 0 if the condition $m-n + N \\ge m$ is satisfied. The condition can be written as $N \\ge n$, which is satisfied!" + "Recall that $Q'_{N}$ is padded with zeros for indices that are at least $m$, therefore the value $\\tilde{Q'}[(j)_{N}] = Q'_{N}[j+N]$ can be 0 as long as $j+N \\ge m$. the condition $m-n + N \\ge m$ is satisfied. Earlier, we showed that $j + N \\ge m - n + N$. Therefore, we just need to make sure $m - n + N$ is at least $m$. \n", + "\n", + "$$m - n + N \\ge m \\implies N \\ge n$$ \n", + "\n", + "And $N \\ge n$ is satisfied!" ] }, { @@ -219,7 +215,15 @@ "id": "03cbd90d-a7e4-466d-af15-81d9a1c4ce6c", "metadata": {}, "source": [ - "This shows that the sliding dot product of $Q$ and $T$ can be computed via $n$-Circular convolution rather than $(n+m-1)$-circular convolution" + "We just showed that the sliding dot product of $Q$ and $T$ can be computed via $n$-Circular convolution. Computing this in the frequency domain will have the time complexity of $O(nlogn)$." + ] + }, + { + "cell_type": "markdown", + "id": "d21fceee-2bab-4d01-9bee-83f73fe53e08", + "metadata": {}, + "source": [ + "Can we do better? See the next approach!" ] }, { From 24ab0d0a76a0e8e483d06a9188942a0de2744dfa Mon Sep 17 00:00:00 2001 From: NimaSarajpoor Date: Mon, 7 Sep 2026 23:16:19 -0400 Subject: [PATCH 4/8] revise objective section and structure --- docs/overlap_add_math.ipynb | 307 ++++++++++++++++++++++++++---------- 1 file changed, 220 insertions(+), 87 deletions(-) diff --git a/docs/overlap_add_math.ipynb b/docs/overlap_add_math.ipynb index fa43df6..ecc5cf0 100644 --- a/docs/overlap_add_math.ipynb +++ b/docs/overlap_add_math.ipynb @@ -2,18 +2,23 @@ "cells": [ { "cell_type": "markdown", - "id": "e984e902-bf76-43ae-a5d5-55e59357b6bf", + "id": "5c9355b9-0bd4-4740-9721-7394a19f4588", "metadata": {}, "source": [ - "# From Linear Convolution to Overlap-Add Method" + "# Objective" ] }, { "cell_type": "markdown", - "id": "e9725572-a39b-4fb5-992c-5f3d22a80797", + "id": "199eb655-fde9-4ee1-82e4-55bad5892f22", "metadata": {}, "source": [ - "This notebook explains the math behind linear convolution, circular convolution, and the overlap-add method in convolution. The content 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/)" + "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) provided different versions of MASS for computing the sliding dot product. At their core, they all come down to computing one/more convolution(s) via the frequency domain. To understand how convolution comes into the picture and how it is used for computing the sliding dot product, this notebook is created to answer the following questions:\n", + "\n", + "1. **How is the sliding dot product (sdp) related to (linear) convolution?** This can help us discover the relationship between sdp and linear convolution. Convolution is a well-studied area and understanding this relationship allows us to leverage convolution methods to compute the sliding dot product.\n", + "2. **How can we compute the convolution faster?** Once we understand the relationship between sdp and convolution, we can focus on learning/using method for faster convolution.\n", + "3. **Can we reduce the output size, and hence the computational load without losing the sdp?** Once we know how to get faster convolution, we can try to enhance its performance further by avoiding unnecessary calculation.\n", + "4. **What to do if $T$ is very long? Overlap-add method!** This is where we learn a divide-and-conquer algorithm on convolution, a useful method when input $T$ is large" ] }, { @@ -21,7 +26,9 @@ "id": "03f9d2f8-1560-4a5d-93a2-ec82fb422f6f", "metadata": {}, "source": [ - "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 name of arrays in the realm of STUMPY. Also, we would like to introduce $Q'$, which is simply the reverse of $Q$." + "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 name of arrays in the realm of 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 defined by\n", + "\n", + "$$Q'[j] = Q[m-1-j], \\quad 0 \\le j \\le m-1,$$" ] }, { @@ -29,7 +36,15 @@ "id": "3dda604a-3a3e-42f4-90c7-7b01ada16337", "metadata": {}, "source": [ - "# Linear Convolution of Two Finite-Length Sequences" + "# 1. How is the sliding dot product(sdp) related to (linear) convolution?" + ] + }, + { + "cell_type": "markdown", + "id": "877af78c-0026-44c0-8295-0f75db075d97", + "metadata": {}, + "source": [ + "### 1.1 Linear Convolution of Two Finite-Length Sequences $Q'$ and $T$" ] }, { @@ -37,7 +52,7 @@ "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:" + "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:" ] }, { @@ -45,7 +60,7 @@ "id": "ad4a63b5-2d91-4875-b4c2-bc7ce11b6af9", "metadata": {}, "source": [ - "$$ C[idx] = \\sum_{i=-\\infty}^{\\infty}{T[i] \\times Q'[idx-i]}, \\quad 0 \\le idx < n + m - 1$$" + "$$ C[idx] = \\sum_{i=-\\infty}^{\\infty}{T[i] \\times Q'[idx-i]}$$" ] }, { @@ -53,9 +68,11 @@ "id": "8028ee30-3244-4c79-86ca-22e2b8876694", "metadata": {}, "source": [ - "where $T[.]$ and $Q'[.]$ are zeros for any index that is outside of their range. So, $T[i]$ is zero if $i$ is outside of the range $0 \\le i \\le n-1$. Hence, the equation above is equivalent to:\n", + "where $T[.]$ and $Q'[.]$ are zeros for any index that is outside of their range. The linear convolution can be seen as a flip and slide operation, and at each index $idx$, the output 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. For instance, at index $idx=0$, there is only one overlapping element between $T$ and $Q'$, and that $T[0]$ and $Q'[0]$. And, $C[0] = T[0]Q'[0]$. Therefore, **the linear convolution is NOT the same as the sliding dot product** since sliding dot product is for cases when the overlap covers the full query and not just one element. \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 revised as follows:\n", "\n", - "$$ C[idx] = \\sum_{i=0}^{n-1}{T[i] \\times Q'[idx-i]}, \\quad 0 \\le idx < n + m - 1$$" + "$$ C[idx] = \\sum_{i=0}^{n-1}{T[i] \\times Q'[idx-i]}$$" ] }, { @@ -63,12 +80,23 @@ "id": "593a8276-49cb-43d2-9300-8d29e67d3104", "metadata": {}, "source": [ - "Note that $Q'[idx-i]$ can be a non-zero value when $0 \\le idx-i \\le m-1$, or equivalently $i \\le idx \\le i+m-1$. \n", + "Note that $Q'[idx-i]$ can be a non-zero value 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", - "Therefore:\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", - "* $idx \\ge i, \\quad i \\ge 0 \\implies idx \\ge 0$\n", - "* $idx \\le i+m-1, \\quad i \\le n-1 \\implies idx \\le n+m-2$" + "So, when $idx$ changes from $0$ to $n+m-2$, there exist an $i$ such that both $T[i]$ and $Q'[idx-i]$ can be non-zero values. This shows that the linear convolution $C$ can have values at indices $\\set{0, 1, ..., n + m - 2}$. " + ] + }, + { + "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$)" ] }, { @@ -76,20 +104,33 @@ "id": "e638d3d7-cf2a-4355-aff3-a9192400defb", "metadata": {}, "source": [ - "This shows that the linear convolution $C$ has values at indices $\\set{0, 1, ..., n + m - 2}$, and it is 0 otherwise. Furthermore, let's check out the values at different indices:\n", + "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", + "> **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[0]=T[0]Q'[1] + T[1]Q'[0]$\n", "* ...\n", - "* $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$\n", - "* $idx=m \\implies C[m]=T[1]Q'[m-1] + T[2]Q'[m-2] + ... + T[m]Q'[0] = T_{1}.Q$\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}$\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}$\n", "* ...\n", - "* $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+1}.Q$\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", - "As observed, when $m-1 \\le idx \\le n-1$, $C[idx]$ becomes the dot product between $Q$ (the reverse of $Q'$) and a subsequnce of $T$. Therefore, $range(m-1,n)$ gives the sliding dot product between $Q$ and $T$. \n", + "When $m-1 \\le idx \\le n-1$, $C[idx]$ actually becomes the dot product between $Q$ (the reverse of $Q'$) and $T_{idx-(m-1)}$. 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", - "The time complexity of calcularing the linear convolution is $O(nm)$. Can we do better? Yes!" + "So, once we calculate the linear convolution between $Q'$ and $T$, we have the sliding dot product between $Q$ and $T$ for free! The time complexity of calculating the linear convolution is $O(nm)$. Can we do better? " + ] + }, + { + "cell_type": "markdown", + "id": "8475be6a-9515-43ae-a8d7-083b8db3af51", + "metadata": {}, + "source": [ + "# 2. How can we compute the convolution faster than O(nm)?" ] }, { @@ -97,7 +138,7 @@ "id": "52bcf566-3ac0-483d-81a3-57fa3b5766c3", "metadata": {}, "source": [ - "# Option I: Convert to Circular Convolution and use FFT-IFFT" + "### 2.1 Convert linear convolution to circular convolution (Step I)" ] }, { @@ -105,11 +146,11 @@ "id": "2c4176e2-60f3-410f-ad2d-6c1d29cca6cf", "metadata": {}, "source": [ - "Circular convolution is defined between two sequences that are both periodic and their period are the same, say $N$. Their circular convolution is also a periodic sequence, with period $N$, and it can be computed as follows:\n", + "Circular convolution is defined between two sequences that are both periodic and their periods are the same, say $N$. Their circular convolution is also a periodic sequence, with period $N$, and it can be computed as follows:\n", "\n", "$$C_{N}[idx] = \\sum_{i=0}^{N-1} \\tilde{T}[i] \\times \\tilde{Q'}[(idx-i)_{N}], \\quad 0 \\le idx \\le N-1$$\n", "\n", - "where, $\\tilde{Q'}$ and $\\tilde{T}$ are both periodic sequence with period $N$. $C_{N}$ represents one period of N-Circular convolution. This can also be computed via FFT, i.e. $C_{N} = IFFT( FFT(\\tilde{T}_{N}) \\times FFT(\\tilde{Q'}_{N}) )$." + "where, $C_{N}$ represents one period of convolution output. $\\tilde{Q'}$ and $\\tilde{T}$ are both periodic sequence with period $N$, and the index $(idx-i)_{N}$ simply means \"$idx-i$ modulo $N$\", i.e. $(idx-i) \\% N$, and it is always between $0$ and $N-1$." ] }, { @@ -117,45 +158,50 @@ "id": "e10c5304-a1a1-4cbd-9b26-266e229bdc57", "metadata": {}, "source": [ - "If there is a way to compute the linear convolution via circular convolution, then we can use FFT to compute the linear convolution. Actually, there is a way! The linear convolution between $Q'$ and $T$ is equivalent to the N-circular convolution between $\\tilde{Q'}$ and $\\tilde{T}$, where:\n", + "The time complexity of calculating $C_{N}$ using the equation above is $O(N^{2})$. However, as will be shown later in the next section, the time complexity can be reduced to $O(NlogN)$. However, this means nothing for sliding dot product unless we understand how the linear convolution can be computed via circular convolution and how $N$ is related to $n$ and/or $m$. \n", + "\n", + "We first prove that the linear convolution between $Q'$ and $T$ is equivalent to $C_{N}$ between $\\tilde{Q'}$ and $\\tilde{T}$ if:\n", "\n", - "* $N \\ge n + m - 1$, to cover the $n+m-1$ data points of linear convolution\n", - "* $\\tilde{Q'}_{N}$ is $Q'$ but zero-padded with $N-m$ zeros\n", - "* $\\tilde{T}_{N}$ is $T$ but zero-padded with $N-n$ zeros" + "* $N = n + m - 1$\n", + "* $\\tilde{Q'}_{N}$, i.e. one period of $\\tilde{Q'}$, is $Q'$ but with $N-m$ zero padding\n", + "* $\\tilde{T}_{N}$, i.e. one period of $\\tilde{T}$, is $T$ but with $N-n$ zero padding" ] }, { "cell_type": "markdown", - "id": "b6529137-fcdc-4d65-aa93-da3746739f94", + "id": "5ce64782-66d0-413b-a116-9ac3555966b6", "metadata": {}, "source": [ - "So:\n", + "**Proof**\n", "\n", - "$$\n", + "> **Note:**
\n", + "> Keep in mind that $N = n + m - 1$, and:\n", + "> \n", + "> $$\n", "\\tilde{Q'}_{N}[i] = \\begin{cases} \n", " Q'[i] & \\text{if } 0 \\le i \\le m-1 \\\\\n", " 0 & \\text{if } m \\le i \\le N-1\n", " \\end{cases}\n", "$$\n", - "\n", - "$$\n", + "> \n", + "> $$\n", "\\tilde{T}_{N}[i] = \\begin{cases} \n", " T[i] & \\text{if } 0 \\le i \\le n-1 \\\\\n", " 0 & \\text{if } n \\le i \\le N-1\n", " \\end{cases}\n", - "$$" - ] - }, - { - "cell_type": "markdown", - "id": "5ce64782-66d0-413b-a116-9ac3555966b6", - "metadata": {}, - "source": [ - "Let's prove this! We start with breaking down the circular convolution equation:\n", + "$$\n", + "\n", + "\n", + "\n", + "We start with rewriting the circular convolution equation by breaking the summation into two parts: \n", "\n", "$$C_{N}[idx] = \\sum_{i=0}^{n-1} \\tilde{T}[i] \\times \\tilde{Q'}[(idx-i)_{N}] + \\sum_{i=n}^{N-1} \\tilde{T}[i] \\times \\tilde{Q'}[(idx-i)_{N}]$$\n", "\n", - "In the second term, $i \\ge n$. Therefore, $\\tilde{T}[i]$ is zero. Hence, the equation above becomes:\n", + "Note that:\n", + "* In the first term, $0 \\le i \\le n-1$, and hence $\\tilde{T}[i]$ is the same as $T[i]$.\n", + "* In the second term, $n \\le i \\le N-1$, and hence $\\tilde{T}[i]$ is zero.\n", + "\n", + "Hence, the equation becomes:\n", "\n", "$$C_{N}[idx] = \\sum_{i=0}^{n-1} T[i] \\times \\tilde{Q'}[(idx-i)_{N}]$$" ] @@ -165,65 +211,155 @@ "id": "a7e66371-4cd4-4673-957c-b7c548dcbfde", "metadata": {}, "source": [ - "Let's compare the linear convolution in eq(2) and the N-circular convolution in eq(7). To show that they give the same value, we need to compare $Q'[idx-i]$ with $\\tilde{Q'}[(idx-i)_{N}]$.\n", + "Recall that the linear convolution equation is:\n", + "\n", + "$$ C[idx] = \\sum_{i=0}^{n-1}{T[i] \\times Q'[idx-i]}$$\n", + "\n", + "So, to show that $C_{N}$, from circular convolution, gives the same value as the linear convolution, we only need to show that $\\tilde{Q'}[(idx-i)_{N}]$ and $Q'[idx-i]$ give the same value when $N=n+m-1$. Let $j$ denote $idx-i$. So, all we need to do is to prove that $\\tilde{Q'}[(j)_{N}]$ and $Q'[j]$ 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", - "Let's define: $j = idx-i$. So, all we need to do is to compare $Q'[j]$ with $\\tilde{Q'}[(j)_{N}]$.\n", + "$$max(j) = max(idx-i) = max(idx) + max(-i) = max(idx) - min(i) = (N-1) - 0 = N-1$$ \n", "\n", - "**Case A: $j \\ge 0$**
\n", + "So, $j$ changes from 0 to $N-1$. hence, $j$ modulo $N$ is still $j$. Therefore: $\\tilde{Q'}[(j)_{N}]=Q'[j]$. Proof is now complete for this case.\n", "\n", - "In this case, $max(j) = max(idx) - min(i) = (N-1) - 0 = N-1 \\implies \\tilde{Q'}[(j)_{N}]=Q'_{N}[j]$\n", + "**Case II: $j < 0$**
\n", "\n", + "In this case, $Q'[j]=0$. So, all we need to do is to prove that $\\tilde{Q'}[(j)_{N}]$ is zero as well. \n", "\n", - "**Case B: $j < 0$**
\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", - "In this case, $Q'[j]=0$.\n", + "Therefore,\n", "\n", - "In this case, $\\tilde{Q'}[(j)_{N}] = Q'_{N}[j+N]$. Note that: \n", - "$$j+N \\ge min(j)+N \\ge (min(idx)-max(i)) + N \\ge (0 - (n-1)) + N$$ \n", + "$$min(j) = min(idx-i) = min(idx) + min(-i) = min(idx) - max(i) = 0 - (n-1)$$\n", "\n", - "Since $Q'_{N}$ is padded with zeros for indices that are at least $m$, therefore the value $\\tilde{Q'}[(j)_{N}] = Q'_{N}[j+N]$ can be 0 if the condition $-(n-1)+N \\ge m$ is satisfied. The condition can be written as $N \\ge n + m - 1$, which is satisfied!" + "So, in this case, $-(n-1) \\le j \\le -1$. \n", + "\n", + "Since $\\tilde{Q'}[(j)_{N}] == \\tilde{Q'}[(j+N)_{N}]$, we can compute the latter instead of the former. To do that, we first need to figure out the range for $j+N$.\n", + "\n", + "$$ -(n-1) \\le j \\le -1 \\implies -(n-1)+N \\le j+N \\le -1+N $$\n", + "\n", + "Let's plug the value of $N=n+m-1$ into $-(n-1)+N$:\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'}[(j)_{N}] == \\tilde{Q'}[(j+N)_{N}] = \\tilde{Q'}[j+N]$. Since the index $j+N$ is $\\ge m$, the value $\\tilde{Q'}[j+N]$ becomes 0 as the element resides in the zero-padding part. Proof is now complete for this case." ] }, { "cell_type": "markdown", - "id": "181cc549-9b2a-48ac-9077-71174cb49f1c", + "id": "66a707c1-8004-4fe6-b954-ffe039eda441", "metadata": {}, "source": [ - "# Option I: Enhanced\n", - "\n", - "Let's revisit the linear convolution, shown in eq (2). Recall that the length of output is $n + m - 1$. The sliding dot product between $Q$ and $T$ is in $range(m-1, n)$ though. Therefore, if the goal is to calculate the sliding dot product, the computation can stop at index $n$. In other words, $idx$ can change from 0 to $n$ (exclusive). How about the N-Circular convolution in this case? If we set $N$ to $n$ instead of $n+m-1$, does the outcome have the correct answer for the slice in $range(m-1,n)$? We need to show that eq(2) and eq(7) give the same result for $m-1 \\le idx \\le n-1$ when $N=n$.\n", - "\n", - "**Case A: $j \\ge 0$**
\n", - "\n", - "Same as before!\n", - "\n", - "**Case B: $j < 0$**
\n", - "In this case, $Q'[j]=0$.\n", + "We just proved that the linear convolution can be computed in the form of circular convolution when the arrays have certain lengths and are zero-padded properly." + ] + }, + { + "cell_type": "markdown", + "id": "8648ff47-4d67-4bec-a91e-4ec6d799cc0e", + "metadata": {}, + "source": [ + "### 2.2 Compute Circular Convolution in the Frequency Domain (Step II)" + ] + }, + { + "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 $n$ and $m$ are the length of $T$ and $Q'$, respectively. However, the computation via the 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. So, instead of computing the circular convolution in time domain, i.e.\n", "\n", - "In this case, $\\tilde{Q'}[(j)_{N}] = Q'_{N}[j+N]$. Note that: \n", - "$$j+N \\ge min(j)+N \\ge (min(idx)-max(i)) + N \\ge ((m-1) - (n-1)) + N \\ge m - n + N$$ \n", + "$$C_{N}[idx] = \\sum_{i=0}^{N-1} \\tilde{T}[i] \\times \\tilde{Q'}[(idx-i)_{N}], \\quad 0 \\le idx \\le N-1$$\n", "\n", - "Recall that $Q'_{N}$ is padded with zeros for indices that are at least $m$, therefore the value $\\tilde{Q'}[(j)_{N}] = Q'_{N}[j+N]$ can be 0 as long as $j+N \\ge m$. the condition $m-n + N \\ge m$ is satisfied. Earlier, we showed that $j + N \\ge m - n + N$. Therefore, we just need to make sure $m - n + N$ is at least $m$. \n", + "we can compute it by taking it to the frequency domain:\n", "\n", - "$$m - n + N \\ge m \\implies N \\ge n$$ \n", + "$$C_{N} = IFFT\\left(\n", + "FFT(\\tilde{T}_{N})\n", + "\\cdot\n", + "FFT(\\tilde{Q'}_{N})\n", + "\\right)$$\n", "\n", - "And $N \\ge n$ is satisfied!" + "and its time complexity becomes $O(NlogN)$." ] }, { "cell_type": "markdown", - "id": "03cbd90d-a7e4-466d-af15-81d9a1c4ce6c", + "id": "64b8a7ee-f1dd-4494-a9a4-4ba3293a3912", "metadata": {}, "source": [ - "We just showed that the sliding dot product of $Q$ and $T$ can be computed via $n$-Circular convolution. Computing this in the frequency domain will have the time complexity of $O(nlogn)$." + "So, if I can 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 than before!" ] }, { "cell_type": "markdown", - "id": "d21fceee-2bab-4d01-9bee-83f73fe53e08", + "id": "181cc549-9b2a-48ac-9077-71174cb49f1c", "metadata": {}, "source": [ - "Can we do better? See the next approach!" + "# 3. Can we reduce the output size, and hence the computational load without losing the sdp?\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", + "Recall that the length of output is $n + m - 1$. The sliding dot product between $Q$ and $T$ is in $range(m-1, n)$ though. 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 we accordingly set $N$ to $n$ instead of $n+m-1$, does the slice in $range(m-1,n)$ still reflect the values of sdp? In other words, we need to show that\n", + "\n", + "$$ C[idx] = \\sum_{i=0}^{n-1}{T[i] \\times Q'[idx-i]}, \\quad m-1 \\le idx \\le n-1$$\n", + "\n", + "and\n", + "\n", + "$$C_{N=n}[idx] = \\sum_{i=0}^{n-1}T[i] \\times \\tilde{Q'}[(idx-i)_{n}], \\quad m-1 \\le idx \\le n-1$$\n", + "\n", + "give the same result. To prove this, we just need to show that $Q'[idx-i]$ and $\\tilde{Q'}[(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", + "$$ j \\le 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'}[(j)_{n}]$ is the same as $Q'[j]$. 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'}[(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", + "* $\\tilde{Q'}[(j)_{n}] = \\tilde{Q'}[(j+n)_{n}]$. Note that: \n", + "$$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", + "\n", + "So, when $j<0$, then: $m \\le j+n \\le n-1$. Therefore $\\tilde{Q'}[(j+n)_{n}]$ simply becomes $\\tilde{Q'}[j+n]$, and the value is zero as the index falls into the padded zeros. Therefore, $Q'[j]=\\tilde{Q'}[j+n]=0$. Proof is now complete for this case. " + ] + }, + { + "cell_type": "markdown", + "id": "03cbd90d-a7e4-466d-af15-81d9a1c4ce6c", + "metadata": {}, + "source": [ + "We just showed that the sliding dot product can be computed via circular convolution with period $N=n$. As noted in previous section, this can be computed in $O(nlogn)$ which is faster than our baseline $O(nm)$ unless $m$ is small." ] }, { @@ -231,8 +367,9 @@ "id": "97fcc328-35f9-4b66-9598-1b5d8f1f97ac", "metadata": {}, "source": [ - "# Option II: Overlap-add method\n", - "This is a divide and conquer algorithm. It applies divide on \"linear convolution\" and conquer each via circular convolution (See: Option I)" + "# 4. What to do if $T$ is a long sequence? Use Overlap-add method!\n", + "\n", + "It applies a divide-and-conquer algorithm on convolution. We first need to learn about the \"linearity\" property in linear convolution and how it allows us to divide the problem into similar problems but with smaller sizes. Then, we learn how we can use circular convolution on those smaller problems, and combine their outputs." ] }, { @@ -240,7 +377,7 @@ "id": "1e45ccf7-1564-4c41-b32e-4de9b1b201a3", "metadata": {}, "source": [ - "## Linearity in Linear Convolution" + "### 4.1 The \"Linearity\" property in Linear Convolution" ] }, { @@ -248,15 +385,15 @@ "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", + "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", + "$$ 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", + "$$ 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] $$" + "$$ C[idx] = C^{(1)}[idx] + C^{(2)}[idx] $$" ] }, { @@ -264,22 +401,18 @@ "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 long $T$ by breaking $T$ into smaller parts." + "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": [] - }, { "cell_type": "code", "execution_count": null, - "id": "acf5516d-e52d-457c-a8c5-f4d1b7ae8f1c", + "id": "fecee5fc-3c10-4c2e-8c66-bce17c3650bd", "metadata": {}, "outputs": [], - "source": [] + "source": [ + "#WIP" + ] } ], "metadata": { From 083c8d8f75e81f85904cac2f5c2cbd1c0a366ee8 Mon Sep 17 00:00:00 2001 From: NimaSarajpoor Date: Mon, 7 Sep 2026 23:23:04 -0400 Subject: [PATCH 5/8] minor changes --- docs/overlap_add_math.ipynb | 45 +++---------------------------------- 1 file changed, 3 insertions(+), 42 deletions(-) diff --git a/docs/overlap_add_math.ipynb b/docs/overlap_add_math.ipynb index ecc5cf0..731676e 100644 --- a/docs/overlap_add_math.ipynb +++ b/docs/overlap_add_math.ipynb @@ -146,11 +146,11 @@ "id": "2c4176e2-60f3-410f-ad2d-6c1d29cca6cf", "metadata": {}, "source": [ - "Circular convolution is defined between two sequences that are both periodic and their periods are the same, say $N$. Their circular convolution is also a periodic sequence, with period $N$, and it can be computed as follows:\n", + "Circular convolution is defined between two sequences that are both periodic and their periods are the same, say $N$. Their circular convolution is also a periodic sequence, with period $N$. The following equation shows how it can be computed for the period in $range(0, N)$.\n", "\n", "$$C_{N}[idx] = \\sum_{i=0}^{N-1} \\tilde{T}[i] \\times \\tilde{Q'}[(idx-i)_{N}], \\quad 0 \\le idx \\le N-1$$\n", "\n", - "where, $C_{N}$ represents one period of convolution output. $\\tilde{Q'}$ and $\\tilde{T}$ are both periodic sequence with period $N$, and the index $(idx-i)_{N}$ simply means \"$idx-i$ modulo $N$\", i.e. $(idx-i) \\% N$, and it is always between $0$ and $N-1$." + "where, $C_{N}$ represents one period of convolution output. $\\tilde{Q'}$ and $\\tilde{T}$ are both periodic sequences, each with period $N$. The index $(idx-i)_{N}$ simply means \"$idx-i$ modulo $N$\", i.e. $(idx-i) \\% N$, and it is always between $0$ and $N-1$." ] }, { @@ -158,7 +158,7 @@ "id": "e10c5304-a1a1-4cbd-9b26-266e229bdc57", "metadata": {}, "source": [ - "The time complexity of calculating $C_{N}$ using the equation above is $O(N^{2})$. However, as will be shown later in the next section, the time complexity can be reduced to $O(NlogN)$. However, this means nothing for sliding dot product unless we understand how the linear convolution can be computed via circular convolution and how $N$ is related to $n$ and/or $m$. \n", + "The time complexity of calculating $C_{N}$ using the equation above is $O(N^{2})$. However, as will be shown later in the next section, the time complexity can be reduced to $O(NlogN)$. However, this means nothing for sliding dot product unless we understand how circular convolution can be used to compute the linear convolution, and how $N$ is related to $n$ and/or $m$. \n", "\n", "We first prove that the linear convolution between $Q'$ and $T$ is equivalent to $C_{N}$ between $\\tilde{Q'}$ and $\\tilde{T}$ if:\n", "\n", @@ -167,45 +167,6 @@ "* $\\tilde{T}_{N}$, i.e. one period of $\\tilde{T}$, is $T$ but with $N-n$ zero padding" ] }, - { - "cell_type": "markdown", - "id": "5ce64782-66d0-413b-a116-9ac3555966b6", - "metadata": {}, - "source": [ - "**Proof**\n", - "\n", - "> **Note:**
\n", - "> Keep in mind that $N = n + m - 1$, and:\n", - "> \n", - "> $$\n", - "\\tilde{Q'}_{N}[i] = \\begin{cases} \n", - " Q'[i] & \\text{if } 0 \\le i \\le m-1 \\\\\n", - " 0 & \\text{if } m \\le i \\le N-1\n", - " \\end{cases}\n", - "$$\n", - "> \n", - "> $$\n", - "\\tilde{T}_{N}[i] = \\begin{cases} \n", - " T[i] & \\text{if } 0 \\le i \\le n-1 \\\\\n", - " 0 & \\text{if } n \\le i \\le N-1\n", - " \\end{cases}\n", - "$$\n", - "\n", - "\n", - "\n", - "We start with rewriting the circular convolution equation by breaking the summation into two parts: \n", - "\n", - "$$C_{N}[idx] = \\sum_{i=0}^{n-1} \\tilde{T}[i] \\times \\tilde{Q'}[(idx-i)_{N}] + \\sum_{i=n}^{N-1} \\tilde{T}[i] \\times \\tilde{Q'}[(idx-i)_{N}]$$\n", - "\n", - "Note that:\n", - "* In the first term, $0 \\le i \\le n-1$, and hence $\\tilde{T}[i]$ is the same as $T[i]$.\n", - "* In the second term, $n \\le i \\le N-1$, and hence $\\tilde{T}[i]$ is zero.\n", - "\n", - "Hence, the equation becomes:\n", - "\n", - "$$C_{N}[idx] = \\sum_{i=0}^{n-1} T[i] \\times \\tilde{Q'}[(idx-i)_{N}]$$" - ] - }, { "cell_type": "markdown", "id": "a7e66371-4cd4-4673-957c-b7c548dcbfde", From f492bdf06d6119875847d784a66bd636c32272e4 Mon Sep 17 00:00:00 2001 From: NimaSarajpoor Date: Mon, 28 Sep 2026 21:34:08 -0400 Subject: [PATCH 6/8] update notebook --- docs/overlap_add_math.ipynb | 491 ++++++++++++++++++++++++++++++------ 1 file changed, 417 insertions(+), 74 deletions(-) diff --git a/docs/overlap_add_math.ipynb b/docs/overlap_add_math.ipynb index 731676e..59c6938 100644 --- a/docs/overlap_add_math.ipynb +++ b/docs/overlap_add_math.ipynb @@ -1,11 +1,19 @@ { "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": [ - "# Objective" + "# 0. Objective" ] }, { @@ -13,12 +21,11 @@ "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) provided different versions of MASS for computing the sliding dot product. At their core, they all come down to computing one/more convolution(s) via the frequency domain. To understand how convolution comes into the picture and how it is used for computing the sliding dot product, this notebook is created to answer the following questions:\n", + "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 relationship between sdp and linear convolution. Convolution is a well-studied area and understanding this relationship allows us to leverage convolution methods to compute the sliding dot product.\n", - "2. **How can we compute the convolution faster?** Once we understand the relationship between sdp and convolution, we can focus on learning/using method for faster convolution.\n", - "3. **Can we reduce the output size, and hence the computational load without losing the sdp?** Once we know how to get faster convolution, we can try to enhance its performance further by avoiding unnecessary calculation.\n", - "4. **What to do if $T$ is very long? Overlap-add method!** This is where we learn a divide-and-conquer algorithm on convolution, a useful method when input $T$ is large" + "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" ] }, { @@ -26,7 +33,7 @@ "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 name of arrays in the realm of 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 defined by\n", + "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,$$" ] @@ -36,7 +43,9 @@ "id": "3dda604a-3a3e-42f4-90c7-7b01ada16337", "metadata": {}, "source": [ - "# 1. How is the sliding dot product(sdp) related to (linear) convolution?" + "# 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." ] }, { @@ -44,7 +53,7 @@ "id": "877af78c-0026-44c0-8295-0f75db075d97", "metadata": {}, "source": [ - "### 1.1 Linear Convolution of Two Finite-Length Sequences $Q'$ and $T$" + "## 1.1 The Linear Convolution of Two Finite-Length Sequences $Q'$ and $T$" ] }, { @@ -68,9 +77,12 @@ "id": "8028ee30-3244-4c79-86ca-22e2b8876694", "metadata": {}, "source": [ - "where $T[.]$ and $Q'[.]$ are zeros for any index that is outside of their range. The linear convolution can be seen as a flip and slide operation, and at each index $idx$, the output 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. For instance, at index $idx=0$, there is only one overlapping element between $T$ and $Q'$, and that $T[0]$ and $Q'[0]$. And, $C[0] = T[0]Q'[0]$. Therefore, **the linear convolution is NOT the same as the sliding dot product** since sliding dot product is for cases when the overlap covers the full query and not just one element. \n", + "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 it is slided 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", - "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 revised as follows:\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. We will talk about sdp in the next section.\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 from $(-\\infty, \\infty)$ to $[0, n-1]$, i.e.\n", "\n", "$$ C[idx] = \\sum_{i=0}^{n-1}{T[i] \\times Q'[idx-i]}$$" ] @@ -80,15 +92,17 @@ "id": "593a8276-49cb-43d2-9300-8d29e67d3104", "metadata": {}, "source": [ - "Note that $Q'[idx-i]$ can be a non-zero value 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", + "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", + "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", - "So, when $idx$ changes from $0$ to $n+m-2$, there exist an $i$ such that both $T[i]$ and $Q'[idx-i]$ can be non-zero values. This shows that the linear convolution $C$ can have values at indices $\\set{0, 1, ..., n + m - 2}$. " + "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}$. \n", + "\n", + "Therefore, the length of linear convolution is $n + m - 1$." ] }, { @@ -96,7 +110,7 @@ "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$)" + "## 1.2 Sliding dot product (of $Q$ and $T$) as a slice of linear convolution (of $Q'$ and $T$)" ] }, { @@ -108,21 +122,30 @@ "\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[0]=T[0]Q'[1] + T[1]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}$\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}$\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", + "* $\\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$, $C[idx]$ actually becomes the dot product between $Q$ (the reverse of $Q'$) and $T_{idx-(m-1)}$. 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", + "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! The time complexity of calculating the linear convolution is $O(nm)$. Can we do better? " + "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." ] }, { @@ -130,7 +153,19 @@ "id": "8475be6a-9515-43ae-a8d7-083b8db3af51", "metadata": {}, "source": [ - "# 2. How can we compute the convolution faster than O(nm)?" + "# 2. How can we compute the convolution faster?\n", + "\n", + "In this section, we make three attempts, and each attempt is to help us 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." ] }, { @@ -138,7 +173,7 @@ "id": "52bcf566-3ac0-483d-81a3-57fa3b5766c3", "metadata": {}, "source": [ - "### 2.1 Convert linear convolution to circular convolution (Step I)" + "### 2.1.1 First Step: Convert linear convolution to circular convolution" ] }, { @@ -146,11 +181,14 @@ "id": "2c4176e2-60f3-410f-ad2d-6c1d29cca6cf", "metadata": {}, "source": [ - "Circular convolution is defined between two sequences that are both periodic and their periods are the same, say $N$. Their circular convolution is also a periodic sequence, with period $N$. The following equation shows how it can be computed for the period in $range(0, N)$.\n", + "**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 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", - "$$C_{N}[idx] = \\sum_{i=0}^{N-1} \\tilde{T}[i] \\times \\tilde{Q'}[(idx-i)_{N}], \\quad 0 \\le idx \\le N-1$$\n", + "where, the subscript $N$ represents one period. This operation is called \"Circular Convolution\".\n", "\n", - "where, $C_{N}$ represents one period of convolution output. $\\tilde{Q'}$ and $\\tilde{T}$ are both periodic sequences, each with period $N$. The index $(idx-i)_{N}$ simply means \"$idx-i$ modulo $N$\", i.e. $(idx-i) \\% N$, and it is always between $0$ and $N-1$." + "\n", + "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$." ] }, { @@ -158,13 +196,11 @@ "id": "e10c5304-a1a1-4cbd-9b26-266e229bdc57", "metadata": {}, "source": [ - "The time complexity of calculating $C_{N}$ using the equation above is $O(N^{2})$. However, as will be shown later in the next section, the time complexity can be reduced to $O(NlogN)$. However, this means nothing for sliding dot product unless we understand how circular convolution can be used to compute the linear convolution, and how $N$ is related to $n$ and/or $m$. \n", - "\n", - "We first prove that the linear convolution between $Q'$ and $T$ is equivalent to $C_{N}$ between $\\tilde{Q'}$ and $\\tilde{T}$ if:\n", + "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", - "* $N = n + m - 1$\n", - "* $\\tilde{Q'}_{N}$, i.e. one period of $\\tilde{Q'}$, is $Q'$ but with $N-m$ zero padding\n", - "* $\\tilde{T}_{N}$, i.e. one period of $\\tilde{T}$, is $T$ but with $N-n$ zero padding" + "* **Condition X:** $N = n + m - 1$\n", + "* **Condition Y:** $\\tilde{T}_{N}$ is $T$, padded with $N-n$ zeros\n", + "* **Condition Z:** $\\tilde{Q'}_{N}$ is $Q'$, padded with $N-m$ zeros" ] }, { @@ -172,11 +208,39 @@ "id": "a7e66371-4cd4-4673-957c-b7c548dcbfde", "metadata": {}, "source": [ - "Recall that the linear convolution equation is:\n", + "#### **Proof:** \n", "\n", - "$$ C[idx] = \\sum_{i=0}^{n-1}{T[i] \\times Q'[idx-i]}$$\n", + "Let's write down the equations for the linear convolution and the circular convolution below:\n", "\n", - "So, to show that $C_{N}$, from circular convolution, gives the same value as the linear convolution, we only need to show that $\\tilde{Q'}[(idx-i)_{N}]$ and $Q'[idx-i]$ give the same value when $N=n+m-1$. Let $j$ denote $idx-i$. So, all we need to do is to prove that $\\tilde{Q'}[(j)_{N}]$ and $Q'[j]$ give the same value for different values of $j$.\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", + "And 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 X:** $N = n + m - 1$\n", + "* **Condition Y:** $\\tilde{T}_{N}$ is $T$, padded with $N-n$ zeros\n", + "* **Condition Z:** $\\tilde{Q'}_{N}$ is $Q'$, padded with $N-m$ zeros\n", + " \n", + "**According to the condition X: $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 Y: $\\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", @@ -189,11 +253,16 @@ "\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$. hence, $j$ modulo $N$ is still $j$. Therefore: $\\tilde{Q'}[(j)_{N}]=Q'[j]$. Proof is now complete for this case.\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 Z, $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 Z, $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]=0$. So, all we need to do is to prove that $\\tilde{Q'}[(j)_{N}]$ is zero as well. \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", @@ -204,17 +273,21 @@ "\n", "$$min(j) = min(idx-i) = min(idx) + min(-i) = min(idx) - max(i) = 0 - (n-1)$$\n", "\n", - "So, in this case, $-(n-1) \\le j \\le -1$. \n", + "So: $-(n-1) \\le j \\le -1$. \n", "\n", - "Since $\\tilde{Q'}[(j)_{N}] == \\tilde{Q'}[(j+N)_{N}]$, we can compute the latter instead of the former. To do that, we first need to figure out the range for $j+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) \\le j \\le -1 \\implies -(n-1)+N \\le j+N \\le -1+N $$\n", + "$$ -(n-1)+N \\le j+N \\le -1+N $$\n", "\n", - "Let's plug the value of $N=n+m-1$ into $-(n-1)+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'}[(j)_{N}] == \\tilde{Q'}[(j+N)_{N}] = \\tilde{Q'}[j+N]$. Since the index $j+N$ is $\\ge m$, the value $\\tilde{Q'}[j+N]$ becomes 0 as the element resides in the zero-padding part. Proof is now complete for this case." + "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", + "---" ] }, { @@ -222,7 +295,7 @@ "id": "66a707c1-8004-4fe6-b954-ffe039eda441", "metadata": {}, "source": [ - "We just proved that the linear convolution can be computed in the form of circular convolution when the arrays have certain lengths and are zero-padded properly." + "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 show how this allows us to compute the linear convolution, and hence the sdp, faster!" ] }, { @@ -230,7 +303,7 @@ "id": "8648ff47-4d67-4bec-a91e-4ec6d799cc0e", "metadata": {}, "source": [ - "### 2.2 Compute Circular Convolution in the Frequency Domain (Step II)" + "### 2.1.2 Second Step: Compute Circular Convolution using Fast Fourier Transform (FFT)" ] }, { @@ -238,27 +311,21 @@ "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 $n$ and $m$ are the length of $T$ and $Q'$, respectively. However, the computation via the 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. So, instead of computing the circular convolution in time domain, i.e.\n", + "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 $n$ and $m$ are the length of $T$ and $Q'$, 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}[i] \\times \\tilde{Q'}[(idx-i)_{N}], \\quad 0 \\le idx \\le N-1$$\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", - "\\cdot\n", + "\\times\n", "FFT(\\tilde{Q'}_{N})\n", "\\right)$$\n", "\n", - "and its time complexity becomes $O(NlogN)$." - ] - }, - { - "cell_type": "markdown", - "id": "64b8a7ee-f1dd-4494-a9a4-4ba3293a3912", - "metadata": {}, - "source": [ - "So, if I can 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 than before!" + "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!!" ] }, { @@ -266,7 +333,9 @@ "id": "181cc549-9b2a-48ac-9077-71174cb49f1c", "metadata": {}, "source": [ - "# 3. Can we reduce the output size, and hence the computational load without losing the sdp?\n", + "## 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", @@ -274,15 +343,15 @@ "\n", "\n", "\n", - "Recall that the length of output is $n + m - 1$. The sliding dot product between $Q$ and $T$ is in $range(m-1, n)$ though. 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 we accordingly set $N$ to $n$ instead of $n+m-1$, does the slice in $range(m-1,n)$ still reflect the values of sdp? In other words, we need to show that\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 we accordingly set $N$ to $n$ instead of $n+m-1$, does the slice in range $[m-1,n)$ still reflects the values of sdp? In other words, we need to show that the linear convolution and circular convolution can give the same value for slice in range $[m-1, n)$ when $N=n$. In other words,\n", "\n", - "$$ C[idx] = \\sum_{i=0}^{n-1}{T[i] \\times Q'[idx-i]}, \\quad m-1 \\le idx \\le n-1$$\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'}[(idx-i)_{n}], \\quad m-1 \\le idx \\le n-1$$\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. To prove this, we just need to show that $Q'[idx-i]$ and $\\tilde{Q'}[(idx-i)_{n}]$ have the same value when $m-1 \\le idx \\le n-1$.
\n", + "should 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", @@ -295,12 +364,14 @@ "* $m-1 \\le idx \\le n-1$\n", "* $0 \\le i \\le n-1$\n", "\n", - "$$ j \\le max(j)=max(idx-i)=max(idx)+max(-i)=max(idx) - min(i)=(n-1)-0=n-1$$\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'}[(j)_{n}]$ is the same as $Q'[j]$. Proof is now complete for this case.\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'}[(j)_{n}]$ becomes zero as well in this case.\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", @@ -309,10 +380,15 @@ "* $0 \\le i \\le n-1$\n", "\n", " \n", - "* $\\tilde{Q'}[(j)_{n}] = \\tilde{Q'}[(j+n)_{n}]$. Note that: \n", - "$$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", + "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", - "So, when $j<0$, then: $m \\le j+n \\le n-1$. Therefore $\\tilde{Q'}[(j+n)_{n}]$ simply becomes $\\tilde{Q'}[j+n]$, and the value is zero as the index falls into the padded zeros. Therefore, $Q'[j]=\\tilde{Q'}[j+n]=0$. Proof is now complete for this case. " + "Proof is now complete for this case. \n", + "\n", + "---" ] }, { @@ -320,7 +396,21 @@ "id": "03cbd90d-a7e4-466d-af15-81d9a1c4ce6c", "metadata": {}, "source": [ - "We just showed that the sliding dot product can be computed via circular convolution with period $N=n$. As noted in previous section, this can be computed in $O(nlogn)$ which is faster than our baseline $O(nm)$ unless $m$ is small." + "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 Real FFT (RFFT) and Inverse Real FFT (IRFFT) when inputs are real-valued arrays. This should improve the performance of computing the ciruclar convolution when input sequences have only real values." ] }, { @@ -328,9 +418,9 @@ "id": "97fcc328-35f9-4b66-9598-1b5d8f1f97ac", "metadata": {}, "source": [ - "# 4. What to do if $T$ is a long sequence? Use Overlap-add method!\n", + "# 3. How to compute convolution faster when $T$ is very long?\n", "\n", - "It applies a divide-and-conquer algorithm on convolution. We first need to learn about the \"linearity\" property in linear convolution and how it allows us to divide the problem into similar problems but with smaller sizes. Then, we learn how we can use circular convolution on those smaller problems, and combine their outputs." + "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 and how it allows us to divide the problem into similar problems. Then, we learn how we can use circular convolution on those smaller problems, and combine their outputs." ] }, { @@ -338,7 +428,7 @@ "id": "1e45ccf7-1564-4c41-b32e-4de9b1b201a3", "metadata": {}, "source": [ - "### 4.1 The \"Linearity\" property in Linear Convolution" + "## 3.1 The \"Linearity\" property in Linear Convolution" ] }, { @@ -366,14 +456,267 @@ ] }, { - "cell_type": "code", - "execution_count": null, + "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 split $T$ at an index $s$ to get two sequences $T^{(1)}$ and $T^{(2)}$ as follows:\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 $T = T^{(1)} + T^{(2)}$." + ] + }, + { + "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 adding $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$. " + ] + }, + { + "cell_type": "markdown", + "id": "44608755-9b19-418a-b690-37da32ec4365", + "metadata": {}, + "source": [ + "Note that we have not improved the efficency 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 in improving 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 becomes 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": 12, + "id": "c042721c-dc3c-4e8d-a10b-0f7fcd07be8d", + "metadata": {}, "outputs": [], "source": [ - "#WIP" + "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", + "s = 2 * len(Q) # overlap-add\n", + "\n", + "ref = direct_convolve(Q, T)\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 optimal chunk size $s$" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "10b4039b-1f79-40b7-a30c-3fd67ae654b4", + "metadata": {}, + "outputs": [], + "source": [] } ], "metadata": { From d7bc922460d8a1b04360d3cc015fc14e0559e81b Mon Sep 17 00:00:00 2001 From: NimaSarajpoor Date: Mon, 28 Sep 2026 22:50:53 -0400 Subject: [PATCH 7/8] minor changes --- docs/overlap_add_math.ipynb | 75 ++++++++++++++++++------------------- 1 file changed, 37 insertions(+), 38 deletions(-) diff --git a/docs/overlap_add_math.ipynb b/docs/overlap_add_math.ipynb index 59c6938..42c86f4 100644 --- a/docs/overlap_add_math.ipynb +++ b/docs/overlap_add_math.ipynb @@ -77,12 +77,12 @@ "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 it is slided 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", + "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. We will talk about sdp in the next section.\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 from $(-\\infty, \\infty)$ to $[0, n-1]$, i.e.\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]}$$" ] @@ -100,9 +100,7 @@ "* 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}$. \n", - "\n", - "Therefore, the length of linear convolution is $n + m - 1$." + "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$." ] }, { @@ -145,7 +143,9 @@ "* 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." + "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)" ] }, { @@ -155,7 +155,7 @@ "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 improve the efficiency of computing sdp." + "In this section, we make three attempts, and each attempt is to help us eventually improve the efficiency of computing sdp." ] }, { @@ -181,14 +181,11 @@ "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 shows how the values of one period, for range $[0, N)$, are computed.\n", + "**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. This operation is called \"Circular Convolution\".\n", - "\n", - "\n", - "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$." + "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$." ] }, { @@ -198,9 +195,9 @@ "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 X:** $N = n + m - 1$\n", - "* **Condition Y:** $\\tilde{T}_{N}$ is $T$, padded with $N-n$ zeros\n", - "* **Condition Z:** $\\tilde{Q'}_{N}$ is $Q'$, padded with $N-m$ zeros" + "* **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" ] }, { @@ -217,18 +214,18 @@ "$$\\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", - "And 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 X:** $N = n + m - 1$\n", - "* **Condition Y:** $\\tilde{T}_{N}$ is $T$, padded with $N-n$ zeros\n", - "* **Condition Z:** $\\tilde{Q'}_{N}$ is $Q'$, padded with $N-m$ zeros\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 X: $N=n+m-1$**
\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 Y: $\\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", + "**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", @@ -255,8 +252,8 @@ "\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 Z, $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 Z, $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", + "* 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", @@ -295,7 +292,7 @@ "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 show how this allows us to compute the linear convolution, and hence the sdp, faster!" + "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!" ] }, { @@ -311,7 +308,7 @@ "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 $n$ and $m$ are the length of $T$ and $Q'$, 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", + "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", @@ -343,7 +340,7 @@ "\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 we accordingly set $N$ to $n$ instead of $n+m-1$, does the slice in range $[m-1,n)$ still reflects the values of sdp? In other words, we need to show that the linear convolution and circular convolution can give the same value for slice in range $[m-1, n)$ when $N=n$. In other words,\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", @@ -351,7 +348,7 @@ "\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", - "should 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", + "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", @@ -410,7 +407,7 @@ "\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 Real FFT (RFFT) and Inverse Real FFT (IRFFT) when inputs are real-valued arrays. This should improve the performance of computing the ciruclar convolution when input sequences have only real values." + "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 ciruclar convolution when input sequences have only real values." ] }, { @@ -418,9 +415,9 @@ "id": "97fcc328-35f9-4b66-9598-1b5d8f1f97ac", "metadata": {}, "source": [ - "# 3. How to compute convolution faster when $T$ is very long?\n", + "# 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 and how it allows us to divide the problem into similar problems. Then, we learn how we can use circular convolution on those smaller problems, and combine their outputs." + "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." ] }, { @@ -472,11 +469,12 @@ "source": [ "### 3.2.1 The overlap-add method: Core logic\n", "\n", - "Suppose $T$ is $\\{t_{0}, t_{1}, ..., t_{n-1}\\}$. I can split $T$ at an index $s$ to get two sequences $T^{(1)}$ and $T^{(2)}$ as follows:\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 $T = T^{(1)} + T^{(2)}$." + "Note that, at any index $i$, $T[i] = T^{(1)}[i] + T^{(2)}[i]$." ] }, { @@ -495,7 +493,7 @@ "id": "a014c77d-fdb1-4f42-ab13-cd6af26ae451", "metadata": {}, "source": [ - "The elements of $C$, the convolution between $T$ and $Q'$, can be computed by adding $C^{(1)}$ and $C^{(2)}$:\n", + "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", @@ -516,7 +514,7 @@ "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$. " + "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\"." ] }, { @@ -648,7 +646,7 @@ }, { "cell_type": "code", - "execution_count": 12, + "execution_count": 1, "id": "c042721c-dc3c-4e8d-a10b-0f7fcd07be8d", "metadata": {}, "outputs": [], @@ -692,9 +690,10 @@ "# input\n", "T = np.random.rand(1234)\n", "Q = np.random.rand(123)\n", - "s = 2 * len(Q) # overlap-add\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)" @@ -707,7 +706,7 @@ "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 optimal chunk size $s$" + "* Compute cost of operation and find an optimal value for the chunk size $s$" ] }, { From 64826b48ae2d4fe0aae8016f656a0fb591491885 Mon Sep 17 00:00:00 2001 From: NimaSarajpoor Date: Mon, 28 Sep 2026 22:57:50 -0400 Subject: [PATCH 8/8] fixed minor typos --- docs/overlap_add_math.ipynb | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/docs/overlap_add_math.ipynb b/docs/overlap_add_math.ipynb index 42c86f4..083953d 100644 --- a/docs/overlap_add_math.ipynb +++ b/docs/overlap_add_math.ipynb @@ -407,7 +407,7 @@ "\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 ciruclar convolution when input sequences have only real values." + "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." ] }, { @@ -522,7 +522,7 @@ "id": "44608755-9b19-418a-b690-37da32ec4365", "metadata": {}, "source": [ - "Note that we have not improved the efficency 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 in improving the efficiency! \n", + "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", @@ -535,7 +535,7 @@ "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 becomes zero. Therefore, we can just compute the convolution for non-zero elements and include its contribution to its corresponding slice in $C$. \n", + "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",