A New Version of the Fast Gauss Transform

Leslie Greengard and Xiaobai Sun

Abstract. The evaluation of the sum of N Gaussians at M points in space arises as a computational task in diffusion, fluid dynamics, finance, and, more generally, in mollification. The work required for direct evaluation grows like the product NM, rendering large-scale calculations impractical. We present an improved version of the fast Gauss transform [L. Greengard and J. Strain, SIAM J. Sci. Stat. Comput. 12, 79 (1991)], which evaluates the sum of N Gaussians at M arbitrarily distributed points in O(N +M) work, where the constant of proportionality depends only on the precision required. The new scheme is based on a diagonal form for translating Hermite expansions and is significantly faster than previous versions.

1991 Mathematics Subject Classification: 65R10 , 44A35, 35K05 Keywords and Phrases: diffusion, fast algorithms, Gauss transform

1 Introduction

Many problems in mathematics and its applications involve the Gauss transform

G_{δ}f(x) = (πδ)^{−d/2} ∫_{Γ} e^{−|x−y|^{2}/δ}f(y) dy  (δ > 0)  (1)

of a function f, where Γ is some subset of R^{d}. This is, of course, the exact solution to the Cauchy problem

u_{t}(x, t) = Δu(x, t),  t > 0
u(x, 0) = f(x),  x ∈ R^{d}

at time t = δ/4 and corresponds to a mollification of the function f. Similar transforms occur in solving initial/boundary value problems for the heat equation by means of potential theory [3, 10, 11] and in nonparametric statistics [4, 17].

In the present paper, we will focus our attention on the discrete Gauss transform

G(x) = ∑_{j=1}^{N} q_{j} e^{−|x−s_{j}|^{2}/δ} ,  (2)

where the coefficients q_{j} and “source” locations s_{j} are given, and we wish to evaluate the expression (2) at a large number of “target” points x_{j}.

If the number of target points is denoted by M, we can define the rectangular transform matrix G by the formula

G_{ij} = e^{−|x_{i}−s_{j}|^{2}/δ}.  (3)

Direct application of this matrix to the vector q = (q_{1}, …, q_{N})^{T}. requires O(NM) work, which makes large scale calculations prohibitively expensive.

To overcome this obstacle, Greengard and Strain developed a fast Gauss transform [9], which requires only O(N + M) work, with a constant prefactor which depends on the physical dimension d and the desired precision. The amount of memory required is also proportional to N + M, so that the algorithm is asymptotically optimal in terms of both work and storage. In this scheme, the sources and targets can be placed anywhere; methods based on the fast Fourier transform (FFT), by contrast, are restricted to a regular grid and require O(N log N) operations. For the case where the variance δ is not constant:

G(x) = ∑_{j=1}^{N} q_{j} e^{−|x−s_{j}|^{2}/δ_{j}},  (4)

a generalization of the fast Gauss transform has been developed by Strain [18], but we will consider only the simpler case (2) here.

The fast Gauss transform is an analysis-based fast algorithm. Like the closely related fast multipole methods for the Laplace and Helmholtz equations [1, 5, 7, 13, 8, 14, 15, 16], it achieves a speedup in computation by using approximation theory to attain a specified, albeit arbitrarily high, precision. The FFT, on the other hand, is exact in exact arithmetic. It is an algebra-based fast algorithm which uses symmetry properties to reduce the computational work.

2 The original fast Gauss transform

The starting point for the fast algorithms of [9, 18] is the generating function for Hermite polynomials [2, 12]

e^{2xs−s^{2}} = ∑_{n=0}^{∞} \frac{s^{n}}{n!} H_{n}(x) ,

where

H_{n}(x) = (−1)^{n}e^{x^{2}}D^{n}e^{−x^{2}} x ∈ R

and D = d/dx. A small amount of algebra leads to the expansion

e^{−(x−s)^{2}/δ} = ∑_{n=0}^{∞} \frac{1}{n!} (\frac{s − s_{0}}{\sqrt{δ}} )^{n} h_{n} (\frac{x − s_{0}}{\sqrt{δ}} ) ,

where

h_{n}(x) = (−1)^{n}D^{n}e^{−x^{2}}

and s_{0} is an arbitrary point.

This formula describes the Gaussian field e^{−(x−s)^{2}/δ} at the target x due to the source at s as an Hermite expansion centered at s_{0}. The higher dimensional analog of (5) is obtained using multi-index notation. Let x and s lie in d-dimensional Euclidean space R^{d}, and consider the Gaussian

e^{−|x−s|^{2}} = e^{−(x_{1}−s_{1})^{2}−...−(x_{d}−s_{d})^{2}}.

For any multi-index α = (α_{1}, α_{2}, …, α_{d}) and any x ∈ R^{d}, we define

|α| = α_{1} + α_{2} + … + α_{d}
α! = α_{1}!α_{2}! … α_{d}!
x^{α} = x_{1}^{α_{1}}x_{2}^{α_{2}} … x_{d}^{α_{d}}
D^{α} = ∂_{1}^{α_{1}}∂_{2}^{α_{2}} … ∂_{d}^{α_{d}}

where ∂_{i} is differentiation with respect to the ith coordinate in R^{d}. If p is an integer, we say α ≥ p if α_{i} ≥ p for 1 ≤ i ≤ d.

The multidimensional Hermite polynomials and Hermite functions are defined by

H_{α}(x) = H_{α_{1}}(x_{1}) … H_{α_{d}}(x_{d})
h_{α}(x) = e^{−|x|^{2}}H_{α}(x) = h_{α_{1}}(x_{1}) … h_{α_{d}}(x_{d})  (5)

where |x|^{2} = x_{1}^{2} + … + x_{d}^{2}. The Hermite expansion of a Gaussian in R^{d} is then simply

e^{−|x−s|^{2}} = ∑_{α≥0} \frac{(x − s_{0})^{α}}{α!} h_{α}(s − s_{0}) .  (6)

Lemma 2.1 ([9], 1991) Let N_{B} sources s_{j} lie in a box B with center s_{B} and side length \sqrt{δ}. Then the Gaussian field due to the sources in B,

G(x) = ∑_{j=1}^{N_{B}} q_{j} e^{−|x−s_{j}|^{2}/δ},  (7)

is equal to a single Hermite expansion about s_{B}:

G(x) = ∑_{α≥0} A_{α} h_{α} (\frac{x − s_{B}}{\sqrt{δ}} ) .

The coefficients A_{α} are given by

A_{α} = \frac{1}{α!} ∑_{j=1}^{N_{B}} q_{j} (\frac{s_{j} − s_{B}}{\sqrt{δ}} )^{α} .  (8)

The error E_{H}(p) due to truncating the series after p^{d} terms satisfies the bound:

|E_{H}(p)| = | ∑_{α≥p} A_{α}h_{α} (\frac{x − s_{B}}{\sqrt{δ}} ) | ≤ 2.75^{d} Q_{B} (\frac{1}{p!})^{d/2} (\frac{1}{2})^{(p+1)d/2} (9)

where Q_{B} = ∑ |q_{j}|.

Lemma 2.2 ([9], 1991) Let N_{B} sources s_{j} lie in a box B with center s_{B} and side length \sqrt{δ} and let x be a target point in a box C with center x_{C}. Then the corresponding Hermite expansion

G(x) = ∑_{α≥0} A_{α} h_{α} (\frac{x − s_{B}}{\sqrt{δ}} ) .

can be expanded as a Taylor series of the form

G(x) = ∑_{β≥0} B_{β} (\frac{x − x_{C}}{\sqrt{δ}} )^{β} .

The coefficients B_{β} are given by

B_{β} = \frac{(−1)^{|β|}}{β!} ∑_{α≥0} A_{α} h_{α+β} (\frac{s_{B} − x_{C}}{\sqrt{δ}} ) .  (10)

The error E_{T} (p) due to truncating the series after p^{d} terms satisfies the bound:

|E_{T} (p)| = | ∑_{β≥p} B_{β} (\frac{x − x_{C}}{\sqrt{δ}} )^{β} | ≤ 2.75^{d} Q_{B} (\frac{1}{p!})^{d/2} (\frac{1}{2})^{(p+1)d/2}  (11)

These are the only tools required to construct a simple fast algorithm for the evaluation of

G(x_{i}) = ∑_{j=1}^{N} q_{j}e^{−|x_{i}−s_{j}|^{2}/δ}  (12)

for 1 ≤ i ≤ M, using O(M + N) work. By shifting the origin and rescaling δ if necessary, we can assume (as a convenient normalization) that the sources s_{j} and targets x_{i} all lie in the unit box B_{0} = [0, 1]^{d}.

Algorithm

Step 1. Subdivide B_{0} into smaller boxes with sides of length \sqrt{δ} parallel to the axes. Assign each source s_{j} to the box B in which it lies and each target x_{i} to the box C in which it lies. The source boxes B and the target boxes C may, of course, be the same.

Step 2. Given ϵ, use Lemma 2.1 to create an Hermite expansion for each source box B with p^{d} terms satisfying:

G(x) = ∑_{B} ∑_{s_{j}∈B} q_{j}e^{−|x−s_{j}|^{2}/δ}

= ∑_{B} ∑_{α≤p} A_{α}(B)h_{α} (\frac{x − s_{B}}{\sqrt{δ}} ) + O(ϵ)

where

A_{α}(B) = \frac{1}{α!} ∑_{s_{j}∈B} q_{j} (\frac{s_{j} − s_{B}}{\sqrt{δ}} )^{α} .  (13)

The amount of work required for this step is of the order p^{d}N.

Consider now a fixed target box C. For each x_{j} ∈ C, we need to evaluate the total field due to sources in all boxes of type B. Because of the exponential decay of the Gaussian field, however, it is easy to verify that, if we include only the sources in the nearest (2r + 1)^{d} boxes, we incur an error bounded by Qe^{−r^{2}}, where Q = ∑_{j=1}^{N} |q_{j}|. Given a desired precision ϵ, we can always choose r so that this truncation error is bounded by Qϵ. With r = 4, for example, we get single precision accuracy (ϵ = 10^{−7}) and with r = 6, we get double precision (ϵ = 10^{−14}). We denote the nearest (2r +1)^{d} boxes as the interaction region for box C, denoted by IR(C).

Step 3. For each target box C, use Lemma 2.2 to transform all Hermite expansions in source boxes within the interaction region into a single Taylor expansion. Thus, we approximate G(x) in C by

G(x) = ∑_{B} ∑_{s_{j}∈B} q_{j}e^{−|x−s_{j}|^{2}/δ}

= ∑_{β≤p} C_{β} (\frac{x − x_{C}}{\sqrt{δ}} )^{β} + O(ϵ)

where

C_{β} = \frac{(−1)^{|β|}}{β!} ∑_{B∈IR(C)} ∑_{α≤p} A_{α}(B) h_{α+β} (\frac{s_{B} − x_{C}}{\sqrt{δ}} ) ,  (14)

and the coefficients A_{α}(B) are given by (13). Because of the product form (5) of h_{α+β}, the computation of the p^{d} coefficients C_{β} involves only O( d p^{d+1}) operations for each box B. Therefore, a total of O((2r + 1)^{d} d p^{d+1}) work per target box C is required. Finally, evaluating the appropriate Taylor series for each target x_{i} requires O(p^{d}M) work. Hence this algorithm has net CPU requirements of the order

O((2r + 1)^{d} d p^{d+1}N_{box}) + O(p^{d}N) + O(p^{d}M) ,

where the number of boxes N_{box} is bounded by min(δ^{−d/2}, N + M). The work is cleanly decoupled into three parts; O(p^{d}N) to form Hermite expansions, O(p^{d}M) to evaluate Taylor series, and a constant term depending on the number of box-box interactions and the cost of transforming Hermite expansions into Taylor series.

Remark: A proper implementation of the fast Gauss transform is a bit more complex. For example, if a box contains only a few sources, it is more efficient to compute their influence directly than to use expansions.

Suppose now that the source boxes are denoted by B_{1}, B_{2}, …, B_{S}, that the target boxes are denoted by C_{1}, C_{2}, …, C_{T}, that N_{j} sources lie in box B_{j}, that M_{j} targets lie in box C_{j}, and that the points are ordered so that

{s_{1}, …, s_{N_{1}}} ⊂ B_{1}
{s_{N_{1}+1}, …, s_{N_{1}+N_{2}}} ⊂ B_{2}
…
{s_{N−N_{S}+1}, …, s_{N}} ⊂ B_{S},
{x_{1}, …, x_{M_{1}}} ⊂ C_{1}
{x_{M_{1}+1}, …, x_{M_{1}+M_{2}}} ⊂ C_{2}
…
{x_{M−M_{T} +1}, …, x_{M}} ⊂ C_{T} .

Then the approximation G_{ϵ} to the discrete Gauss transform matrix (3) can be written in the factored form

G_{ϵ} = D · E · F.  (15)

Here, F is a block diagonal matrix of dimension S × S. The jth diagonal block F(j) ∈ R^{p^{d}×N_{j}} satisfies

F(j)_{n,m} = \frac{1}{α_{n}!} (\frac{s_{m} − s_{B_{j}}}{\sqrt{δ}} )^{α_{n}} ,

where s_{B_{j}} is the center of box B_{j} and the p^{d} Hermite expansion coefficients are ordered in some fashion from n = 0, …, p^{d}. D is very similar. It is a block diagonal matrix of dimension T × T, with the jth diagonal block D(j) ∈ R^{M_{j}×p^{d}} satisfying

D(j)_{n,m} = \frac{(−1)^{β_{m}}}{β_{m}!} (\frac{x_{n} − x_{C_{j}}}{\sqrt{δ}} )^{β_{m}}

where x_{C_{j}} is the center of box C_{j}, and the p^{d} Taylor expansion coefficients are ordered from n = 0, …, p^{d} in the same fashion as the Hermite series. Note that, if the sources and targets coincide, then D is the transpose of F.

The mapping E is a sparse block matrix of dimension T × S, with up to (2r + 1)^{d} nonzero entries per row. The nonzero entries E(ij) are matrices of dimension p^{d} × p^{d}, corresponding to a conversion of the Hermite series for box S_{j} into a Taylor series for box T_{i}, assuming S_{j} is in the interaction region IR(T_{i}). The matrix entries are dense.

E(ij)_{nm} = h_{α_{n}+β_{m}} (\frac{s_{B_{j}} − x_{C_{i}}}{\sqrt{δ}} ) .

Given this notation, Step 2 of the fast Gauss transform described above corresponds to multiplying the vector {q_{1}, q_{2}, …, q_{N}} by F. Step 3 of the fast Gauss transform corresponds to multiplying the output of Step 2 by E to create all the Taylor expansions. The result is then multiplied by D to evaluate the Taylor series at all target locations.

Remark: The factorization (15) reveals the structure of G_{ϵ}. When δ is large enough, only one box is created and the rank of G_{ϵ} is bounded by p^{d} (the order of the factor E). When δ is very small, the dimensions of E grow, but it becomes sparse and structured.

3 Diagonal form for translation operators

Our new version of the fast Gauss transform is based on replacing Hermite and Taylor expansions with an expansion in terms of exponentials (plane waves). The starting point is the Fourier relation

e^{−|x−s|^{2}/δ} = (\frac{1}{2\sqrt{π}} )^{d} ∫_{R^{d}} e^{−|k|^{2}/4}e^{ik·(x−s)/\sqrt{δ}} dk  (16)

which is easily seen to satisfy the estimate

| e^{− \frac{|x−s|^{2}}{δ}} − (\frac{1}{2\sqrt{π}} )^{d} ∫_{|k|≤K} e^{− \frac{|k|^{2}}{4}} e^{i \frac{k·(x−s)}{\sqrt{δ}}} dk | ≤ \{ \begin{matrix}e^{− \frac{K^{2}}{4}} & for d = 1, 2 \\ Ke^{− \frac{K^{2}}{4}} & for d = 3. \end{matrix}

Setting K = 7.5, the truncation error from ignoring high frequency contributions is approximately 10^{−7}. Ssetting K = 12, the truncation error is approximately 10^{−14}. It still remains to discretize the Fourier integral in (16) within the range determined by K. The trapezoidal rule is particularly appropriate here since it is rapidly convergent for functions which have decayed at the boundary. Note, however, that the integrand is more and more oscillatory as x − s grows. Fortunately, we only need accurate quadrature when s is within the interaction region of x, so that |x − s|/\sqrt{δ} ≤ 5 for seven digit precision and |x − s|/\sqrt{δ} ≤ 7 for fourteen digit precision. It is easy to verify that p = 12 equispaced modes in the interval [0, 7.5] are sufficient to reduce the quadrature error to 10^{−7} when |x − s|/\sqrt{δ} ≤ 5 and that p = 24 equispaced modes in the interval [0, 12] are sufficient to reduce the quadrature error to 10^{−14} when |x − s|/\sqrt{δ} ≤ 7.

Thus, for a source box B with center s_{B}, we replace the Hermite series of Lemma 2.1 with

G(x) = ∑_{s_{j}∈B} q_{j}e^{−|x−s_{j}|^{2}/δ}

= ∑_{β≤p} C_{β}e^{i \frac{Kβ·(x−s_{B})}{p\sqrt{δ}}} + O(ϵ),

where

C_{β} = (\frac{K}{2p\sqrt{π}} )^{d} e^{−|β|^{2}K^{2}/(4p^{2})} ∑_{j=1}^{N_{B}} q_{j} e^{−i \frac{Kβ·(s_{j} −s_{B})}{p\sqrt{δ}}} .

There are two reasons to prefer this form. First, the translation operator described in Lemma 2.2 becomes diagonal.

Corollary 3.1 Let N_{B} sources s_{j} lie in a box B with center s_{B} and side length \sqrt{δ} and let x be a target point in a box C with center x_{C}. Then the plane wave expansion

G(x) = ∑_{β≤p} C_{β}e^{i \frac{Kβ·(x_{i}−s_{B})}{p\sqrt{δ}}} + O(ϵ),

can be expanded about x_{C} as

G(x) = ∑_{β≤p} D_{β}e^{i \frac{Kβ·(x_{i}−x_{C})}{p\sqrt{δ}}} + O(ϵ).

The coefficients D_{β} are given by

D_{β} = C_{β} e^{i \frac{Kβ·(x_{C} −s_{B})}{p\sqrt{δ}}} .  (17)

In terms of matrix factorization, we have

G_{ϵ} = D′ · E′ · F′.  (18)

In this formulation, the diagonal blocks of F′ and D′ are given by

F′(j)_{n,m} = (\frac{K}{2p\sqrt{π}} )^{d/2} e^{\frac{−|βn|K}{2p}} e^{−i \frac{Kβn·(sm−s_{B_{j}} )}{p\sqrt{δ}}} .

D′(j)_{n,m} = (\frac{K}{2p\sqrt{π}} )^{d/2} e^{\frac{−|βn|K}{2p}} e^{i \frac{Kβm·(xn−x_{C_{j}} )}{p\sqrt{δ}}} .

As in the original algorithm, note that if the sources and targets coincide, then D′ is the adjoint of F′. The nonzero entries E′(ij) are now diagonal matrices of dimension p^{d} × p^{d}, with entries defined in (17). The net cost of all translations per target box is reduced from O((2r + 1)^{d} d p^{d+1}) work to O((2r + 1)^{d} p^{d}) work.

The second (and more important) reason to prefer the new form is that the number of translations can be dramatically reduced. We describe the modification to the algorithm in the one-dimensional case. For this, imagine that we are sweeping across all boxes from left to right and that, at present, a target box C_{j} has accumulated all plane wave expansions from source boxes within its interaction region (Fig. 1(a)). The net expansion can be shifted to the center of C_{j+1} using Corollary 3.1. By adding in the contribution from the box marked by + and subtracting the contribution from the box marked by −, we have the correct plane wave expansion for box C_{j+1} (Fig. 1(b)). Thus, (2r + 1) translations are replaced by three. In d-dimensions, the cost O((2r + 1)^{d} p^{d}) work can be reduced to O(3d p^{d}), by sweeping across each dimension separately.

4 Conclusions

We have presented a new version of the fast Gauss transform, which uses plane wave expansions to diagonalize the translation of information between boxes. The

C^{j}C^{j+1}

Figure 1: After shifting the expansion from box C_{j} to box C_{j+1}, one needs only to subtract the contribution from the box marked − and add the contribution from the box marked +. (The interaction regions are indicated by the square brackets).

approach is similar to the new diagonal forms used in fast multipole methods for the Laplace and Helmholtz equations [8, 6, 13]. When the present improvements have been incorporated into existing fast Gauss transform codes, the resulting scheme should provide a powerful kernel for one, two and three-dimensional calculations.

References

[1] J. Carrier, L. Greengard, and V. Rokhlin, A Fast Adaptive Multipole Algorithm for Particle Simulations, Siam J. Sci. Stat. Comput., 9 (1988), pp. 669–686.

[2] H. Dym and H. P. McKean, Fourier Series and Integrals, Academic Press, San Diego, 1972.

[3] A. Friedman, Partial Differential Equations of Parabolic Type, Prentice-Hall, New Jersey, 1964.

[4] S. Geman and C. Hwang, Nonparametric Maximum Likelihood Estimation by the Method of Sieves, Ann. Statist., 10 (1982), pp. 401–414.

[5] L. Greengard, Fast algorithms for classical physics, Science, 265 (1994), pp. 909-xxx.

[6] L. Greengard, J. Huang, V. Rokhlin and S. Wandzura, Accelerating fast multipole methods for low frequency scattering, CMCL Report 1998-003, New York University.

[7] L. Greengard and V. Rokhlin, A fast algorithm for particle simulations, J. Comput. Phys., 73 (1987), pp. 325–348.

[8] L. Greengard and V. Rokhlin, A new version of the fast multipole method for the Laplace equation in three dimensions, Acta Numerica, 6 (1987), pp. 229–269.

[9] L. Greengard and J. Strain, The Fast Gauss Transform, SIAM J. Sci. Stat. Comput., 12 (1991), pp. 79–94.

[10] L. Greengard and J. Strain, A Fast Algorithm for the Evaluation of Heat Potentials, Comm. on Pure and Appl. Math., 43 (1990), pp. 949–963.

[11] R. B. Guenther and J. W. Lee, Partial Differential Equations of Mathematical Physics and Integral Equations, Prentice-Hall, Englewood Cliffs, New Jersey, 1988.

[12] E. Hille, A Class of Reciprocal Functions, Ann. Math., 27 (1926), pp. 427–464.

[13] T. Hrycak and V. Rokhlin, An improved fast multipole algorithm for potential fields, Department of Computer Science Research Report 1089, Yale University, 1995.

[14] V. Rokhlin, Rapid solution of integral equations of classical potential theory, J. Comput. Phys. 60 (1985), pp. 187–207.

[15] V. Rokhlin, Rapid solution of integral equations of scattering theory in two dimensions, J. Comput. Phys., 86 (1990), pp. 414–439.

[16] V. Rokhlin, Diagonal forms of translation operators for the Helmholtz equation in three dimensions, Appl. and Comput. Harmonic Analysis 1 (1993), pp. 82–93.

[17] B. W. Silverman, Density Estimation for Statistics and Data Analysis, Chapman and Hall, London, 1986.

[18] J. Strain, The Fast Gauss Transform with variable scales, SIAM J. Sci. Stat. Comput. 12 (1991), pp. 1131–1139.

Leslie Greengard Courant Institute New York University New York, NY 10012-1110 USA greengard@cims.nyu.edu Xiaobai Sun Department of Computer Science Duke University Durham, NC 27708-0129 USA xiaobai@cs.duke.edu