The KDE Procedure

Fast Fourier Transform

As shown in the last subsection, kernel density estimates can be expressed as a submatrix of a certain convolution. The fast Fourier transform (FFT) is a computationally effective method for computing such convolutions. For a reference on this material, see Press et al. (1988).

The discrete Fourier transform of a complex vector $\mb{z}=(z_{0},\ldots ,z_{N-1})$ is the vector $\mb{Z}=(Z_{0},\ldots ,Z_{N-1})$, where

\[ Z_{j}=\sum _{l=0}^{N-1}z_{l}e^{2\pi i l j/N}, \quad j=0,\ldots ,N-1 \]

and i is the square root of –1. The vector $\mb{z}$ can be recovered from $\mb{Z}$ by applying the inverse discrete Fourier transform formula

\[ z_{l}=N^{-1}\sum _{j=0}^{N-1}Z_{j}e^{-2\pi i l j/N}, \quad l=0,\ldots ,N-1 \]

Discrete Fourier transforms and their inverses can be computed quickly using the FFT algorithm, especially when N is highly composite; that is, it can be decomposed into many factors, such as a power of 2. By the discrete convolution theorem, the convolution of two vectors is the inverse Fourier transform of the element-by-element product of their Fourier transforms. This, however, requires certain periodicity assumptions, which explains why the vectors $\bK $ and $\bC $ require zero-padding. This is to avoid wrap-around effects (see Press et al.; 1988, pp. 410–411). The vector $\bK $ is actually mirror-imaged so that the convolution of $\bC $ and $\bK $ will be the vector of binned estimates. Thus, if S denotes the inverse Fourier transform of the element-by-element product of the Fourier transforms of $\bK $ and $\bC $, then the first g elements of S are the estimates.

The bivariate Fourier transform of an $N_{1}\times N_{2}$ complex matrix having $ (l_{1}+1,l_{2}+1)$ entry equal to $z_{l_{1}l_{2}}$ is the $N_{1}\times N_{2}$ matrix with $ (j_{1}+1,j_{2}+1)$ entry given by

\[ Z_{j_{1}j_{2}}=\sum _{l_{1}=0}^{N_{1}-1} \sum _{l_{2}=0}^{N_{2}-1} z_{l_{1}l_{2}} e^{2\pi i(l_{1}j_{1}/N_{1}+l_{2}j_{2}/N_{2})} \]

and the formula of the inverse is

\[ z_{l_{1}l_{2}}=(N_{1}N_{2})^{-1} \sum _{j_{1}=0}^{N_{1}-1} \sum _{j_{2}=0}^{N_{2}-1} Z_{j_{1}j_{2}} e^{-2\pi i(l_{1}j_{1}/N_{1}+l_{2}j_{2}/N_{2})} \]

The same discrete convolution theorem applies, and zero-padding is needed for matrices $\bC $ and $\bK $. In the case of $\bK $, the matrix is mirror-imaged twice. Thus, if S denotes the inverse Fourier transform of the element-by-element product of the Fourier transforms of $\bK $ and $\bC $, then the upper-left $g_{X}\times g_{Y}$ corner of S contains the estimates.