Pavel Panchekha


Share under CC-BY-SA.

The Gaussian and the Sine

My favorite math problem is the following:

A Gaussian and a sinusoid have very similar shapes. The big difference is that the Gaussian has one “hump”, while the sinusoid has infinitely many. Prove the following formalization of this intuition:

Let \(G(x, \mu)\) be the standard Gaussian \(e^{-(x-\mu)²/2}\). Then the infinite sum \[\sum_{k\in\mathbb{Z}} G(x, k \pi)\] of shifted Gaussians is just a shifted sinusoid \[\alpha \cos(x) + \beta.\]

The theorem demanded is not outrageous. The Gaussian is a variant on the exponential function, as is the sinusoid. There are some geometric series to immediately make use of. The theorem seems within reach of some clever algebra.

Nonetheless, it makes sense to test the claim on a computer.

The two curves diverge because I have added up only three shifted Gaussians, not infinitely many of them; the red curve is the sum of Gaussians, and the blue curve is the properly-shifted sinusoid. They clearly match; the error is small enough to dismiss as floating-point artifacts.

The computer affirms what we already suspected: the theorem must be a few clever algebraic manipulations away. But the truth is not so.

In fact, the theorem is false, and the graph above a vicious lie. The graph hides an error term of order 0.0003; to see why, we need to go deeper into the problem.

It’s simplest to talk about equivalence to a sinusoid in frequency space. So, let’s try to understand our sum of Gaussians in the frequency domain.

Our sum \(\sum G(x, k \pi)\) adds shifted versions of a function together; thus it is a convolution. In this case, it is a convolution of the Gaussian \(G(x, 0)\) with the Dirac comb \(W(\pi x)\). To compute the Fourier transform \(t\): \[\mathcal{F} \{ \sum_{k \in \mathbb{Z}} G(x, k \pi) \} = \mathcal{F} \{ G(x, 0) * W(\pi x) \},\] we can use the fact that the Fourier transform of a convolution is the product of Fourier transforms: \[t = \mathcal{F} \{ G(x, 0) \} \mathcal{F}\{W(\pi x)\}.\]

However, both of these component Fourier transforms are easy to compute; both the Gaussian \(G\) and the Dirac comb \(W\) are fixed points of the Fourier transform. So we have \[t = G(x, 0) W(2 x).\]

\(t\) is the sum of delta functions, scaled appropriately. There is one delta function at 0, which represents the constant offset; one at \(\pm 2\), which represent a sinusoid with wavelength \(\pi\); and then delta functions at every other even integer.

Those delta functions at other even integers are errors in our representation of the sum of Gaussians as a sinusoid. But did we not notice their effect?

The answer is only that those spikes are scaled by very small numbers. How small? Well, \(G(x, 0)\) falls off very rapidly, so even the spikes at \(\pm4\) are scaled by \(G(4, 0) = e^{-4^2 / 2} = e^{-8}\). Thus our error term is many orders of magnitude smaller than the “signal”, which is why we didn’t notice it before.

If you carry out the calculation with a bit more rigor (I dropped a multiple of \(\pi\) or two), you can easily find the exact size of the error term. If you add a second sinusoid to the expression, the error term is smaller than single-float precision, so your computer can even more confidently swear to the “truth” of the theorem.