FFT Convolution for Equally Spaced Grids

Why convolution helps on a grid

At each look in a group sequential calculation, the survivor density is propagated through a Gaussian transition and then restricted to states that have not crossed a stopping boundary. On an equally spaced grid, the transition depends on the difference between the destination and source states. That translation structure makes the update a convolution: instead of evaluating every source-to-destination contribution separately, an FFT can evaluate the full linear convolution efficiently.

The speed-up is most useful when the grid is reasonably long and the transition kernel is dense. It does not remove the need to approximate a continuous density on a grid, represent the stopping boundary, or control probability outside a truncated computational range. Those are separate sources of numerical error.

Model and discrete update

For clarity, consider a state process in information time with Gaussian increments. Conditional on the current state \(X_i=x\), suppose the next state satisfies

X_{i+1}\mid X_i=x \sim N(x+\theta\Delta t_i,\Delta t_i),

where \(\Delta t_i>0\) is the information-time increment and \(\theta\) is the drift under the probability model being evaluated. This transition assumes independent Gaussian increments, with variance equal to the information-time increment. The null calculation uses the null drift; a power calculation uses the drift specified under the alternative.

Let \(q_i(x)\) be the survivor density at look \(i\): it describes probability density among paths that have not crossed a stopping boundary before that look. The unrestricted transition density from \(x\) to \(y\) is

K_i(y-x)=\frac{1}{\sqrt{\Delta t_i}}\, \phi\!\left(\frac{y-x-\theta\Delta t_i}{\sqrt{\Delta t_i}}\right),

where \(\phi\) is the standard normal density. Before applying the next look's stopping rule, the propagated density is

\widetilde q_{i+1}(y)=\int q_i(x)K_i(y-x)\,dx.

The survivor update then retains \(\widetilde q_{i+1}(y)\) only at states that survive the boundary rule. This tutorial concerns accelerating the transition integral; boundary handling remains part of the sequential calculation.

On a grid with spacing \(h\), write \(x_j=x_0+jh\), and let \(q_j\) approximate \(q_i(x_j)\). For destination \(y_k\), a Riemann-sum approximation is

\widetilde q_{k}\approx h\sum_j q_j K_i(y_k-x_j).

If the source and destination grids have the same spacing and origin, \(y_k-x_j=(k-j)h\). The kernel is therefore a function of the index difference \(k-j\), which is exactly the structure required for discrete convolution.

Definition: A discrete convolution of arrays \(a_j\) and \(b_r\) is the sequence \(c_s=\sum_j a_j b_{s-j}\). For the grid update, take \(a_j=q_j\) and \(b_r=K_i(rh)\), with the drift included in the kernel argument: \(b_r=K_i(rh)\) as defined above. The propagated density at destination index \(k\) is \(h c_k\), subject to the chosen finite ranges and indexing.

From the transition sum to an FFT

The FFT computes a discrete Fourier transform efficiently. The key identity is that multiplication of Fourier transforms corresponds to convolution of the original arrays. One must distinguish linear convolution, which sums only actual source-kernel pairs, from the circular convolution produced by a discrete Fourier transform of fixed length.

Suppose the source array has \(N\) entries, indexed \(0,\ldots,N-1\), and the retained kernel has \(2M+1\) entries for lags \(-M,\ldots,M\). Store the kernel in an array \(b\) using \(b_{r+M}=K_i(rh)\). The full linear convolution has length \(N+2M\). Its entry at index \(k+M\) supplies the transition sum for destination index \(k\):

c_{k+M}=\sum_{j=0}^{N-1}q_j b_{k-j+M}, \qquad \widetilde q_k\approx h c_{k+M}.

For an FFT of length \(P\), append zeros to both arrays so that

P\ge N+(2M+1)-1=N+2M.

Then multiply their length-\(P\) FFTs and apply the inverse FFT. The required linear-convolution values occupy the first \(N+2M\) positions of the result; the index shift above selects the desired destination values. In practice, a convenient FFT length is often the next power of two meeting the lower bound, though other efficient lengths are also possible.

Derivation and the role of padding

The finite convolution has terms \(q_j b_\ell\) only when both \(j\) and \(\ell\) are valid array indices. For the linear convolution, its output index \(s\) is the sum \(j+\ell\). The smallest possible sum is \(0\), and the largest is \((N-1)+((2M+1)-1)=N+2M-1\); hence there are \(N+2M\) output positions.

An FFT of length \(P\), however, treats each input as periodic with period \(P\). Without enough padding, terms whose linear-convolution indices exceed \(P-1\) wrap around and are added to early output indices. If \(P\) is at least the full linear-convolution length, no nonzero output wraps into the retained result. This is why zero-padding is part of the calculation, not merely a speed preference.

The Gaussian kernel itself is not finite. Keeping lags only from \(-M\) to \(M\) is an additional approximation: it discards transitions beyond the retained lag range. Padding prevents circular aliasing of the retained kernel; it does not recover the Gaussian tail that was omitted.

Worked Example: A Small Linear Convolution

Consider a deliberately small grid calculation with \(h=1\), source values \(q=[0.2,0.5,0.3]\), and a three-point kernel \(b=[0.2420,0.3989,0.2420]\), stored for lags \(-1,0,1\). These kernel values resemble a unit-variance Gaussian density at those lags; the short kernel is used here only to make the indexing arithmetic visible.

The source length is \(N=3\), and the kernel length is \(3\), so the full linear convolution has length \(3+3-1=5\). Its entries are

[0.2,0.5,0.3]*[0.2420,0.3989,0.2420] =[0.04840,\ 0.20078,\ 0.32045,\ 0.24067,\ 0.07260].

For example, the middle entry is \(0.2(0.2420)+0.5(0.3989)+0.3(0.2420)=0.32045\). The destination indices \(k=0,1,2\) correspond to convolution entries \(k+M=1,2,3\), since \(M=1\). Multiplying by \(h=1\), the three propagated density approximations are therefore \(0.20078,\ 0.32045,\ 0.24067\). An FFT of length at least 5 produces this linear convolution after zero-padding. Because this example retains only three kernel values, it should not be treated as a well-resolved Gaussian update.

Grid padding and boundary truncation

Two different forms of truncation matter in a practical implementation. First, the computational state space may be restricted to a finite interval, as in the earlier material on truncating the Gaussian state space. Probability that moves outside that interval is not represented unless the calculation explicitly accounts for it. Second, the transition kernel may be cut off at a finite lag to make its array manageable. That removes the tails of the Gaussian transition even if the state-space interval is wider.

Choose the state-space range wide enough that the omitted state probability is negligible for the desired accuracy, and choose the kernel lag range wide enough that its omitted transition probability is negligible. These choices are related but not interchangeable. Also, after propagating the density, apply the specified stopping-boundary rule on the destination grid. Values beyond a stopping boundary represent stopped paths, not merely FFT wraparound.

Worked Example: Selecting a Gaussian Kernel Range

Suppose a transition has variance \(\Delta t_i=0.25\), so its standard deviation is \(0.5\), and its mean increment is \(\theta\Delta t_i=0.05\). Let the grid spacing be \(h=0.05\), and retain lags \(r\) for which \(|rh-0.05|\le 2.0\).

The retained interval is four standard deviations on either side of the transition mean, because \(2.0/0.5=4\). Under the continuous Gaussian transition, the probability outside that interval is

P(|X_{i+1}-X_i-0.05|>2.0) =2\Phi(-4)\approx 0.00006334.

The lag condition gives \(-39\le r\le41\): the endpoints correspond to \(rh=-1.95\) and \(2.05\), whose differences from \(0.05\) are \(-2.00\) and \(2.00\). Thus the stored kernel has \(81\) entries. The continuous tail calculation indicates the scale of the omitted probability; the actual discrete quadrature and finite state grid should still be checked, because their errors are not identical to this continuous tail probability.

Worked Example: Why an Unpadded FFT Is Wrong

Take \(a=[1,2,3]\) and \(b=[4,5,6]\). Their linear convolution is \([4,13,28,27,18]\), with length \(3+3-1=5\). If both arrays are transformed at length \(3\), the result is a circular convolution. The last two linear terms wrap around and are added to the first two positions:

[4+27,\ 13+18,\ 28]=[31,\ 31,\ 28].

Those values are not the desired linear convolution. Padding each input to length at least 5 (for example, to length 8), multiplying the FFTs, and taking the inverse FFT recovers \([4,13,28,27,18]\). In a density update, wraparound can look like probability arriving from the opposite edge of the grid, so it may be mistaken for a boundary or state-space effect.

Practical checks and common failures

An FFT update is reliable only when its indexing and approximation choices match the transition calculation. Useful checks include comparing a few destinations against the direct sum, confirming that the kernel is centered at the correct drift increment, and verifying that the post-update density is nonnegative up to small floating-point error. Keep units straight: \(K_i\) has units of inverse state, \(q_j\) is a density, and the factor \(h\) converts the sum into a density approximation.

Insufficient zero-padding. If the FFT length is shorter than the full linear-convolution length, circular wraparound contaminates the result. Check \(P\ge N+2M\).
Incorrect kernel shift. A drift increment that is omitted or given the wrong sign shifts the transition. Check that the kernel is evaluated at \(rh-\theta\Delta t_i\) in its normal-density formula.
Missing grid-width factor. The discrete sum approximates an integral, so include \(h\). Omitting it changes the scale of the density and its integrated mass.
Confusing the two truncations. Zero-padding prevents circular aliasing; it does not prevent loss beyond the finite state interval or restore discarded Gaussian kernel tails.
Overinterpreting tiny negative values. FFT roundoff can produce very small negative numbers. First compare their size with the density scale and check the direct calculation. Clipping substantial negatives without investigation can hide an indexing error.

A useful diagnostic is to compare selected FFT outputs with the direct vectorized sum \(h\sum_j q_jK_i(y_k-x_j)\), using exactly the same source values, kernel range, and destination indices. Agreement checks the FFT implementation, but does not validate the grid spacing, state-space range, or boundary representation. Those approximations require their own sensitivity checks.

When to use FFT convolution

Direct vectorized summation is often preferable for small grids, a sparse transition kernel, or a calculation where only a few destination values are needed. It is also easy to inspect term by term. For a dense kernel on a grid with \(N\) points, direct summation typically requires a number of operations proportional to \(N^2\); an FFT-based convolution uses a padded length \(P\) and requires work proportional to \(P\log P\). The FFT can therefore be faster on larger grids, although the crossover depends on implementation and hardware.

The FFT is not automatically more accurate. Both approaches use the same grid approximation if given the same discretized inputs, and the FFT additionally requires careful padding and indexing. Unequal source and destination spacings, a kernel that varies with location, or a non-Gaussian transition may remove the simple convolution structure. The method is an acceleration for a translation-invariant transition on equally spaced grids, not a replacement for checking that those conditions hold.

Key takeaway: On an equally spaced grid, a Gaussian survivor-density transition is a discrete convolution multiplied by the grid spacing. Use enough zero-padding for the full linear convolution, account separately for omitted Gaussian tails and finite state-space truncation, and apply the stopping rule after propagation.

Check Your Understanding

Answer each question, then compare your reasoning with the solutions below.

A source array has 120 entries and a retained kernel has 41 entries. What is the minimum FFT length that guarantees the full linear convolution without wraparound?
Why does zero-padding not eliminate error from truncating a Gaussian kernel?
For a grid spacing of \(0.02\), what factor converts a discrete sum of density times transition-density values into a Riemann approximation to the transition integral?
A density update has small negative FFT outputs. Give one appropriate check before deciding how to handle them.

Solutions. (1) The full convolution length is \(120+41-1=160\), so \(P\) must be at least 160. (2) Padding prevents nonzero entries of the retained linear convolution from wrapping around; it cannot add kernel values that were discarded before the FFT. (3) Multiply the sum by \(h=0.02\). (4) Compare selected destinations with the direct sum using the same inputs and indexing, and assess whether the negatives are tiny relative to the density scale; substantial discrepancies require investigating the implementation rather than simply clipping values.