The EPIC Algorithm

From visibility gridding to electric-field gridding

In the optimal map-making, or software-holography, formalism (Morales & Matejek 2009; Morales 2011), each visibility is gridded onto the $uv$ plane with the holographic beam $B_{ab}=W_a{W}_b^{\ast }$. In the aperture plane this kernel is the cross-correlation ${\mathcal{W}}_a\star {\mathcal{W}}_b$ of the two aperture illuminations, and the image is

$$I^D(l,m)={\mathcal{F}}^{-1}[\sum _{a,b}({\mathcal{W}}_a\star {\mathcal{W}}_b)(\mathbf{u}-{\mathbf{u}}_{ab})V_{ab}].$$

Substituting $V_{ab}=⟨E_a{E}_b^{\ast }{⟩}_t$, the sum over antenna pairs factorizes into the correlation of two gridded electric-field distributions:

$$\sum _{a,b}({\mathcal{W}}_a\star {\mathcal{W}}_b)⟨E_a{E}_b^{\ast }⟩=⟨[\sum _a{\mathcal{W}}_a(\mathbf{u}-{\mathbf{u}}_a)E_a]\star [\sum _b{\mathcal{W}}_b(\mathbf{u}-{\mathbf{u}}_b)E_b]{⟩}_t.$$

By the convolution theorem, correlation in the aperture plane corresponds to pixel-wise multiplication in the image plane. This yields the EPIC imaging equation (Morales 2011; Thyagarajan et al. 2017, eq. 13; Krishnan et al. 2023):

$$\boxed{I(l,m)=⟨\left|{\mathcal{F}}^{-1}\right[{\sum }_a{\mathcal{W}}_a(\mathbf{u}-{\mathbf{u}}_a)E_a](l,m){|}^2{⟩}_t}$$

For each F-engine time sample and frequency channel, the algorithm therefore consists of four operations:

1. Gridding. The complex voltage $E_a$ of each antenna is placed on a regular aperture grid at position ${\mathbf{u}}_a={\mathbf{r}}_{a}f/c$ and convolved with the gridding kernel ${\mathcal{W}}_a$.
2. Spatial Fourier transform. A two-dimensional inverse FFT of the grid yields the instantaneous electric-field image $\tilde{E}(l,m)$.
3. Multiplication. Pixel-wise products of the electric-field images are formed.
4. Averaging. The products are accumulated over the integration time.

For two orthogonal polarizations, the pixel-wise (Hadamard) products of the two electric-field images form the coherency matrix of each pixel:

$$\mathbf{C}(l,m)=⟨\left(\begin{array}{cc}{\tilde{E}}_X{\tilde{E}}_X^{\ast }& {\tilde{E}}_X{\tilde{E}}_Y^{\ast }\\ {\tilde{E}}_Y{\tilde{E}}_X^{\ast }& {\tilde{E}}_Y{\tilde{E}}_Y^{\ast }\end{array}\right){⟩}_t.$$

Pseudo-Stokes parameters follow, e.g. $I=X{X}^{\ast }+Y{Y}^{\ast }$ and $V=2\mathrm{Im}(X{Y}^{\ast })$. The sign conventions depend on the orientation of the feeds.

Order of operations

In FX correlation the sequence is multiplication, averaging and then imaging. In EPIC it is imaging, multiplication and then averaging. The gridded electric field is a stochastic quantity with zero mean, so averaging cannot precede multiplication. Every channelized time sample must therefore be gridded and transformed: one transform every 40 µs per channel at LWA-SV. The cost of EPIC is consequently set by the Nyquist rate of the channelized data and is independent of the output cadence (Thyagarajan et al. 2017).

Autocorrelation removal

The squared sum $|{\sum }_a{\mathcal{W}}_a{E}_a{|}^2$ contains, besides the cross terms $a\ne b$, the self terms $a=b$. These are the antenna autocorrelations $|E_a{|}^2$, gridded with ${\mathcal{W}}_a\star {\mathcal{W}}_a$ at zero spacing. They contribute the bulk of the total power and add a smooth pedestal to the image.

They are removed by the following procedure:

1. The per-antenna autocorrelations, $|E_X{|}^2$, $|E_Y{|}^2$ and $E_X{E}_Y^{\ast }$, are accumulated alongside the images.
2. The accumulated autocorrelations are gridded at the grid center with the same kernels.
3. The resulting grid is Fourier transformed.
4. The transform is subtracted from the accumulated image.

The procedure costs one additional transform per integration. In epic++ (from version 0.8) it is performed within the imaging kernel.

Gridding kernels and grid correction

The gridding kernel ${\mathcal{W}}_a$ determines the primary-beam weighting applied to each antenna. With antenna-specific kernels EPIC correctly images telescopes composed of non-identical antennas. Thyagarajan et al. (2017) showed that assuming identical kernels for non-identical antennas biases flux estimates by more than $5\sigma$.

The kernel must also satisfy two practical requirements:

• It must suppress aliasing of emission from beyond the image boundary.
• It must remain compact, because the cost of gridding scales with the kernel area.

The kernel used in epic++ has the following properties:

• Kernel function. A zeroth-order prolate spheroidal wave function, the standard anti-aliasing kernel for convolutional gridding.
• Scaling. The kernel is scaled to the effective aperture of an LWA dipole, $A_{\mathrm{eff}}\approx 0.13{\lambda }^2$ (LWA Memo 65).
• Support. Typically $5\times 5$ cells.
• Normalization. The kernel is integrated over each grid cell and normalized to unit sum for each antenna, which conserves flux.

The image-plane taper imposed by the kernel is removed by multiplying the image by $1/|\mathcal{F}[\bar{\mathcal{W}}]{|}^2$, where $\bar{\mathcal{W}}$ is the antenna-averaged kernel.

Grid geometry

For an $N\times N$ transform, the image-plane pixel scale and the aperture-plane cell size $\delta u$ (in wavelengths) are related by $\Delta l=1/(N\delta u)$. In epic++ the cell size is derived from the requested angular resolution at zenith, ${\theta }_{\mathrm{res}}$:

$$\delta u=\frac{1}{2N\sin ({\theta }_{\mathrm{res}}/2)},\delta u\approx 0.42\lambda \text{for }N=128,{\theta }_{\mathrm{res}}=1.056^{\circ }.$$

A cell size below $\lambda /2$ ensures that the image spans $|l|\le 1$, so the full hemisphere is imaged without aliasing of the horizon. The conventional LWA all-sky pixel scale is $130^{\circ }/N$ (Memo 012).

Antenna positions are expressed in wavelengths, ${\mathbf{u}}_a={\mathbf{r}}_{a}f/c$. Each channel therefore has its own set of grid positions and kernel weights. The pixel scale in $(l,m)$ is fixed by $\delta u$ and is common to all channels.

The coplanar approximation

The derivation assumes that all antennas lie in a common plane and neglects the $w$ term, $e^{-2\pi iw(\sqrt{1-l^2-m^2}-1)}$. LWA-SV is nearly coplanar, and epic++ applies a per-antenna phase corresponding to the height $z_a$ of each stand. Over the full hemisphere, residual wide-field errors nevertheless reach approximately 10% within the beam (Kent et al. 2020).

Two corrections have been investigated:

• w-projection. The $w$ correction is incorporated into the per-antenna gridding kernels.
• Direct Fourier transform. The FFT is replaced by a direct Fourier transform, implemented as a dense matrix multiplication (Kent et al. 2020). This is exact, at a cost of $\mathcal{O}(N_a{N}_{\mathrm{pix}})$.

Computational cost

Per channel, per polarization product and per F-engine time sample, the costs of EPIC and FX correlation are approximately (Thyagarajan et al. 2017)

$$C_{\mathrm{EPIC}}\approx \underset{\text{gridding}}{\underbrace{N_a{N}_k}}+\underset{\text{FFT}}{\underbrace{\beta (4N_g){\log }_2(4N_g)}}+\underset{\text{multiplication}}{\underbrace{4N_g}},C_{\mathrm{FX}}\approx \frac{N_a(N_a-1)}{2}.$$

The symbols are:

• $N_k$, the number of cells in the gridding kernel;
• $N_g$, the number of grid cells (the factor of four accounts for zero-padding);
• $\beta \approx 5$, a constant characterizing the FFT.

Imaging the full antenna field of view requires a cell size comparable to the antenna size, so $N_g\approx b_{\max }^2/A_{\mathrm{ant}}$. The cost of EPIC is therefore governed by the area of the aperture in units of the antenna area, not by the number of antennas. For a densely filled aperture, $N_g\sim N_a$, and EPIC scales as $\mathcal{O}(N_a{\log }_2{N}_a)$.

Thyagarajan et al. (2017) applied this model to a range of radio telescopes:

Favored architecture	Telescopes
EPIC	HERA-331, HERA-6769, CHIME, HIRAX, SKA1-Low stations, expanded LWA configurations
FX	MWA, LOFAR, small HERA configurations
Near break-even	LWA1, LWA-SV

The output data rates differ more strongly than the arithmetic:

$$r_{\mathrm{EPIC}}\sim \frac{4N_g}{\Delta t}\frac{\Delta B}{\Delta f},r_{\mathrm{X}}\sim \frac{N_a(N_a-1)}{\Delta t}\frac{\Delta B}{\Delta f}.$$

At LWA-SV the arithmetic costs of the two approaches are comparable. Short integration times favor EPIC, because the volume of the image output is independent of $N_a$, and images at every time sample are intrinsic to the method.

On GPUs, an implementation of EPIC as separate gridding, transform and multiplication kernels is limited by memory bandwidth (Kent et al. 2019; Memo 009). FX correlation, by contrast, is limited by arithmetic throughput. The fused implementation described in GPU Kernel addresses this limitation.

Comparison of architectures

	FX correlator	EPIC
Correlation	All antenna pairs	Pixel-wise multiplication after FFT
Cost per sample	$\mathcal{O}(N_a^2)$	$\mathcal{O}(N_g{\log }_2{N}_g)$
Output per integration	$N_a^2$ visibilities	$N_{\mathrm{pix}}$ pixels
Time averaging	After multiplication	After multiplication; every sample is imaged
Calibration	Applied to visibilities, possibly offline	Applied to voltages before gridding
Antenna layout	Arbitrary	Arbitrary; most efficient for filled apertures
Optimal regime	Sparse antenna layouts, long integrations	Densely packed large-$N$ apertures, short integrations