Memory-efficient $w$-projection with the fast Gauss transform

The aim of convolutional gridding is to take visibilities sampled on arbitrary $u,v,w$ coordinates and interpolate them onto a regular grid so that a 2D fast Fourier transform (FFT) can be used to estimate the sky brightness distribution (Taylor et al., 1999, p. 134ff). Mathematically, the grid is evaluated at a point $(u_{c},v_{c})$ by multiplying each visibility by a shifted convolution function according to:

For a fixed $w$, equation 2 is the convolution of the anti-aliasing function with a complex Gaussian. If we choose a Gaussian for the anti-aliasing function $A$, the resulting convolution function is also Gaussian. To simplify notation, we will consider the radially symmetric problem by introducing a new parameter $t^{2}=u^{2}+v^{2}$, and can now write the Fresnel term as:

The gridding function is given by the convolution $G(t)*A(t)$, which is another Gaussian, whose variance is equal to the sum of the variances of the original Gaussians, i.e., to within a scaling factor:

The width of the envelope is given by

The convolution function is, therefore, a complex Gaussian, which can be expressed as the product of real-valued Gaussian envelope (the $\delta_{R}$ term in equation 10) with a complex chirp (the $\delta_{I}$ term). Clearly, as $|w|$ increases, the width of the Gaussian envelope $\delta_{R}$ increases.

For a given visibility indexed by $j$ the form of the convolution function can be calculated from its $w$ coordinate using:

Up to this point, the $u,v,w$ co-ordinates and the Gaussian widths have been in units of wavelength. For the remainder of this paper, width will often be quoted in terms of pixels. We use the terms pixel and $uv$-cell interchangeably. To convert from wavelengths to pixels, a size of $uv$-cell is required. The size of the $uv$-cell, (in wavelengths) is set by the inverse of the desired field of view (in radians) i.e. $\phi=1/\theta_{\rm f.o.v}$.

Therefore, a Gaussian width in pixels can be calculated from the width in wavelength by:

In this section we introduce the Fast Gauss Transform (FGT) with variable scales (Strain, 1991) (hereafter S91), which is an approximate technique for calculating the sum of Gaussians. We conform to the language of S91, which uses the term ‘source’ to describe a Gaussian function with a given position, amplitude and width. We will describe the application of the FGT onto the gridding problem in section 3.2.

The FGT works by partitioning the evaluation region into square boxes, with each box containing a set of Taylor coefficients (Fig. 1). For each Gaussian ‘source’, the Taylor coefficients for the boxes surrounding the source position are updated, with the size of the updated region being defined by the ‘width’ of the source, i.e. wider Gaussians update more boxes. When all sources have been applied, the Taylor coefficients are evaluated at any number of target locations. The algorithm is parameterised by three numbers: $r$, which sets the box size, $p$, which defines number of Taylor coefficients to consider, and $\delta$ which sets the smallest size of Gaussian that will be considered. Reducing the box size, or increasing the number of Taylor coefficients increases the accuracy of the sum.

Now that we have a Gaussian convolution function (Equation 14), we can solve the sum of Equation 1 the FGT with source-dependant scales. In the language of gridding visibilities, the evaluation region is the $uv$ plane. The ‘sources’ are the visibilities, with the source position being given by the coordinates of the visibility on the $uv$ plane. The width of the source is related to the $w$ coordinate of the visibility as in Equation 15, with the minimum width when $w=0$, being essentially the width of the anti-aliasing function. When all the visibilities have been applied to the Taylor coefficients, the Taylor coefficients are evaluated on a regular grid, for each $uv$ cell on the $uv$ plane.

For a detailed description of the algorithm, we refer the reader to S91. For the purposes of explaining our implementation, we outline the algorithm in the following text, using similar notation to S91. This description applies to the one-dimensional case, but S91 describes the relationships to extend it to the multi-dimensional case (using multi-index notation), for which the 2D case is relevant here. Additionally, this description is only valid for real Gaussians. The extension to complex Gaussians will be considered in the next section.

The algorithm proceeds as follows:

1. 1.

   Given values of $r$ and $p$ calculate:

   $$\epsilon=\frac{r^{p}_{p}}{1-r_{p}},$$

   (17)

   where

   $\displaystyle r_{p}$

   $\displaystyle=$

   $\displaystyle r\sqrt{\frac{e}{p+1}}<1$

   (18)

2. 2.

   Partition the $uv$ plane into the set of boxes $B$ with centers $t_{B}$, and side length $L=r\sqrt{2\delta}$ parallel to the axes, where $\delta=\min_{j}\delta_{j}=\delta_{A}$ is the smallest Gaussian width to be considered.

3. 3.

   For each visibility indexed by $j$, calculate the range $R_{j}=\sqrt{-\delta_{j}\log\epsilon}$. Also calculate the order $q\leq p$ satisfying

   $$\frac{r^{q}_{q,j}}{1-r_{q,j}}\leq\epsilon$$

   (19)

   where $r_{q,j}=r\sqrt{e\delta/\delta_{j}(q+1)}$.

4. 4.

   For each box within $R_{j}$ of the visibility position $t_{j}$, update the Taylor coefficients according to:

   $$C_{\beta}=\frac{1}{\beta!}\sum_{\delta_{j}\geq b}v_{j}\left(\frac{\delta}{\delta_{j}}\right)^{|\beta|/2}h_{\beta}\left(\frac{t_{j}-t_{B}}{\sqrt{\delta}}\right)$$

   (20)

   where $\beta$ runs from $0$ to $q$, $v_{j}$ is the visibility weight, and $h_{\beta}$ is the Hermite function (see S91).

5. 5.

   Once all the visibilities have been applied, evaluate the Taylor coefficients at each cell in the $uv$ plane, by finding the box in which the $uv$ cell is situated, and calculating the sum:

   $$F(t)=\sum_{\beta\leq p}C_{\beta}\left(\frac{t_{C}-t_{B}}{\sqrt{\delta}}\right)^{\beta}$$

   (21)

   where $t_{C}$ is the position of the $uv$ cell, and $\beta$ runs from $0$ to $p$.

The error in the evaluation of the sum is bounded by:

$r$ and $p$ are chosen to make $|E_{p}|\leq V\epsilon$.

The formulation of S91 is explicitly for real Gaussians, i.e. with $\delta_{j},v_{j}\in\Re$, but in our case $\delta_{j},v_{j}\in\mathds{C}$. A full derivation of the algorithm for complex values is outside the scope of this work . The important questions are whether the series of Equation 20 converges (and for what values of a complex $\delta$) and what is a suitable range ($R_{j}$).

Empirically, we have determined that:

1. 1.

   Using a range of $R_{j}=\sqrt{-\delta_{R,j}\log\epsilon}$ provides good results. This also makes sense because it is the amplitude of the envelope that determines the range.

2. 2.

   For calculating the order $q$ (Equation 19) we used

   $$\delta_{j}=max(\delta_{A},\delta_{I})$$

   (24)

   Other possibilities (e.g. $\delta_{j}=\delta_{R}$) reduced $q$ too aggressively as $w$ increased, which resulted in very large errors.

In our case, the minimum value of $\delta$ occurs when $w=0$ which also coincides with $\Im(\delta)=0$ so the box size is easily set as equal to $L=r\sqrt{2\delta_{A}}$.

S91 describes a number of optimisations that can be used in the FGT with source-dependent scales:

1. 1.

   S91 introduces a single set of Taylor coefficients, positioned at the center of the $uv$ plane, to capture Gaussians with the largest scales. Visibilities larger than an upper scale size $\delta_{j}>b$ can be added to this single set, rather than the sets associated with the boxes. S91 recommends the upper scale size be at least 10% of the size of the evaluation region, in order to allow for practical precision in evaluating the large number of Taylor coefficients. In practice, most array configurations would have relatively few baselines with such a large $w$ that their convolution functions would span $>$ 10% of the $uv$ plane, so this optimisation has limited impact.

2. 2.

   S91 introduces a lower scale $a$ below which the Gaussian sum is directly evaluated. If there are only a few Gaussians with scales smaller than $a$, we can increase the box size, which reduces the memory and arithmetic capacity required, while incurring minimal penalty from direct evaluation. In practice, any interferometric observation is likely to be arranged such that the mean value of $w$ is approximately zero. Therefore, there are unlikely to be outlying, small Gaussians that would allow for increasing the box size. Once again, this optimisation has limited impact.

3. 3.

   S91 proposes doing direct evaluation for boxes containing only a few evaluation points. As we are evaluating on a regular grid, we have the same number of evaluation points inside each box, so applying this optimisation would degenerate into standard gridding.

4. 4.

   S91 proposes a method for avoiding needlessly accurate computation, i.e. for large scales (which are smoother than small ones) more boxes are updated, but smaller number of Taylor coefficients (i.e. the order $q$ described in Section 3.2) per box can be updated. This optimisation is useful and derives the most benefit for visibilities with large scales ($|w|$ large, see section 4.). We have incorporated this optimisation in Equation 19.

5. 5.

   S91 notes that the equation to calculate the order (Equation 19) overestimates the error by at sometimes two orders of magnitude. Memory bandwidth and operations can be saved by ‘cheating’ on the value of $q$ by decrementing $q$ by some value. This strategy is described in the section 3.4.1.

6. 6.

   The fact that each uv-cell contains a number of Taylor terms affords an additional axis available for parallelisation. The strategy is described in section 3.4.2.

Degridding is used as part of the major cycle of Cotton-schwab clean (Rau & Cornwell, 2011). Degridding involves taking the dot-product of the $uv$ plane with the convolution function shifted to the location of a visibility at an arbitrary location (i.e. not at the center of a $uv$ cell). It is the dual of gridding.

Degridding can be performed analogously to gridding using the fast Gauss transform with target-dependent scales (S91). All the same comments about optimisations for gridding (Section 3.4) apply as for degridding.

The key conceptual differences between classical gridding and FGT gridding are described in Table 1

In this section we derive the operations count, memory bandwidth and storage required to perform gridding and degridding with stored convolution funcitons. In order to put the two methods on equal footing in terms of imaging performance, we will use Gaussian anti-aliasing functions for the $w$-projection, rather than the commonly-used prolate spheroidal wavefunctions. This means that the imaging performance of the two approaches is essentially the same (except for the errors in truncating the Taylor series of the FGT) and the support size of the convolution functions is also easily computed.

We assume standard $w$ projection computes the convolution functions in advance and stores them in memory. We assume a convolution function of 1-D size for the $k$th $w$-plane $M_{\rm k}=k\Delta w$, where $k$ is an integer, and $\Delta w$ difference in size between different $w$ planes. Usually the cached convolution function is oversampled by a factor $\kappa$ of 4–8 in order to accurately grid visibilities whose coordinates are not exactly in the center of a $uv$ cell.

We assume that we truncate the convolution function when the real envelope reaches a value $\epsilon$, which is at a distance $t_{\epsilon}=\sqrt{-\delta_{R}\log\epsilon}$ from the centre of the convolution function. Therefore, the 1D support size in pixels, of the $k$th $w$-plane is (by substituting Equation 12)

A typical error pattern is shown in Figure 3. The absolute error in the gridded visibility contains two components. The circular component is due to the truncation to finite support, i.e. boxes outside the circle were not updated. The vertical and horizontal lines are due to truncation to a finite number of Taylor terms, with the largest error where the grid is evaluated near the boundary of adjacent boxes. The box structure of the gridding process can also be clearly seen in the FFT error, which also shows a box structure.

The peak error in the FFT of the visibility (i.e. the image plane), as a function of the $w$ for a range of algorithm parameterisations are shown in Figure 4. To begin with, we consider the case for $L~\sim 1$. Somewhat surprisingly, Equation 17, which describes the error in gridding a single visibility on $uv$ plane, also predicts the error the image plane for many parameterisations and values of $w$. For $p_{d}=0$, and small $w$, the errors are substantially smaller than the predicted $\epsilon$. As $w$ increases, $q^{\prime}$ is reduced, and the error generally matches the predicted $\epsilon$. For some parameterisations (e.g. $r=0.5$, $p=7$) there is a jump in error above predicted $\epsilon$ at a particular $w$.

For the $L~2$ case, for small $w$, the errors are at or below the theoretical limits in most cases, however in many cases the errors are considerably worse than for $L=1$. For example, for $w$ small in some cases (e.g. $L=2.5$, $p=18$, $p_{d}=0$) the error increases far beyond the theoretical limit. For large $w$, the errors are several orders of magnitude worse than predicted.

The measured processing time for the FGT vs. standard gridding is shown in Figure 5. Our implementation operates substantially slower than standard gridding for all parameterisations and ranges of $w$, and the operations rate is well matched by the our theoretical model for the operations count.

The results from testing our implementation are somewhat conflicting. On the plus side, Figure 4 clearly shows that, for small values of $w$, Equation 19 is overestimating the truncation error, and the value of $q^{\prime}$ being used is too large. This was the original motivation for introducing cheating, and it is clear that the error is still overestimated even for $p_{d}=3$ in some cases.

Unfortunately, for larger values of $w$, the fact that the error jumps above the predicted $\epsilon$, even with no cheating, implies that our empirical extension to the FGT to complex Gaussians, (in particular Equation 19, and Equation 24) is not correct. As a result of this incorrectness, our current cheating algorithm has limited value. Even with $p_{d}=1$, most parameterisations contain jumps in error of factors of few to 100 over the range of $w$. A different mapping between $w$ and $q^{\prime}$ is required to maintain a constant error across $w$.

Increasing the box size above $L>2.3$ appears to have catastrophic affects on the error for small $w$. We suspect this is due either to numerical instability in evaluating high-order polynomials, or more likely, the Hermite expansions simply fail to converge when the Gaussian width is less than some function of the box size.

Our implementation’s processing time is clearly compute-bound as our computing times match our theoretical model for the operations count almost exactly in most instances. The theoretical model suggests the operations count is dominated by the cost of computing the complex exponential, and we expect that is the case, although other overheads may also be responsible. Sadly, we were not able to approach the memory-bound performance that originally motivated this work.

In spite of the disappointing performance of our implementation, we hope that some of the ideas from this work could be extended in future. In particular:

• •

  The idea of having tuneable error is attractive. For arrays with high sidelobes in their synthesised beams (i.e. poor $uv$ coverage), large gridding errors can be made, as long as they remain below the clean threshold for the given major cycle.

• •

  The required memory bandwidth of our algorithm is only weakly dependant on $|w|$, so for arrays with very long baselines, our algorithm may be suitable.

• •

  The Hermite expansions were originally proposed by Greengard & Strain (1991) for real Gaussians. It may be that faster-converging expansions can be found for complex Gaussians - particularly in the case for large $w$ where the real envelope is smooth and the complex chirp contains many peaks.

• •

  The form of the convolution function is smooth at the centre and oscillates more rapidly as $t$ increases, with the fastest oscillations being damped at the edges by the envelope. Currently, our approach is to use the same order for all boxes, which leads to the largest errors being made away from the centre of the convolution function (Figure 3) where the oscillations are the fastest and largest. There may, therefore, be some value in varying the order $q^{\prime}$ as a function of the distance from the centre, thereby applying the majority of the computational effort where the errors are likely to be worst.

• •

  We implemented the polynomial evaluations in a naive manner. More sophisticated methods exist (Knuth, 1997, p. 485ff) that could reduce the operations count and improve numerical stability.

• •

  To calculate the complex exponential, we used the CEXP function from the standard C library, which is accurate but slow. Faster, approximate methods are available¹¹1http://gruntthepeon.free.fr/ssemath/ (e.g. (Schraudolph, 1998). For an implementation whose run time is dominated by CEXP such as ours, these methods could improve speed by a factor of a few, albeit at the cost of some accuracy.

• •

  Gaussian anti-aliasing functions are not well suited to imaging, as the alias rejection is not very good. The modern state of the art is to use prolate spheroidal wavefunctions (Schwab, 1980), that have better alias rejection. Unfortunately, the Hermite functions are not a suitable basis for Taylor expansion of the prolate spheroidal wavefunctions. If suitable analytic multipole expansions of the convolution of the complex chirp, and the prolate spheroidal wavefunctions can be found, then they can be applied with the fast multiple method (Greengard & Rokhlin, 1987; Greengard, 1988) to obtain a similar result as our FGT, but with better anti-aliasing properties.

• •

  The fast multipole method could also be used to compensate for the primary beam in the $uv$-plane (AW-projection, (Bhatnagar et al., 2008)). Once again, this would require analytic multipole expansions of the required convolution functions. While this may not be possible with an arbitrary primary beam, it may be more tractable if we assume Gaussian primary beams.

• •

  For problems where the compute time of $w$-projection is dominated by the time to compute convolution functions (by Fourier transforming a phase screen in the image plane), using standing gridding with Gaussian anti-aliasing functions could be preferable, as they can be directly and efficiently computed in the $uv$ domain.

• •

  Easily parallelised algorithms will clearly be important for large arrays such as SKA, where the computational work is likely to be spread over many thousands of nodes. Once the gridding problem has been distributed over the usual axes of frequency, pointing direction and polarisation, it is possible that problem may still be too large to be efficiently computed by a single node. For the gridding operation, the option of distributing and parallelising over Taylor coefficients adds an additional axis which an be exploited when distributing work among processors, over and above the traditional axes.