A naive and easy formula to calculate convolution power

Question

In mathematics (in particular, functional analysis) convolution is a mathematical operation on two functions ($f$ and $g$) to produce a third function that expresses how the shape of one is modified by the other. Convolution of two function $f, \ g$ is defined as follows

$$(f*g)(t)\triangleq \ \int _{-\infty }^{\infty }f(t-\tau )g(\tau )\,d\tau$$

And discrete case of convolution of $f,g$ is given by

$$(f*g)[n]=\sum_{m=-\infty}^{\infty}f[n-m]g[m].$$

In case of $f=g$ the convolution power of $f^{*2}$ taking place. Convolution power of function $x$ is the n-fold iteration of the convolution with itself. Thus if $x$ is a function on Euclidean space $R^d$ and $n$ is a natural number, then the convolution power is defined by $$x^{*n}=\underbrace{x*x*\cdots*x}_{n \ \mathrm{times}}.$$ Assume we have the real function $f$, for example, let be the power function $f=x^k$ where $k$ is natural. Is there a generzlized reccurence formula (i.e by power rule $f^{*n}=f*f^{*n-1}$ represent the $n$-th power convolution via $n-1$-th convolution for every $n\geq 0$) for calculation both discrete and infinitesimal $n$-th convolution power of $f$ ?

Answer 1

A more explicit formula may be given by the Fourier transform. Let me start from the continuous case. Assuming $f\in L^1(\mathbb{R})$ (otherwise $f\star f$ may not even be defined), then we have $$\hat{f}(\xi)=\int_{-\infty}^{+\infty}f(x)e^{-ix\xi}\,dx \qquad \xi \in \mathbb{R}$$ Then it is well-known that the Fourier transform turns convolution into product, so that if $f,g\in L^1$ then $$\widehat{f\star g}=\hat{f}\cdot \hat{g}\Rightarrow f\star g=\mathcal{F}^{-1}(\hat{f}\cdot \hat{g}) $$ Where $\mathcal{F}^{-1}$ denotes the inverse Fourier transform. In particular, for $f\in L^1$ we have $$f^{*n}(x)=\mathcal{F}^{-1}(\hat{f}^n)(x)=\frac{1}{2\pi}\int_{-\infty}^{+\infty}\hat{f}(\xi)^ne^{ix\xi}\,d\xi \qquad x\in \mathbb{R}$$ and this is an explicit formula, provided that you can compute the two integrals.

In the discrete case everything still works, except the Fourier transform of a function $f:\mathbb{Z}\to \mathbb{R}$ with $\sum_{n=-\infty}^{+\infty}|f(n)|<\infty$ is in this case its Fourier series expansion $$\hat{f}(\theta)=\sum_{n=-\infty}^{+\infty}f(n)e^{in\theta}\qquad \theta \in [0,2\pi] $$ and the inverse is the expression of the Fourier coefficients, so that $$f^{*n}(m)=\frac{1}{2\pi}\int_{0}^{2\pi}\hat{f}(\vartheta)^ne^{-im\vartheta}\,d\vartheta\qquad m\in \mathbb{Z} $$

Answer 2

There exists a formula, I think it was first presented in De Pril, N. (1985). Recursions for Convolutions of Arithmetic Distributions.

Given $f:\{0,1,2,\ldots\}\to\mathbb R^+$ with $f[0]>0$, and $g=\underbrace{f\circledast f \circledast \cdots \circledast f}_{n \text{ times}}$ one can exploit generating functions to get \begin{equation} g[k] = \begin{cases} \left(f[0]\right)^n & k=0\\ \frac{1}{f[0]} \sum_{j=1}^k \left(\frac{n+1}{k}j-1\right)f[j] g[k-j] & k\ge 1 \end{cases} \end{equation}

This lets you compute the values of $g$ incrementally: first $g[0]$, then $g[1]$ based on $g[0]$, then $g[2]$ based on $g[0]$ and $g[1]$, and so on.

Comment on the practical use of this method: Let's consider the case where the support of $f$ is bounded, say $\{0,1,2,\ldots,N\}$. I'm afraid the computational cost of using the formula above is the same if not higher than naive repeated convolution of $f$ with itself $n$ times (disregarding $k$, that would be $\mathcal O (N^2)$ in both cases). However if, for instance, you already know that $g[k]$ becomes negligibly small for $k$ greater than some small $k^*$, then you can truncate the recursion at $k^*$ and save some computation.