does the "convolution theorem" apply to weaker algebraic structures?

Question

The Convolution Theorem is often exploited to compute the convolution of two sequences efficiently: take the (discrete) Fourier transform of each sequence, multiply them, and then perform the inverse transform on the result.

The same thing can be done for convolutions in the quotient ring Z/pZ via the analogous Number Theoretic Transform.

Does this procedure generalize to other algebraic structures? Arbitrary rings? Semirings? Fast convolutions over the semiring (min,+) would be particularly useful.

I'm led to suspect this because both sorting and FFT can be computed using the same butterfly-like network with simple operations like "min", "max", "+", and "*" at the nodes of the network, and sorting can be thought of as a kind of convolution.

Answer 1

In general, it is a major open question in discrete algorithms as to which algebraic structures admit fast convolution algorithms and which do not. (To be concrete, I define the $(\oplus,\otimes)$ convolution of two $n$-vectors $[x_0,\ldots,x_{n-1}]$ and $[y_0,\ldots,y_{n-1}]$, to be the vector $[z_0,\ldots,z_{n-1}]$ with $$z_k = (x_0 \otimes y_k) \oplus (x_1 \otimes y_{k-1}) \oplus \cdots \oplus (x_k \otimes y_0).$$ Here, $\otimes$ and $\oplus$ are the multiplication and addition operations of some underlying semiring.)

For any $\otimes$ and $\oplus$, the convolution can be computed trivially in $O(n^2)$ operations. As you note, when $\otimes = \times$, $\oplus = +$, and we work over the integers, this convolution can be done efficiently, in $O(n \log n)$ operations.

But for more complex operations, we do not know efficient algorithms, and we do not know good lower bounds. The best algorithm for $(\min,+)$ convolution is $n^2/2^{\Omega (\sqrt{\log n})}$ operations, due to combining my recent APSP paper

Ryan Williams: Faster all-pairs shortest paths via circuit complexity. STOC 2014: 664-673

and

David Bremner, Timothy M. Chan, Erik D. Demaine, Jeff Erickson, Ferran Hurtado, John Iacono, Stefan Langerman, Perouz Taslakian: Necklaces, Convolutions, and X + Y. ESA 2006: 160-171

A substantially subquadratic algorithm for $(\min,+)$ convolution would (to my knowledge) imply a subcubic algorithm all-pairs shortest paths in general graphs, a longstanding open problem. The above ESA06 reference also gives a $O(n^2 (\log \log n)^2/\log n)$ algorithm for a "(median,+) convolution".

The situation is subtle. It's not clear when convolution over a semiring is easy and when it's hard. For instance, the $(\min,\max)$ convolution can be computed in subquadratic time: I believe that $O(n^{3/2} \log n)$ operations suffice. This can be obtained from adapting the $(\min,\max)$ matrix multiplication algorithm in my work with Vassilevska and Yuster on all-pairs bottleneck paths. Basically you reduce the problem to computing $\sqrt{n}$ instances of a $(+,\times)$ ring convolution.

Answer 2

You can do infimal convolution using the Legendre transform. You can certainly implement a reasonably fast algorithm for infimal convolution of convex piecewise linear functions though I don't know about the fastest known algorithms that don't use the Legendre transform. The Legendre transform is also called the Fenchel conjugate in optimisation papers so that's another good term to google on.

Answer 3

Look to the Chinese remainder theorem over rings. The discrete Fourier transform is just one example.

The discrete Fourier transform over a field $\mathbb{F}$ (for example $\mathbb{C}$ or $\mathbb{F}_{q}$) using an $n^{th}$ primitive root of unity $\alpha$ (in $\mathbb{C}$, $\alpha=e^{2\pi i/n}$) of a vector $\left[a_{n-1}\ a_{n-2}\ \ldots\ a_0\right]$ is equivalent to evaluating the polynomial $\sum^{(n-1)}_{j=0}a_jx^k $ at $\alpha^k$ for $k=0,1,\ldots, (k-1)$: $$\mathcal{F}\left( \sum^{(n-1)}_{j=0}a_jx^k\right)= \left[\sum^{(n-1)}_{j=0}a_j(\alpha^k)^j\right]_{j=0}^{n-1}$$ and the inverse DFT can be equated to (or easily derived from) Lagrangian interpolation.

Interpolation is a special case of the Chinese remainder theorem over $\mathbb{F}[x]$. Reducing a polynomial $f(x)$ modulo $(x-\alpha)$ is identical to evaluating $f(x)$ at $\alpha$ ($x\equiv \alpha \bmod{(x-\alpha)}$). Interpolation of the set of distinct points $\left\{(\alpha_k,\beta_k),\ k=0,1,\ldots,(n-1)\right\}$ is equivalent to finding the polynomial $f(x)\in \mathbb{F}[x]/\left(\prod_{k=0}^{n-1}(x-\alpha_k)\right)\mathbb{F}[x]$ such that $f(x)\equiv \beta_k\bmod{(x-\alpha_k)}$:$$\begin{array}{rl}f(x)&\equiv\sum_{k=0}^{n-1}\left(\prod_{j\neq k}(x-\alpha_j)\right)\left(\left(\prod_{j\neq k}(x-\alpha_j)\right)^{-1}\beta_k\bmod{(x-\alpha_k)}\right)\\ &\equiv \sum_{k=0}^{n-1}\left(\frac{\prod_{j\neq k}(x-\alpha_j)}{\prod_{j\neq k}(\alpha_k-\alpha_j)}\right)\beta_k \bmod{\prod_{k=0}^{n-1}(x-\alpha_k)}\end{array}$$ The second line of that equation, minus the mod, is the Lagrangian interpolation formula. The mod just means that any polynomial of the form $f(x)+h(x)\prod_{k=0}^{n-1}(x-\alpha_k)$, where $h(x)$ is any polynomial over $\mathbb{F}$, will pass through these same points.

For the DFT the final modulus product is $\prod_{k=0}^{n-1}(x-\alpha^n)=(x^n-1)$ and polynomial multiplication modulo $(x^n-1)$ is equivalent to the convolution of two $n$ long vectors.

The butterfly effect of the FFT comes from using a $2^n$-th primitive root of unity, $\alpha$. A similar effect (though not as elegant) can be achieved using $3^n$-th primitive roots of unity (or any other $p^n$) and combining specified sets of three ($p$).

What the Chinese remainder theorem is really giving you is an isomorphism between the factor rings $R/\bigcap_{k=0}^{n-1}I_k\cong \bigotimes_{k=0}^{n-1}R/I_k$ where $I_k<R$ are ideals in $R$ with special properties. If $R$ is a Euclidian domain, such as $\mathbb{Z}$ or $\mathbb{F}[x]$, these special properties amount to the ideals being generated by relatively prime elements. For example over $\mathbb{Z}$ the ideals would look like $I_k=p_k\mathbb{Z}$ with $GCD(p_k,p_j)=1$ for all distinct pairs of $(k,j)$, and the intersection of the ideals would be $\bigcap I_k=\left(\prod p_k\right)\mathbb{Z}$. Any operations which done modulo the product can be performed in parallel modulo each factor.

Answer 4

I stumbled onto this question after encountering the problem in Bayesian inference (specifically, performing max-product inference on random variables where $M = L + R$ where $L$ and $R$ (or, an an alternate form of the problem, $M$ and $L$) have known discrete distributions and the distribution on $M$ (or $R$, in the alternative problem) is sought. The probability mass function (PMF) on $M$ that is achieved by maxing out $L$ and $R$ will be the result of the max-convolution between the PMFs of $L$ and $R$.

After playing around with it a bit, I noticed that you can rewrite the problem by generating $u^{(m)}$ for each index $m$ of the result (where $u^{(m)}[\ell] = L[\ell] R[{m-\ell}]$). When the vectors are nonnegative (in my case, this was true because they are PMFs), then you can perform the $\max_\ell u^{(m)}_\ell$ with the Chebyshev norm:

$$ M[m] = \max_\ell u^{(m)}_\ell \\ = \lim_{p \to \infty} \| u^{(m)} \|_p \\ = \lim_{p \to \infty} {\left( \sum_\ell {\left( u^{(m)}[\ell] \right)}^{p} \right)}^{\frac{1}{p}} $$

This in turn can be approximated by substituting the $p^*$-norm, where $p^*$ is a sufficiently large constant (essentially, this exploits the fact that the $p^*$-norm does have inverse operations, whereas the $\max$ does not, but the $p^*$-norm can still approximate the $\max$): $$ \approx {\left( \sum_\ell {\left( u^{(m)}[\ell] \right)}^{p^*} \right)}^{\frac{1}{p^*}}\\ = {\left( \sum_\ell {L[\ell]}^{p^*} ~ {R[{m-\ell}]}^{p^*} \right)}^{\frac{1}{p^*}}\\ = {\left( \sum_\ell {\left(L^{p^*}\right)}[\ell] ~ {\left(R^{p^*}\right)}[{m-\ell}] \right)}^{\frac{1}{p^*}}\\ = {\left( L^{p^*} ~*~ R^{p^*} \right)}^{\frac{1}{p^*}}[m] $$

where $*$ is standard convolution, which can be performed in $n \log(n)$ via FFT. A short paper illustrating the approximation for large-scale probabilistic inference problems is in press at the Journal of Computational Biology (Serang 2015 arXiv preprint).

Afterward a Ph.D. student, Julianus Pfeuffer, and I hacked out a preliminary bound on the absolute error ($p^*$-norm approximations of the Chebyshev norm are poor when $p^*$ is small, but on indices where the normalized result $\frac{M[m]}{\max_{m'} M[m']}$ is very small, large $p^*$ can be numerically unstable). Julianus and I worked out a modified method for cases when the dynamic range of the result $M$ is very large (if not, then the simple method from the first paper works fine). The modified method operates piecewise over $\log(\log(n))$ different $p^*$ values (including runtime constants, this amounts to $\leq 18$ FFTs even when $n$ is on the scale of particles in the observable universe, and those 18 FFTs can be done in parallel). (Pfeuffer and Serang arXiv link which has a link to a Python demonstration of the approach) The approach (and the Python demo code) generalize to tensors (i.e., numpy.arrays) by essentially combining the element-wise Frobenius norm with multidimensional FFT.

Oliver Serang