Numerical partial differentiation of a convolution product with FFT

Question

How can one numerically calculate the partial derivatives of a convolution function, particularly when the closed-form or analytical expressions of the derivatives are not readily available? I am dealing with the convolution of two functions, where at least one does not have an easily computable derivative in the convolution. I am considering a numerical approach involving Fourier transforms (FFT) to simplify the convolution operation, followed by a perturbation method for derivative estimation.

Specifically, let's denote the convolution of two functions $f$ and $g$ as $h = f * g$, where $f$ and $g$ depend on parameters $\mathbf{p} = [p_i, p_j, p_k, p_l]$. I am interested in computing the partial derivative of $h$ with respect to a specific parameter, say $p_i$, numerically. Suppose we do the following:

1. Performing FFT on $f$ and $g$ to obtain $\hat{f}$ and $\hat{g}$, respectively.
2. Modifying the parameter $p_i$ by a small amount $h$ to get $p_i + h$ and then recalculating the FFT of the modified function.
3. Using the difference in inverse FFT results to numerically approximate the partial derivative with respect to $p_i$.

Example: For a Lorentzian function $L(x; p_i, p_j)$ and a Gaussian function $G(x; p_k, p_l)$, where $p_i, p_j, p_k, p_l$ are parameters of these functions, how can we accurately and efficiently compute $\frac{\partial}{\partial p_i}(L * G)$ using FFT and the described numerical method?

How would we go about calculating the Hessian elements numerically?

Is this approach viable for accurately calculating partial derivatives, especially for complex convolutions like that of a Lorentzian with a Gaussian, where analytical solutions are challenging? Are there any pitfalls or considerations to be aware of?

Answer 1

This "answer" doesn't address the question specifically, but outlines what might be a possible alternative to the numerical approach in the question.

---

Consider the Lorentzian function

$$L(x)=\frac{1}{\pi} \frac{\frac{1}{2} \Gamma}{(x-x_0)^2+\left(\frac{1}{2} \Gamma\right)^2}\tag{1}$$

with Fourier transform

$$\hat{L}(\omega)=\int\limits_{-\infty}^\infty L(x)\, e^{-2 \pi i \omega x}\, dx=e^{-2 \pi i x_0 \omega-\pi \Gamma |\omega|}\tag{2}$$

and Gaussian function

$$G(x)=\frac{1}{\sigma \sqrt{2 \pi}} e^{-\frac{(x-\mu)^2}{2 \sigma ^2}}\tag{3}$$

with Fourier transform

$$\hat{G}(\omega)=\sqrt{\frac{1}{\sigma^2}}\, \sigma\, e^{-2 \pi i \mu \omega-2 \pi^2 \sigma^2 \omega^2}\tag{4}.$$

---

By the convolution theorem, the Fourier transform of the convolution

$$h(x)=[L(y) * G(y)](x)=\int\limits_{-\infty}^{\infty} L(y)\, G(x-y)\, dy\tag{5}$$

is

$$\hat{h}(\omega)=\hat{L}(\omega)\, \hat{G}(\omega)=\sqrt{\frac{1}{\sigma^2}} \sigma \exp\left(-2 i \pi (\mu+x_0) \omega-2 \pi^2 \sigma^2 \omega^2-\pi \Gamma |\omega|\right)\tag{6}.$$

---

Mathematica gives the inverse Fourier transform

$$h(x)=\int\limits_{-\infty}^\infty \hat{h}(\omega)\, e^{2 \pi i x \omega}\, d\omega=\frac{e^{-\frac{(i \Gamma-2 \mu+2 x-2 x_0)^2}{8 \sigma^2}} \left(1+i\, \text{erfi}\left(\frac{i \Gamma-2 \mu+2 x-2 x_0}{2 \sqrt{2} \sigma}\right)\right)+e^{-\frac{(i \Gamma+2 \mu-2 x+2 x_0)^2}{8 \sigma^2}} \left(1-i\, \text{erfi}\left(\frac{-i \Gamma-2 \mu+2 x-2 x_0}{2 \sqrt{2} \sigma}\right)\right)}{2 \sqrt{2 \pi} \sigma}\tag{7}$$

which I believe can be represented as the "nested" exponential Fourier series

$$h(x)=\underset{N, f\to\infty}{\text{lim}}\left(\sum\limits_{n=1}^N \frac{\mu(2 n-1) }{2 n-1} \left(\frac{3 F(0)}{8}\\+\frac{1}{2} \sum\limits_{k=1}^{2 f (2 n-1)} (-1)^k\, \cos\left(\frac{\pi k}{2 n-1}\right) \left(F\left(\frac{k}{4 n-2}\right) e^{\frac{i \pi k x}{2 n-1}}+F\left(-\frac{k}{4 n-2}\right) e^{-\frac{i \pi k x}{2 n-1}}\right)\\-\frac{1}{8} \sum\limits_{k=1}^{4 f (2 n-1)} (-1)^k\, \left(F\left(\frac{k}{8 n-4}\right) e^{\frac{i \pi k x}{4 n-2}}+F\left(-\frac{k}{8 n-4}\right) e^{-\frac{i \pi k x}{4 n-2}}\right)\right)\right)\tag{8}$$

where $\mu(n)$ is the Möbius function, the evaluation frequency $f$ in the two inner sums over $k$ is assumed to be a positive integer, and $F(\omega)=\hat{h}(\omega)$ defined in formula (6) above.

---

My related MSO question provides information on the derivation of formula (8) above (which corresponds to formula (5) in my related MSO question).

---

Now consider the derivatives

$$\frac{\partial^m\, \hat{h}(\omega)}{\partial\, \Gamma^m}=(-\pi\, |\omega|)^m\, \hat{h}(\omega)\tag{9}$$

$$\frac{\partial^m\, \hat{h}(\omega)}{\partial\, x_0^m}=(-2 \pi i \omega)^m\, \hat{h}(\omega)\tag{10}$$

$$\frac{\partial^m\, \hat{h}(\omega)}{\partial\, \mu^m}=(-2 \pi i \omega)^m\, \hat{h}(\omega)\tag{11}$$

and note that formulas (10) and (11) above give the same result.

---

The $m^{th}$ order derivative of $\hat{h}(\omega)$ with respect to $\sigma$ is a bit more complicated, but the first few are given in Table (1) below, and the general result can be defined recursively.

---

Table (1): $m^{th}$ order derivatives of $\hat{h}(\omega)$ with respect to $\sigma$

$$\begin{array}{cc} m & \frac{\partial^m\,\hat{h}(\omega)}{\partial\, \sigma^m} \\ 1 & -4 \pi^2 \sigma \omega^2\, \hat{h}(\omega) \\ 2 & 4 \pi^2 \omega^2 \left(4 \pi^2 \sigma^2 \omega^2-1\right) \hat{h}(\omega) \\ 3 & \left(48 \pi^4 \sigma \omega^4-64 \pi^6 \sigma^3 \omega^6\right) \hat{h}(\omega) \\ 4 & 16 \pi^4 \omega^4 \left(16 \pi^4 \sigma^4 \omega^4-24 \pi^2 \sigma^2 \omega^2+3\right) \hat{h}(\omega) \\ \end{array}$$

---

Now for example, this seems to suggest one can evaluate $\frac{\partial^m\, h(x)}{\partial\, \Gamma^m}$ using the nested Fourier series representation defined in formula (8) above where $F(\omega)$ is given by formula (9) above. For illustration purposes consider the specific case

$$f(\Gamma)=\frac{\partial h(x)}{\partial\, \Gamma}=\\\frac{i \left(-\sqrt{2 \pi} e^{\frac{1}{8} (\Gamma+i)^2} (\Gamma+i) \left(\text{erfi}\left(\frac{1-i \Gamma}{2 \sqrt{2}}\right)+i\right)+\sqrt{2 \pi} e^{\frac{1}{8} (\Gamma-i)^2} (\Gamma-i) \left(\text{erfi}\left(\frac{1+i \Gamma}{2 \sqrt{2}}\right)-i\right)+8 i\right)}{16 \pi},\quad x=0\land x_0=\frac{1}{2}\land\sigma=1\land\mu=-1\tag{12}.$$

---

Figure (1) below illustrates the nested exponential Fourier series representation defined in formula (8) above in orange overlaid on the blue reference function $f(\Gamma)$ defined in formula (12) above where formula (8) above is evaluated using the evaluation limits $f=4$ and $N=20$ and $F(\omega)=\frac{\partial\, \hat{h}(\omega)}{\partial\, \Gamma}=-\pi |\omega|\, \hat{h}(\omega)$ (using the same parameter assignments specified at the end of formula (12) above).

---

Figure (1): Illustration of representation of $f(\Gamma)$ derived from formula (8) above in orange overlaid on the blue reference function $f(\Gamma)$ defined in formula (12) above.