### Tutorial on Fast Convolution

1. Introduction
2. Convolution
3. Fourier Transform and Convolution
4. Max- and P-Norm Convolution
5. Literature

#### Introduction

This is a tutorial on how to implement the convolution operation or faltung on discretized input. We will outline standard, p-norm, and max-norm convolution and see how to utilize the FFT to reduce the computational runtime from cubic to $\mathcal{O}(n\log n)$ for p-norm and max-norm convolution.

Convolution is a kernel function with various applications in arithmetics, probability theory, signal processing in physics, electrical engineering, and image processing in computer vision. The recent development of convolutional neural networks ( see tutorial on deep learning) for learning image data is one of the most exciting applications.

#### Convolution

Convolution, also known as faltung or superposition integral, is a binary operator which measures the amount of overlap of a continuous function $\phi$ as it is shifted over another function $\psi$. $$\phi *\psi = \int_{-\infty}^{\infty}\phi(\tau)\psi(t-\tau) d\tau$$ For many real-world examples the ranges are finite and $\phi$ and $\psi$ are either unknown or not solvable analytically in an efficient manner. In the first case we may have a set of measured values, in the latter case, we can easily transform the signal into a discrete series by sampling from the continuous function. We therefore continue with the discrete version of the convolution. Assume we have two input vectors $L\in\mathbb{C}^m$ and $R\in\mathbb{C}^n$, the standard convolution $*$ is defined as $$(L*R)[m] \equiv \sum_lL[m]\cdot R[m-l].$$ Many authors split the above convolution operation at position $m$ logically into two parts. In the first part a vector $u$ is constructed as the point-wise product of every non-zero overlap position of input vectors $L$ and reversed $R$ (see Figure bla). $$u[m] = L[0:m]\cdot R[m:-1:0]$$ In the second part vector $u$ is reduced to scalar. The type of utilized reduction operator gives the convolution its precise name: instead of the sum we may plug in a p-norm or a Tchebychev norm.
Reduction functionConvolution Type
$\sum_i u_i$standard
$\left\lVert u \right\rVert _p$p-norm
$\left\lVert u \right\rVert _{\infty} = \max (u)$max-norm or Tchebychev norm
We have $|L|+|R|-1$ overlap shift positions and a maximal overlap size of $\min (|L|, |R|)$. Let $n$ be the size of the larger vector, the costs for constructing a single $u[m]$ are $O(n)$. Constructing $\vec{u}$ for all positions $l\in [0:|L|+|R|-2]$ the naive way is therefore $O(n^2)$ In the standard convolution we compute the scalar product for all overlaps, i.e. we sum up $u[m]$: $$(L*R)[m] = \sum u[m] = \sum_lL[m]\cdot R[m-l]$$ Any linear operation on $u$ will not change the runtime. For the max and p-norm convolution the addition is replaced by the maximum or p-norm function. The exact method is to compute a $u$ vector for each overlap shift of $L$ and $R$ and determine the maximum. But even omitting explicit storage will not help us to improving runtime.

#### Fourier Transform and Convolution

The convolution theorem states that under suitable conditions the Fourier transform of a convolution is the pointwise product of Fourier transforms: $$\mathcal{F}\{L*R\} = \mathcal{F}\{L\}\cdot \mathcal{F}\{R\}$$ Applying the inverse of the Fourier transform yields $$L*R = \mathcal{F}^{-1}\big \{ \mathcal{F}\{L\}\cdot \mathcal{F}\{R\}\big \}$$

Together with zero padding to the next power of 2, we can use the Cooley-Tuckey FFT to compute the convolution in $\mathcal{O}(n\log n)$. The circular convolution of two vectors is then found by taking an FFT of each, and then applying an inverse FFT on their point-wise product.

#### Max- and P-Norm Convolution

In order to utilize the faster runtime of the FFT, note that FFT always fuses information in contrast to a maximum entry selection where information gets lost. On the one hand we need this feature of the FFT, on the other we want that larger entries in $L$ and $R$ to dominate. The trick is to take the power and revert it at the end of the computation. Maximizing $p$ is exactly the Chebyshev norm or supremum norm: $$max(u) = \lvert\lvert u\rvert\rvert_\infty = \lim_{p\rightarrow\infty}|| u||_p$$ with $||u||_p$ defined as $(\sum u^p)^{1/p}$ in the discrete case. We can get close (error in range of machine precision) to the exact result, by chosing a $p^*$ large enough, s.t. $||u||_\infty \in \big [ ||u||_\infty-\epsilon; ||u||_\infty+\epsilon\big ]$ (see Pfeuffer et al. on how to choose p).

The code below demonstrates how to use the faster p-norm convolution to approximate a max convolution: