# Efficient computation of “discrete infimal convolution”

This question arises from an application to graphical models in probability theory, but I have abstracted that part out so only algebra remains. Let $\mathbb{R}$ denote standard field of real numbers and $\mathcal{R} = (\mathbb{R}\cup\{\infty\},\min,+)$ the tropical semiring. Let $\mathbb{Z}_n$ denote the cyclic group of order $n$.

Elements of the group ring $\mathbb{R}\mathbb{Z}_n$ are tuples $(x_0,\ldots,x_{n-1})$ and multiplication of these corresponds to discrete cyclic convolution. The Fast Fourier Transform gives an embedding $\mathbb{R}\mathbb{Z}_n\to\mathbb{C}^n$ (with elementwise sum and product). The FFT and its inverse can be computed in $O(n\log n)$ arithmetic operations, so elements of $\mathbb{R}\mathbb{Z}_n$ can be multiplied in $O(n\log n)$ arithmetic operations.

Define the group semiring $\mathcal{R}\mathbb{Z}_n$ in an analogous way. The product of $x = (x_0,\ldots, x_{n-1})$ and $y = (y_0,\ldots,y_{n-1})$ in this semiring is given by $(x\cdot y)_k = \min _{j \in \mathbb{Z}_n} (x_j + y_{k-j})$. This operation could perhaps be called "discrete cyclic infimal convolution" by analogy with the notion of infimal convolution in convex analysis. I'm not sure whether there is a more standard name -- this one does not pop up in quick google searches.

The naive way of computing a discrete infimal convolution uses $O(n^2)$ operations, just as the naive method for computing a standard cyclic convolution. My question is: is there a way to compute this "discrete cyclic infimal convolution" in $O(n\log n)$ arithmetic operations?

In convex analysis, there is an analog of the Fourier transform which turns infimal convolution into pointwise addition: the convex conjugate (or Fenchel or Legendre transformation). However, this operation only behaves nicely for convex functions, so it is not clear to me how one would translate it to an equivalent tool for $\mathcal{R}\mathbb{Z}_n$, but perhaps there is something there.

I would be interested in answers to the question regardless of whether they go through some analog of the Fourier transform. Also, I don't mind various restrictions such as making $n$ a power of $2$, replacing $\mathbb{R}$ with $\mathbb{Q}$, removing $\infty$ from the definition of the semiring, etc. Really anything on this theme would be helpful. Also any suggestions for better tags would be appreciated; perhaps this is a well-studied area I'm not aware of.

You can approximate this using a fast numerical method (note this answer assumes you are performing max-convolution on an equivalent transformed problem on the ring $(\times, \max)$ rather than $(+, \min)$-- you can convert from this one to the one in the question using logarithms and negation if necessary):

If $M[m]$ is the exact result at index $m$ (i.e., $M = L ~*_\max~ R$), then consider every vector $u^{(m)}$ (for each index $m$ of the result, where $u^{(m)}[\ell] = L[\ell] R[{m-\ell}]$). When the vectors $L$ and $R$ are nonnegative (in my case, this was true because they are full of probabilities), 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}}\\ \approx {\left( \sum_\ell {\left( u^{(m)}[\ell] \right)}^{p^*} \right)}^{\frac{1}{p^*}}\\$$

(where $p^*$ is a large enough constant)

$$= {\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]$$

The standard convolution (denoted $*$) 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 that is numerically stable in cases when the dynamic range of the result $M$ is very large (when the dynamic range is small, 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 linked by the second paper) generalizes to tensors (i.e., numpy.arrays) by essentially combining the element-wise Frobenius norm with multidimensional FFT. There is no other known faster-than-naive method for matrices (and tensors).

Max-convolution is a very fun puzzle-- I'd be interested to hear about your specific application for this in comments (maybe the fast numerical method would work for you?).

Oliver Serang

• Thanks. This is one of the methods I thought about, but it seemed difficult to prove anything about in an appropriate computational model; I think this is related to the numerical instabilities. For example, if the data are all rational numbers, can we set $p^*$ appropriately and compute everything to some finite precision nicely bounded in terms of the data, then use some rounding scheme to recover the exact result? – Noah Stein Jun 29 '15 at 13:50
• For me the problem was mostly a curiosity. Some of my coworkers develop a software package (Dimple) that does belief propagation and other inference tasks, and it occurred to me that it was slightly strange that for the sum-product algorithm certain messages could be computed efficiently using FFT convolutions, but there wasn't a clear analog for the min-sum algorithm because that would require a fast, stable max-convolution algorithm. It's frustrating that numerical stability comes into it when the FFT is so nice on that front (as far as I know). – Noah Stein Jun 29 '15 at 13:53
• Hi Noah, faster max-product belief prop. (small world) is exactly why I created this trick. Note that the speed & novelty compared to extant methods lies in 1) not only using an approximation ($\| \cdot \|_{p^*}$), but also computing the $p^*$-approximation with ordinary FFT 2) it can be made numerically stable (see second paper); an error bound between $\|\cdot \|_\infty$ and $\| \cdot \|_{p^*}$ can be derived (and in practice, the error is very small). The code from the papers is brief and w/ Apache 2 license, so feel free to use it in Dimple (and enjoy a pastry from Flour bakery for me). -O – Oliver Jun 29 '15 at 21:40

Hi Noah,

Nice question! I tried thinking of the discrete fourier transform as a change of basis to the eigenbasis of a circulant matrix, but this did not generalize well tropically.

However, there are work on generalized Legendre-Fenchel transform to non convex/concave functions, including discrete ones, which preserves the convolution-to-sum property. See, for example

Characterization of a simple communication network using legendre transform

or

Slope transforms: theory and application to nonlinear signal processing

As to fast computation of the discrete Legendre-Fenchel transform, see Lucet's thesis:

La transformee de legendre-fenchel etla convexifiee d'une fonction: algorithmes rapides de calcul, analyse et regularite du second ordre

I haven't read the last two papers very carefully, but hopefully they're relevant. As a small observation, you can order either $y$ or $x$ in decreasing order, thus at least one of the Legendre transform is straight forward.

• Thank you for these references. However, looking through them, they all seem to work with versions of the Legendre transform $^\ast$ which retain the standard property that $f^{\ast\ast}$ is the convex hull of $f$ (the function whose epigraph is the convex hull of the epigraph of $f$). As such, there doesn't seem to be a way to go back from $f^\ast + g^\ast$ to the infimal convolution of $f$ and $g$ when $f$ and $g$ are not convex. Am I missing something? – Noah Stein Feb 13 '13 at 13:05
• I just noticed that you are at Berkeley and working with Bernd Sturmfels, so I should mention that I had already emailed him about this yesterday before you replied. Please let me know if you two have any further thoughts about this. – Noah Stein Feb 13 '13 at 16:06
• Good point. They did mention the extended inverse Legendre transform in the first paper (and generally searching for "slope transform deconvolution" returns some papers with claims that it is possible to deconvolve for non-convex/concave functions, but I have yet to find something explicit). There are nice pictures in www2.ensc.sfu.ca/~ljilja/cnl/presentations/takashi/seminar.ps which explain the intuition of the extended Lengedre transform. Given that it keeps the information on all slopes, I would expect deconvolution to be possible. – Ngoc Mai Tran Feb 13 '13 at 23:09
• Ah yes, I learned of this problem from Bernd. – Ngoc Mai Tran Feb 13 '13 at 23:52

After a bit more google searching, I found the 2006 paper Necklaces, Convolutions, and X+Y, which addresses this problem. The nine authors give an $O(n\sqrt{n})$ algorithm in the nonuniform linear decision tree model (I'm having little trouble pinning down the details of this computational model) and an $O\left(\frac{n^3(\log \log n)^3}{(\log n)^2}\right)$ algorithm in the real RAM model. I haven't found any newer results, so it seems that the problem I posed is open.

• Regarding the etiquette of answering my own question in this way: do I accept this answer because it shows the problem is considered open, which is as far as MO is intended to go? Or not because I would still be happy for someone to answer it? – Noah Stein Feb 18 '13 at 15:03
• Here's a paper on infimal convolution with a pretty excellent review of the literature, asserting that the problem is open (as of 2015). (The paper's actual contribution is the same $p$-norm convolution trick that Oliver and others have suggested.) – Bill Bradley Feb 24 '16 at 4:16
• @BillBradley: Thanks. By the way, that paper is by the same Oliver. :-) – Noah Stein Feb 24 '16 at 17:07