# matvis: a matrix-based visibility simulator for fast forward modelling of many-element 21 cm arrays

Piyanat Kittiwisit, Steven Murray, Hugh Garsden, Philip Bull, Michael

Wilensky, Christopher Cain, Aaron Parsons, Jackson Sipple, Tyrone Adams,
James Aguirre, et al.

► To cite this version:

Piyanat Kittiwisit, Steven Murray, Hugh Garsden, Philip Bull, Michael Wilensky, et al. matvis: a matrix-based visibility simulator for fast forward modelling of many-element 21 cm arrays. RAS Techniques and Instruments, 2025, 4, 10.1093/rasti/rzaf001. hal-05392536

HAL Id: hal-05392536

https://hal.science/hal-05392536v1

Submitted on 5 Dec 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

matvis: a matrix-based visibility simulator for fast forward modelling of many-element 21 cm arrays

Piyanat Kittiwisit,¹,²* Steven G. Murray,³,⁴ Hugh Garsden,⁵,⁶ Philip Bull,¹,⁶ Michael J. Wilensky,⁶ Christopher Cain,⁴ Aaron R. Parsons,⁷ Jackson Sipple,⁷ Tyrone Adams,² James E. Aguirre,⁵,⁸ Rushelle Baartman,² Adam P. Beardsley,⁴,⁹ Lindsay M. Berkhout,⁴,¹⁰ Gianni Bernardi,²,¹¹,¹² Tashalee S. Billings,⁸ Judd D. Bowman,⁴ Richard F. Bradley,¹³ Jacob Burba,⁶ Steven Carey,¹⁴ Chris L. Carilli,¹⁵ Kai-Feng Chen,¹⁶,¹⁷ Samir Choudhuri,⁵,¹⁸ Tyler Cox,⁷ David R. DeBoer,¹⁹ Eloy de Lera Acedo,¹⁴ Matt Dexter,¹⁹ Joshua S. Dillon,⁷ Nico Eksteen,² John Ely,¹⁴ Aaron Ewall-Wice,⁷,²⁰ Nicolas Fagnoni,¹⁴ Steven R. Furlanetto,²¹ Kingsley Gale-Sides,¹⁴ Bharat Kumar Gehlot,⁴ Brian Glendenning,¹⁵ Adelie Gorce,²² Deepthi Gorthi,⁷ Bradley Greig,²³ Jasper Grobbelaar,² Ziyaad Halday,² Bryna J. Hazelton,²,²⁵ Jacqueline N. Hewitt,¹⁶,¹⁷ Jack Hickish,¹⁹ Daniel C. Jacobs,⁴ Alec Josaitis,¹⁴ Nicholas S. Kern,¹⁶,¹⁷† Joshua Kerrigan,²⁶ Honggeun Kim,¹⁶,¹⁷ Matthew Kolopanis,⁴ Adam Lanman,²⁶ Paul La Plante,⁷,⁸,²⁷ Adrian Liu,⁷,²⁸ Yin-Zhe Ma,²⁹,³⁰ David H. E. MacMahon,¹⁹ Lourence Malan,² Cresshim Malgas,² Keith Malgas,² Bradley Marero,² Zachary E. Martinot,⁸ Lisa McBride,²²,²⁸ Andrei Mesinger,³ Mathakane Molewa,² Miguel F. Morales,²⁴ Tshegofalang Mosiane,² Chuneeta Devi Nunhokee,³¹,³² Hans Nuwegeld,² Robert Pascua,⁷,²⁸ Yuxiang Qin,³³ Eleanor Rath,¹⁶,¹⁷ Nima Razavi-Ghods,¹⁴ James Robnett,¹⁵ Mario G. Santos,¹,² Peter Sims,⁴ Saurabh Singh,²⁸,³⁴ Dara Storer,²⁴ Hilton Swarts,² Jianrong Tan,⁸ Nithyanandan Thyagarajan,¹⁵,³⁵ Pieter van Wyngaarden,² Zhilei Xu¹⁶ and Haoxuan Zheng¹⁷

*Affiliations are listed at the end of the paper*

Accepted 2025 January 7. Received 2024 December 25; in original form 2023 December 28

ABSTRACT

Detection of the faint 21 cm line emission from the Cosmic Dawn and Epoch of Reionization will require not only exquisite control over instrumental calibration and systematics to achieve the necessary dynamic range of observations but also validation of analysis techniques to demonstrate their statistical properties and signal loss characteristics. A key ingredient in achieving this is the ability to perform high-fidelity simulations of the kinds of data that are produced by the large, many-element, radio interferometric arrays that have been purpose-built for these studies. The large scale of these arrays presents a computational challenge, as one must simulate a detailed sky and instrumental model across many hundreds of frequency channels, thousands of time samples, and tens of thousands of baselines for arrays with hundreds of antennas. In this paper, we present a fast matrix-based method for simulating radio interferometric measurements (visibilities) at the necessary scale. We achieve this through judicious use of primary beam interpolation, fast approximations for coordinate transforms, and a vectorized outer product to expand per-antenna quantities to per-baseline visibilities, coupled with standard parallelization techniques. We validate the results of this method, implemented in the publicly available *matvis* code, against a high-precision reference simulator, and explore its computational scaling on a variety of problems.

**Key words:** software: development, public release, simulations – algorithm – techniques: interferometric – cosmology: diffuse radiation, dark ages, reionization, first stars.

1 INTRODUCTION

The 21 cm line emission from neutral hydrogen during the Cosmic Dawn and Epoch of Reionization (EoR) is currently one of the most challenging observational targets for radio astronomers. According

<footnote>* E-mail: piyanat.kittiwisit@gmail.com</footnote>

<footnote>† NASA Hubble Fellow</footnote>
to theoretical models (e.g. Cohen, Fialkov & Barkana 2018), the total line emission is expected to be quite weak, on the order of tens to hundreds of mK, with fluctuations of $\Delta_{21}^2 \lesssim 10^3$ mK$^2$ (where $\Delta_{21}^2$ is the 'dimensionless' power spectrum at a given scale and redshift). In comparison, the Galactic synchrotron foreground emission is many orders of magnitude brighter, with a much larger variance (Santos, Cooray & Knox 2005). Contributions from extragalactic sources are also considerably brighter than the target 21 cm signal. As such, 21 cm surveys must achieve both high sensitivity, to reach an adequate signal-to-noise ratio to detect the 21 cm fluctuations, and excellent instrumental stability and calibration precision, in order to prevent spurious leakage of some fraction of the bright foregrounds into the 21 cm signal modes.

Even purpose-built 21 cm radio interferometer arrays have complex instrumental responses and systematics that interact with the sky and contribute to the difficulty of characterizing the instrument sufficiently well to mitigate the foreground leakage problem mentioned above. These include the frequency- and direction-dependent polarized primary beam pattern of the antenna, and the time- and frequency-dependent receiver bandpass required to calibrate the data (Barry et al. 2016; Barry & Chokshi 2022), the idiosyncratic properties of the receiver chain in the presence of radio frequency interference (RFI; Wilensky et al. 2023), environmental thermal variations (Gehlot et al. 2021), reflection and retransmission of signals over the air and through cables (Kern et al. 2019, 2020; Josaitis et al. 2022), among others. In addition, the data analysis tools used to mitigate imperfections in the data (e.g. missing channels due to RFI contamination, filtering of cable reflections etc.) often have their own non-ideal behaviours, including the possibility of suppressing or filtering some of the 21 cm signal itself (Cheng et al. 2018; Kolopanis et al. 2019; Pagano et al. 2023).

Physical (electronic and mechanical) approaches to instrumental characterization may be useful, but also have their limitations in terms of understanding the ultimate impact of these systematics on the processed data. For example, primary beams can be characterized through a variety of physical methods, including photogrammetry, overflights by drones carrying standard sources and so on, but may still fail to capture the precise shapes of far sidelobes with sufficient accuracy (Jacobs et al. 2017). Experimental approaches that isolate a single component of the system may also be unable to properly constrain the couplings and interactions between different components that give rise to some effects.

As a result, simulation-based, 'forward-modelling,' approaches, for which realistic mock data are produced to study the complex process of observation and data analysis, are now a vital part of all 21 cm experiments. An important advantage of simulations is that they can be used to develop increasingly complex 'what if' scenarios for the systematic effects that are identified, particularly for effects that would be close to (or even below) the noise level of the available data. Simulations can also be used as part of holistic efforts to 'validate' the data analysis pipeline, allowing ones to explicitly determine levels of signal loss or other spurious features in the processing of the data as each step of the pipeline proceeds (Aguirre et al. 2022). A disadvantage of simulation-based approaches is that they require models to be available for each system component, each systematic effect, and indeed the sources of emission on the sky, which can lack realism. They can also be computationally expensive, particularly when simulations for the whole array are required – not only must we run the data analysis pipeline on a large volume of simulated data but we must also generate and store the simulated data itself.

The first step in the production of the mock data for the simulation-based study is the simulation of the complex visibility measured by

radio interferometer arrays, and the basis of all of visibility simulators is a suitable simplification or special case of the radio interferometer measurement equation (RIME), which, as its name suggests, is the mathematical description of what an interferometer measures. Many visibility simulators have been developed, although only a small number of them are well documented. In general, they can be grouped into four categories:

(i) Simulators that directly evaluate the original RIME formalism from Hamaker (2000, 2006) and Hamaker, Bregman & Sault (1996), which include simulators such as OSKAR (Mort et al. 2010), WODEN (Sullivan et al. 2012), healvis (Lanman & Kern 2019), and pyuvsim (Lanman et al. 2019).

(ii) Simulators that evaluate the updated, generalized, RIME from Smirnov (2011a, b, c, d), which can be used to describe all effects along the propagation path of the sky signals to the telescope, making the simulation results more suitable for testing direction-dependent calibration, for example. Such simulators include MeqTrees (Ordam & Smirnov 2010) and Montblanc (Perkins et al. 2015).

(iii) Simulators that utilize the 'm-mode' spherical harmonic formalism in Shaw et al. (2014, 2015), which is more suitable for wide-field, transit, radio telescopes that use stationary antenna elements to capture the horizon-to-horizon sky signal as it transits across the arrays. Examples include driftscan (Shaw et al. 2023) and RIMEz (Martinot 2022).

(iv) General-purpose radio interferometric data analysis software with visibility simulation capability, including FHD (Sullivan et al. 2012; Barry et al. 2019), CASA (The CASA Team 2022), and SAGECal (Spreeuw et al. 2020a, b, c). These tools usually simulate the visibility as part of the calibration model building step by evaluating either the original or generalized RIME in the Fourier transform domain.

Apart from the flavours of the RIME, all visibility simulators can vary significantly in realism, performance, and importantly accuracy. The realism of the simulation is largely governed by the realism of the models used for each component of the RIME, such as the primary beam response, array layout and configuration, and sky intensity distribution. Inclusion and implementation of specific models of these components by a simulator are largely dictated by its targeted applications and the experiments for which the simulators were originally developed. For example, OSKAR supports the use and customization of beamformed stations, as it was developed as part of the Square Kilometre Array (SKA) study, while WODEN supports shapelet models of bright extended sources that are crucial for calibrating the Murchison Widefield Array (MWA) 21 cm observations (Line et al. 2020). Further instrumental and observational effects are usually added on top of the RIME as post-processing steps; there is no need to simulate RFI pulses through the visibility simulator machinery in most cases, for example.

The performance of the simulators is governed by a few factors: the scale of the simulation, the particular choices of the algorithms and technology used to implement the RIME, and the approximations that are used in evaluating the RIME and its constituent component models. For instance, simulating 'realistic' data volumes of a very large array, consisting of many antennas with wideband receivers, observing for a potentially very long period of time, can soon become prohibitively expensive, especially if multiple runs (e.g. to model different scenarios) are required. The performance of such a simulation can often be improved by employing parallelization (e.g. pyuvsim), vectorization, or GPU acceleration (e.g. OSKAR, WODEN).
One may also utilize approximations, such as analytical forms of the sky intensity and primary beam response that can be evaluated quickly, to improve performance at the cost of realism, but these approximations can also be important sources of inaccuracy. No simulator will be perfectly accurate, given limits to numerical precision and the uncertainty in the models of different components, but certain choices of algorithm and functional basis can give rise to inaccuracies that can potentially exceed the size of the effects that one is trying to model. Balancing realism, performance and accuracy can be complicated and depend on specific use cases as we have touched on. Regardless, it is important that the impact of the approximations and the overall accuracy of the simulator is well understood, validated, and documented for both users and other developers in the astronomical community, particularly when the simulation will be used to validate highly sensitive analyses of real data.

In this paper, we introduce the *matvis* visibility simulator and conduct rudimentary validation of its outputs. *matvis* has been designed to perform fast simulation of the visibility response of many-element, drift-scanning, 21 cm arrays at low frequency, in particular the Hydrogen Epoch of Reionization Array (DeBoer et al. 2017; HERA Berkhout et al. 2024). The primary use case of *matvis* outputs has been the validation of the HERA analysis pipeline via simulated mock data (Aguirre et al. 2022). However, the underlying architecture of the simulator is more general and applicable to any radio interferometric array at any frequency. Similarly, its outputs can also be used in other applications with some limitations as we will discuss. We elaborate on the design principles and use cases of *matvis* in Section 2. In Section 3, we provide an overview of the basic calculation of a visibility simulator, as well as the specific algorithm and the approximations and assumptions used by *matvis*. Section 4 describes the software aspects of the publicly available *matvis* package hosted on GitHub. In Section 5, we perform a rudimentary ‘validation’ of *matvis* against *pyuvsim*, an alternative simulator that has been designed to carefully and correctly evaluate each step of the RIME, at the cost of decreased computational performance. We explore effects of *matvis* coordinate transformation scheme and show that *matvis* can attain accuracy comparable to *pyuvsim*. We also demonstrate the limitation of the point-like sky assumption in this section. Lastly, we study the computational scaling of *matvis* to understand its performance in simulating ever-growing data sets on large arrays in Section 6 and conclude in Section 7.

## 2 DESIGN PRINCIPLES AND USE CASES

The original motivation for *matvis* was the need for a well-documented, well-tested, and fast visibility simulator to produce mock data for validating the analysis pipeline of the Hydrogen Epoch of Reionization Array (HERA; DeBoer et al. 2017; Berkhout et al. 2024). This validation process has been deliberately described in Aguirre et al. (2022), for which its main goals of this process is to quantify potential signal loss from analysis steps through series of tests with mock data. In particular, pipeline validation requires production of mock data with ‘realistic physical properties’ at the scale of the true data. What constitutes these properties is an experiment-specific question, but generally includes reasonable scale-dependent magnitudes for different sky components (diffuse foregrounds, extragalactic point sources, and the 21 cm signal), with special attention to reproducing the spectral and temporal characteristics of the instrument.

This use case leads to the design principles for *matvis* that naturally priorities the need for computational efficiency over detailed

realism. In contrast, a use case of calibrating data would require exquisite realism and accuracy (both in the simulator and the input component models). To meet these requirements, *matvis* introduces a matrix-based approach to evaluate the RIME, which is able to naturally take advantage of GPU acceleration. The primary novelty of *matvis* is that it treats the RIME as an outer product over antennas, enabling the RIME for all baselines and polarizations at a single time and frequency to be computed as a single matrix product. The discretization of the sky required to enable this matrix-based approach is predicated on being able to treat small, local equal-solid-angle, areas as delta-function point sources, essentially approximating the RIME angular integral as a Riemann sum, i.e. the ‘point-source approximation’. This treatment of the sky components has been elaborately studied and validated against analytical calculations in the *pyuvsim* simulator (Lanman et al. 2019; Lanman, Murray & Jacobs 2022), for which we compare and discuss in the paper. A further novelty (and approximation) of *matvis* is that it utilizes a coordinate transformation scheme that combines accurate-but-slow *astropy* (Astropy Collaboration 2013, 2018, 2022) algorithms with simple-and-fast trigonometric-based calculations to significantly improve computational performance for simulations requiring multiple time-snapshots although this impact accuracy of the simulation.

These computational tricks employed by *matvis* does not mean that it completely sacrifices accuracy. In fact, the accuracy of *matvis* is equivalent to that of *pyuvsim* when simulating a single snapshot observation. Both simulators belongs to the first category of visibility simulators that we have described – i.e. those that directly evaluate the RIME. Both assume the point-source approximation for sky models and use *E*-field primary beams (interpolated to given source positions and frequencies) in the calculations. The main differences are that *pyuvsim* (i) directly calculates a full Jones-matrix contribution (Hamaker 2006) from each point source to each baseline, and (ii) uses *astropy* for coordinate transformation, providing the best possible source position accuracy. In contrast, *matvis* uses (i) a matrix-based optimization to improve speed of the calculations, and (ii) a geometric rotation scheme in place of *astropy* for coordinate transformation to trade some accuracy for more speed. The use of geometric rotations in *matvis* does contribute to a lower accuracy than *pyuvsim* when simulating more than one time-step at a time. We will demonstrate and discuss in Section 5. *matvis* is also similar to *pyuvsim* in that it is more general in its applicability and this is reflected in the flexibility with which input components can be specified in the simulator. In terms of performance, however, we show in Section 6 that *matvis* far outperforms *pyuvsim*.

In context with other simulators that we have listed above, *matvis* is likely in a similar performance regime as the GPU-accelerated WODEN and OSKAR. However, both are highly specific. For example, OSKAR is built with simulating the SKA in mind, while WODEN has support only for built-in primary beams, which are limited to some analytic models and the MWA beams. Instead, *matvis*, like *pyuvsim*, accept any external primary beam, and any sky model that can be represented as a list of spatially compact components.

Despite, there are a few limitation to the current version of *matvis* (1.2.1) presented in this paper. First, it is currently not possible to simulate from fully polarized sky models although polarized beam is supported. This is due to the matrix-based algorithm that we will discuss shortly. Secondly, it currently only support simulating drift-scan observations. This limitation is the design choice specific to HERA and most of other 21 cm arrays, which are drift-scan instruments. Modification of our current algorithm
to include supports for polarized sky and tracking observations is possible with some caveates. We will discuss them in Section 3.

## 3 COMPUTATIONAL SCHEME

In this section, we introduce the algorithm used by matvis to simulate visibilities. In common with all visibility simulators, a variety of choices need to be made in terms of how to model the sky (particularly whether and how to decompose it into a set of discrete elements); how to perform the integral over the whole sky required by the RIME; how to model the primary beams of the antenna elements; what coordinate systems to work in; and how to approximate or optimize computationally demanding steps of the process.

### 3.1 Basic framework

The purpose of the visibility simulation is to evaluate some form of the RIME (e.g. Hamaker et al. 1996; Smirnov 2011a) at the desired frequencies and observing times for a given sky model and instrumental model. Let $\Omega$ be celestial coordinates for a point on the sky (e.g. equatorial coordinates in RA and DEC); then the RIME for a snapshot visibility (at a particular time $t$, and frequency $\nu$) can be written in a reasonably general form as (Kohn et al. 2019)

$$V_{ij}(v, t) = \int_{4\pi} \mathbf{A}_i(v, t) \mathbf{C}(v) \mathbf{A}_j^{\dagger}(v, t) \times \exp(-2\pi v i \mathbf{b}_{ij} \cdot \mathbf{n}/c) \mathrm{d}^2\Omega. \quad (1)$$

Here, $V_{ij}$ is the complex visibility for a baseline formed from antennas $i$ and $j$, and we have suppressed dependence on sky coordinate inside the $\mathbf{A}$ and $\mathbf{C}$ terms in the integral for notational clarity. In this expression, $\mathbf{A}_i$ denotes the polarized primary (E-field) beam for antenna $i$, which is an $N_{\text{feed}} \times N_{\text{ax}} \equiv 2 \times 2$ direction-dependent Jones matrix, whose components are the antenna patterns for each feed $p$, and two unit vectors needed to model the response to an incident E-field from spherical polar direction $k \in (\theta, \phi)$. Throughout, we will assume the case in which $p \in \{X, Y\}$, i.e. in which there are two perpendicular feeds, though matvis trivially deals with the case in which there is a single feed. Note that $\mathbf{A}$ is a special case of the full Jones matrix that determines the linear transformation of the electric field that propagates from the sky to the antenna. In the nomenclature of Smirnov (2011a), $\mathbf{A}$ represents the full direction-dependent part of the Jones matrix (sans the phase difference term, which is captured explicitly by the exponential fringe factor in equation 1). In this paper, we consider $\mathbf{A}$ to be fully specified by the primary beam of the instrument multiplied by a horizon mask. However, in principle it could include direction-dependent effects such as ionospheric refraction and Faraday rotation. While such effects are in principle able to be included in simulations with matvis, to do so would require manually altering the input primary beam model, and is not considered a natural use case for matvis at this time. We will thus proceed to refer to $\mathbf{A}$ as the ‘primary beam’ within this paper.

Since we have chosen to represent the RIME as an integral over celestial coordinates, $\mathbf{A}$ is represented as evolving over time, as the earth rotates. $\mathbf{C}$ denotes the coherence matrix, which models the polarized flux distribution on the sky, which we assume is static in celestial coordinates.<fnref n="1" /> In terms of the Stokes parameters, $\mathbf{C}$ is written

as (Smirnov 2011a)

$$\mathbf{C} = \begin{pmatrix} \langle E_{\theta} E_{\theta}^{*} \rangle \langle E_{\theta} E_{\phi}^{*} \rangle \\ \langle E_{\theta}^{*} E_{\phi} \rangle \langle E_{\phi} E_{\phi}^{*} \rangle \end{pmatrix} = \begin{pmatrix} I + Q & U + iV \\ U - iV & I - Q \end{pmatrix}. \quad (2)$$

Finally, the exponential term, which shapes the fringe pattern of the radio interferometer, includes a scalar product between the baseline vector $\mathbf{b}_{ij}$ connecting the antennas, and a unit vector $\mathbf{n}$ towards a particular point on the sky. The integral over solid angle is performed over $4\pi$ steradians, so that in principle the whole sky contributes to each complex visibility. In practice, only emission above the horizon at any particular time actually contributes, and this ‘horizon cutoff’ can be bundled into the primary beam.<fnref n="2" />

From here, all visibility simulators must make at least one approximation: we must make a choice about which discrete basis the sky model should be represented in, so that the integral in equation (1) can be discretized. This is due to the fact that we do not have analytic forms of the sky or primary beam model as a function of sky coordinate. A very natural basis is the ‘point-source approximation’, i.e. a choice of a set of points on the sky to which we assign the total flux density from their surrounding local region, and this is the representation that matvis adopts. Note that while this form of discretization is natural, it is by no means the only one. Discretizations in other kinds of bases exist, for example spherical harmonics (Shaw et al. 2014, 2015; Martinot 2022).

Furthermore, matvis makes a slightly more stringent point-source approximation for each sky pixel; it explicitly assumes that the full flux density from a given region defined as a ‘pixel’ really does come from an infinitesimal point source located within the region, at which the primary beam and fringe term are evaluated. This point-source representation is in fact perfect if the sky purely consists of unresolved point-sources. However, if the sky contains extended sources, this approximation will inherently miss fluctuations on scales smaller than the sampling points (whether these fluctuations are intrinsic to the sky, or they come from the beam or fringe). In principle, this limitation can be overcome by following the Nyquist sampling principle, i.e providing input point sources dense enough on the sky that the smallest resolvable scales of the array are sampled (Lanman et al. 2022). We investigate the ramifications of this approximation on the accuracy of matvis in detail in Section 5.3.

Summarily, matvis uses a ‘point source’ representation, treating the sky as a set of point sources with arbitrary frequency spectra binned into discrete frequency channels. It does not enforce any particular ‘pixelization’, other than the condition that the sum of all pixels’ flux densities equals the total flux density of the sky. A common example of such a pixelization for diffuse emission is HEALPIX<fnref n="3" /> (Górski et al. 2005), but we emphasize that matvis does not require a HEALPIX representation, or even an equal-area pixel decomposition.<fnref n="4" />

<footnote n="1">¹Of course, we could have represented the RIME in topocentric coordinates, in which case $\mathbf{A}$ would have been static, and $\mathbf{C}$ would evolve, but this representation is conceptually simpler for our algorithm presented later.</footnote>

<footnote n="2">²The approach of bundling horizon effects into the beam is actually useful in terms of implementation because in practice the horizon may not always be perfectly flat (Bassett et al. 2021), and therefore should be left flexible for a user to modify.</footnote>

<footnote n="3">³http://healpix.sf.net/</footnote>

<footnote n="4">⁴Note that several well-used simulators adopt this same essential discretization and associated approximation, including pyuvsim (Lanman et al. 2019), CASA (The CASA Team 2022) miriad (Sault, Teuben & Wright 1995). The latter two, which are more focused on narrow-field-of-view instruments and imaging, also adopt the specific pixelization of a regular $(l, m)$ grid, and express their intensities in the commonly used units of Jy pixel$^{-1}$, which is only possible when all pixels have the same solid angle.</footnote>
From here on, we will refer to each component of the sky as a ‘source’ or ‘point source’ with the understanding that it may either be a truly physical point source, or an approximation of a region of solid angle as such. With this explicit assumption of a point-like sky model in place, we can write the RIME, at a particular frequency, as

$$V_{ij}^{pq}(v_a, t) = \sum_{nkk'} A_i^{pk}(v_a, \theta_{hrz}^{(n)}(t)) \cdot A_j^{qk'*}(v_a, \theta_{hrz}^{(n)}(t)) \\ \times \mathbf{C}_{kk'}^{(n)}(v_a, \theta_{eq}^{(n)}) \exp(-2\pi i v_a \tau_{ij}^{(n)}). \quad (3)$$

This gives the complex visibility $V_{ij}^{pq}$ for antenna pair i, j and polarization pq. The sum is over all point sources in the sky model, and pairs of E-field axes. $A_i^{pk}$ is the complex primary beam response for antenna i on feed p to the E-field direction k. The primary beam is evaluated at $\theta_{hrz}^{(n)}(t)$, the apparent angular position of each point source n in a horizontal (topocentric) coordinate system at time t. The phase factor is a function of a baseline- and source position-dependent delay,

$$\tau_{ij}^{(n)} = \mathbf{X}_{hrz}^{(n)}(t) \cdot \mathbf{b}_{ij}/c, \quad (4)$$

where $\mathbf{X}_{hrz}^{(n)}(t)$ is the Cartesian unit vector pointing towards the source at an observed time t in the same horizontal coordinate system as the baseline vector $\mathbf{b}_{ij}$ between the centres of antennas i and j, and c is the speed of light. We have chosen to write the source position unit vectors in the coordinate frames and representations that are more practical to evaluate for the different terms. We also note that the only approximation made in arriving at equation (3) has been the ‘point source approximation’ that discretizes the spatial integral.

## 3.2 Matrix-based computation of the RIME

To obtain the visibilities, equation (3) must be calculated for all discrete frequencies, and at all requested discrete times, for all antenna pairs i, j. This becomes computationally prohibitive for large arrays. For reference, the completed HERA array will have 350 antennas (up to ≈ 61 000 baselines), 1535 frequency channels, and so observing with a ≈ 10 s cadence over, say, 8 h per night, we can see that equation (3) must be evaluated around ~ 2.7 × 10¹¹ times to simulate even a single night of observing in full. Multiplying this by the number of individual sources in the sky model (e.g. the ≈ 300 000 point sources in the GLEAM catalogue), this grows to ~ 8 × 10¹⁷. Clearly, the computational problem is large in scale.

A common approach to calculating equation (3) is to perform a summation over frequencies, times and baselines, and for each, computing the primary beam product $A_i A_j^*$ and the fringe term $\exp(-2\pi i v_a \tau_{ij})$ for each source, before finally summing over the sources (this approach is taken, for example, by pyuvsim). This has two non-optimal ramifications for computational performance: first, naively this involves redundant interpolations of each antenna beam, since each antenna enters ≈ $N_{ants}$ times throughout the loop over baselines. This can be avoided, even within this common baseline-based approach, by first computing all antenna-beam interpolations, and then drawing from these pre-computed values within the baseline loop. The second inefficiency in this approach is that the final sum over sources is a relatively small computational task (essentially a vector-vector inner product). It is well-known that modern computational algorithms for matrix-matrix multiplication are far more efficient than the equivalent loop over vector-vector inner products: for a square $N \times N$ matrix M, the computation of $P = MM$ via looping over all $N^2$ elements of P and computing $P_{ij} = \sum_k M_{ik} M_{kj}$ scales as $N^3$, whereas modern matrix-matrix

multiplication algorithms can scale as well as ~ $N^{2.3}$ (Le Gall & Urrutia 2018; Vassilevska Williams et al. 2023).

The fundamental approach of the matvis algorithm is to express equation (3) as a tensor product (and, in some special cases, a matrix product) of tensors that are sized proportional to $N_{ants}$ instead of $N_{bls}$. The tensor product ‘expands’ the calculated quantities into a full set of cross-correlations for all baselines. This significantly reduces the number of repeated calculations, especially for a large array, and makes good use of the vectorization capabilities of modern linear algebra frameworks and compute architectures (including GPUs). The downside is that all antenna-antenna correlations must be performed to make use of this approach, so it is not possible to be selective about which baselines are formed. For problems where the number of baselines to be simulated is close to the maximum for the array, $N_{bls} \approx N_{ants}(N_{ants} - 1)/2$, this approach will be highly efficient.

There are two main tricks required to implement this approach. The first trick is to re-write the exponential fringe term as a product of antenna-based factors. Letting $\mathbf{y}_n = 2\pi v_a \mathbf{X}_{hrz}^{(n)}(t)/c$ for compactness, we write

$$e^{-i\mathbf{y}_n(t)\mathbf{b}_{ij}} = e^{-i\mathbf{y}_n(t)\mathbf{x}_i} e^{+i\mathbf{y}_n(t)\mathbf{x}_j} \quad (5)$$

$$\equiv \mathbf{F}_{in}(v_a, t) \mathbf{F}_{jn}^*(v_a, t). \quad (6)$$

The second trick is to note that since the $2 \times 2$ matrix $\mathbf{C}_n$ for any particular source n is Hermitian, it can be written

$$\mathbf{C}_n = \mathbf{M}_n \mathbf{M}_n^\dagger. \quad (7)$$

We give the form of M in Appendix B. With this in place, we can re-write equation (3) as the tensor product

$$V_{ij}^{pq}(v_a, t) = \sum_{k,k',n} \mathbf{A}_{ipkn} \mathbf{F}_{in} \mathbf{C}_{kk'n} \mathbf{F}_{jn}^* \mathbf{A}_{jqk'n}, \quad (8)$$

$$= \sum_{k,k',k'n} \mathbf{A}_{ipkn} \mathbf{F}_{in} \mathbf{M}_{kk'n} \mathbf{M}_{k'k'n}^* \mathbf{F}_{jn}^* \mathbf{A}_{jqk'n}, \quad (9)$$

where each of the factors in the sum is understood to be dependent on the frequency and time. Here i and j index antennas, p and q index the (up to) two feeds on each antenna, k and k' index the E-field component, and n indexes the sky source/component. Now, let

$$\mathbf{Z}_{ip,nk''} = \sum_k \mathbf{A}_{ipkn} \mathbf{F}_{in} \mathbf{M}_{kk''n}, \quad (10)$$

be a matrix (note the comma between ip and nk'' in the indices, indicating that we are rolling two indices into one axis). Then, we simply have that

$$V_{ij}^{pq}(v_a, t) = \mathbf{Z}\mathbf{Z}^\dagger, \quad (11)$$

i.e. the full visibility can be expressed as a matrix product of a $N_{feed} N_{ant} \times N_{ax} N_{src}$ matrix with its own conjugate transpose. This matrix product produces the necessary two beam factors (with the Hermitian conjugate), squares the square root of the coherency matrix, and forms the necessary phase factor.

The required memory for this matrix product is quite reasonable, as there is simply one copy of an array that is of size $4N_{ant}N_{source}$, and matrix-matrix products are extremely fast with modern algorithms and architectures, especially on the GPU.

We may go a little further if we assume that the sky is unpolarized. In this case, the coherency matrix can be approximated as

$$\mathbf{C}(v_a, t) = \begin{pmatrix} I(v_a) & 0 \\ 0 & I(v_a) \end{pmatrix}, \quad (12)$$
in which case we simply have

$$ \mathbf{M} = \begin{pmatrix} \sqrt{I(v_a)} & 0 \\ 0 & \sqrt{I(v_a)} \end{pmatrix}. \qquad (13) $$

Having zero off-diagonal terms means that the sum over k in equation (10) is no longer required, and we can simply write

$$ \mathbf{Z}_{ip,nk''} = \mathbf{A}_{ipkn} \mathbf{F}_{in} \sqrt{I_n}, \qquad (14) $$

which saves a rather large sum. The combination of equations (11) and (14) is what matvis currently implement in its algorithm.

There are a few important observations to be made about this final approximation. First, we again note that we have assumed that all sky sources are unpolarized. This approximation has served well for the initial intended purpose of the matvis software-producing simulations for the validation of the HERA analysis pipeline. We note that even though the sky is approximated as unpolarized, the output visibilities still have non-zero cross-polarization, due to the inclusion of polarized beam sensitivity on each feed. Thus far, this unpolarized formulation is the only one implemented in the matvis software (driven by these motivations), but it will be trivial to add full polarization support according to equation (9) in future versions.

Secondly, note that the peak memory usage for this algorithm is almost completely controlled by the number of antennas and feeds, as they are combined as an outer product. Both frequency and time are neither summed over nor interact with other frequencies and times, making them ideal as outer-loops for the algorithm. The source and field-component axes are summed over, and therefore may be 'chunked' to suit any memory needs, at the expense of some vectorization. Thirdly, we re-iterate that the formulation of equation (11) requires that all baselines are simulated at once. This can be a disadvantage when simulating arrays with high number of redundancies although the optimization implemented here still leads to several orders of magnitude speed up in comparison to a simulator that directly evaluates per-baseline visibilities in equation (3) such as pyuvsim (see Section 6). Nevertheless, we believe that further optimization to this algorithm is possible, and we describe the idea in Appendix A to be fleshed out in future work.

## 3.3 Coordinate transformations

One of the key contributors to the accuracy of a visibility simulator, as well as its computational performance, is the manner in which the necessary coordinate transformations are carried out. On the one hand, there is a sky model, which is most naturally represented in a celestial coordinate system, e.g. equatorial coordinates. On the other, there are the instrumental terms, particularly the primary beam response and the set of baseline vectors, which for stationary antennas like HERA, the MWA, and SKA, are most naturally represented in a horizontal coordinate system connected to the location of the array. The complex exponential (fringe pattern) term in the RIME mixes the two, as it involves both the source position and baseline vectors. It is therefore always necessary to transform between celestial and horizontal coordinates at some point in the RIME calculation.

One way to perform these coordinate transforms is to use simplified trigonometric calculations that (i) converts from equatorial (RA and Dec, α and δ) coordinates to a Cartesian system via geometric vector projection, and then (ii) performs a solid body rotation into a horizontal coordinate system. The first step is,

$$ X_{eq}^{(n)} = (\cos \alpha^{(n)} \cos \delta^{(n)}, \sin \alpha^{(n)} \cos \delta^{(n)}, \sin \delta^{(n)}). \qquad (15) $$

This vector of coordinates only needs to be calculated once for each source n in the catalogue. The rotation matrix for (ii) is the same for

all point sources for a given LST, and so there is also a computational saving from only needing to calculate it once per LST. Defining the horizontal coordinate system as an ENU (East-North-Up) topocentric system, we obtain the rotation matrix,

$$ \mathbf{R}(t) = \begin{pmatrix} \sin(H) & \cos(H) & 0 \\ -\sin(l) \cos(H) & \sin(l) \sin(H) & \cos(l) \\ \cos(l) \cos(H) & -\cos(l) \sin(H) & \sin(l) \end{pmatrix}, \qquad (16) $$

where H is the hour angle of a reference point on the celestial sphere as a function of LST, and l is the latitude of the array on the Earth. Then, the horizontal coordinates of the sources are simply,

$$ X_{hrz}^{(n)}(t) = \mathbf{R}(t) \cdot X_{eq}^{(n)}. \qquad (17) $$

Although this is fast, as it relies only on basic operations, and can be used to give a reasonable representation of the necessary coordinate transformations, it is also inaccurate owing to the lack of astrometric corrections. A number of time-dependent corrections for precession, aberration, and so on, are generally required to recover the apparent positions of celestial sources within arcsecond accuracy. In the case of matvis, the fringe term, i.e. the exponent term in equation (14) which requires time-dependent source positions in horizontal coordinates, is particularly sensitive to positional errors, and we have found in our tests that the accuracy of the equatorial-to-horizontal coordinate transform is often the main driver of the precision of the simulator as a result.

The easy solution is to use astropy as a full set of astrometric corrections has been implemented in its coordinate transformation module. However, these corrections can be expensive to compute, particularly for large catalogues containing many sources. During the early state of the development of matvis, we have found that employing a full astropy-based coordinate transforms would result in the coordinate transformation step being the bottleneck, significantly limiting the performance of matvis.<fnref n="5" />

As such, we implement a two-step coordinate transformation scheme, which in our case are applied to the point sources, to find their locations in Cartesian horizontal coordinates at each LST. In short, we perform a pre-computation procedure to 'perturb' the input equatorial coordinates of the sources so that steps (i) and (ii) in the trigonometric transforms (i.e. equations 15 and 16) result in a more accurate observed source location when applied. These fictitious coordinates only need to be computed once, and significantly improve the accuracy of the results. Algorithm 1 details how we compute them.

When steps (i) and (ii) are subsequently applied to the perturbed equatorial coordinates for other observation times in the simulation, we find a significant improvement to the accuracy of the observed source locations compared with the unperturbed ones. This correction degrades with increasing separation in time from the reference time, although it will improve again once the separation is longer than ≈ 12 h, and reach the best improvement again at slightly less

<footnote n="5">⁵During the course of the publication of this paper, we have significantly generalized the implementation of sky-coordinate calculations within matvis. In version from 1.3.0 onwards, the user is not restricted to using rigid rotations for the sky coordinates, but instead any coordinate transform may be used. In particular, the accurate transforms from astropy are built-in, along with new GPU implementations of the same, as well as schemes similar to the fast rigid-rotators with reference-time corrections described in this paper. The coordinate transformation schemes discussed in this paper are still applicable to version 1.2.1 and below, and have been used for all tests in Section 5. The details of the improvements mentioned here will be discussed in a future publication.</footnote>
**Algorithm 1** Algorithm for matvis’s ‘fast’ coordinate transformation scheme.

(i) Given a reference time $t_{\text{ref}}$, which can be different that the observation time of the simulation, compute the apparent local sidereal time at the location of the telescope via *astropy*.

(ii) Compute $\mathbf{R}'(t_{\text{ref}})$, the ENU-to-topocentric rotation matrix (inverse of equation 16) at the reference time with the LST obtained from the above step and the latitude of the telescope.

(iii) Convert the input equatorial source positions to horizontal coordinates in altitude and azimuth, $\theta_{\text{hrz}}$ through the *astropy*.s SkyCoord transform to method, using $t_{\text{ref}}$ as the observation time.

(iv) Perform a geometric vector projection of $\theta_{\text{hrz}}$ into Cartesian coordinate to obtain $\mathbf{X}_{\text{hrz}}^{(n)}(t_{\text{ref}})$, essentially an equivalent of equation 15.

$$ (v) \text{ Compute } \mathbf{X}_{\text{eq,perturbed}}^{(n)} = \mathbf{R}'(t_{\text{ref}}) \cdot \mathbf{X}_{\text{hrz}}^{(n)}(t_{\text{ref}}) $$

(vi) Do the inverse of equation 15 to infer the RA and Dec from $\mathbf{X}_{\text{eq,perturbed}}^{(n)}$.

than 24 h from the reference time, because the correction is tied to the source positions. This also means that the source positions will match those from the *astropy* transformation and thus errors from source will be minimum when simulating a snapshot observation at a single time. We show how the errors from this correction can affect our use cases of *matvis* in Section 5.

## 3.4 Beam interpolation

One of the most important factors in the accuracy of the simulator is how accurately the antenna beams are represented. This is made more complicated by the fact that the beams must be evaluated at floating sky-coordinates defined by the discrete sky sources, as they drift in time. *matvis* allows two main forms of antenna beam to be input: analytic beams, and simulated beams. The former are easily evaluated at any coordinate (and frequency) by their nature, and typically in a very efficient manner. The latter presents more challenges; simulated beams are typically simulated on a discrete angular grid that cannot coincide with the floating sky-coordinates at which the beam must be evaluated. Therefore, they must be interpolated.

*matvis* makes use of the UVBeam class from the `pyuvdata` package (Hazelton et al. 2017), developed and maintained by the Radio Astronomy Software Group (RASG<fnref n="6" />), to handle the beam information and interpolation of regularly gridded, HEALPIX gridded, and CST-simulated beams. The UVBeam class also implements a number of interpolation methods (with more being added over time), which *matvis* directly exploits in its CPU implementation. The order of interpolation in UVBeam is selectable, and can use either ` numpy摘自` `RectBivariateSpline` or the `map出口` `coordinates` 出现从 `scipy` 内联。For example, this can be bilinear, bicubic, or quintic, on a regular grid of altitude and azimuth. Thus, while accurate beam interpolation is important for accurate simulations, *matvis* out-sources this complexity, and arbitrarily accurate and/or performant interpolation schemes can be adopted.

This picture is not quite as simple for the GPU implementation. There are currently no GPU-accelerated beam interpolators within the UVBeam class, so beam interpolation on the GPU is implemented directly within *matvis*. The simplest such interpolation is bilinear, and thus far this is the only scheme implemented for the GPU version of *matvis*.

Bilinear interpolation is likely to be more approximate than bicubic or higher orders of interpolation, but we note this limitation is *not* fundamental to the *matvis* algorithm – it is simply a practical limitation of the current version of the code. It will not be overly difficult to add higher order interpolation schemes in the future. In any case, the degradation in accuracy incurred by bilinear interpolation is strongly dependent on the inherent structure in the beam model (i.e. more structure at smaller scales will incur higher inaccuracy) and the resolution of the raw beam input. The latter is controlled by the user to some degree – certainly for particular experiments, it is possible to re-simulate the beam with an EM solver (e.g. FEKO or CST) at higher resolution if required. If this is not possible for practical reasons, the beam can simply be interpolated to a higher regular resolution grid with higher order interpolation before simulation with *matvis*, and then the bilinear interpolation applied to that ‘upsampled’ beam model. Thus, even the particular practical limitation of the current GPU incarnation of *matvis* is not *really* a limitation but can be overcome by pre-processing.

Regardless, for the particular purpose of validating the HERA analysis pipeline, which was the primary motivation for development of *matvis*, we find that bilinear interpolation on our CST-modelled beam with resolution of 1 deg² (Fagnoni et al. 2021) is suitable as it does not hinder the power spectrum derived from the simulations. We will discuss this in Section 5.

One point to note concerning this interpolation is that it is essential that the regular coordinates of the simulated lattice be defined on a spherical grid which has the zenith as a pole. This is due to the fact that many common E-field beam responses have a discontinuity in their phase (recalling that the response is complex) for one of their components at zenith, which is properly interpolated only if the zenith is considered a pole and interpolated ‘away from’.

We describe how the beam information (as well as other required inputs) are passed to *matvis* in the next section.

## 3.5 Summary of algorithm

Having described the mathematical framework in the preceding sections, we now summarize the algorithm of *matvis* as it stands in its present implementation.<fnref n="7" />

First, we note that each step is to be performed at a single frequency. Thus, frequency forms our outer-most loop, and indeed it is envisaged that each frequency would typically be performed as a separate job on a high-performance computing scheduler, i.e. as an *embarrassingly parallel* task. Our second, internal, loop is over times.

The following inputs are expected by the basic *matvis* algorithm:

(i) A set of antenna locations, $\mathbf{D} = \{\mathbf{x}_1, \dots, \mathbf{x}_{N_{\text{ant}}}\}$ in Cartesian East-North-Up (ENU) coordinates as an $N_{\text{ant}} \times 3$ matrix.

(ii) A beam model, $A_{ijk}(v, \theta)$ for each *unique* antenna, $i$, feed $j$, and E-field component $k \in (\theta, \phi)$, that may be evaluated (or interpolated) to any particular set of topocentric coordinates. Note that two physically separate antennas that nevertheless are assumed to have the same beam need only be provided once.

(iii) An integer-valued vector $a_i$ whose value for antenna $i$ is the index of its corresponding beam model, $i$.

(iv) A matrix of sky model source locations in Cartesian equatorial coordinates (ECI), $\mathbf{X}_{\text{eq}} = \{\mathbf{X}_{\text{eq}}^{(0)}, \dots, \mathbf{X}_{\text{eq}}^{(N_{\text{src}})}\} \in \mathbb{R}^{3 \times N_{\text{src}}}$, where each source is given by equation (15).

(v) A positive, real-valued vector of source intensities, $\mathbf{I} \in \mathbb{R}_{+}^{N_{\text{src}}}$, where $\mathbf{I}_n = 2\mathbf{C}_{00}^{(n)} = 2\mathbf{C}_{11}^{(n)}$.

<footnote n="6">⁶https://radioastronomysoftwaregroup.github.io/</footnote>

<footnote n="7">⁷That is version 1.2.1 at the time of the publication of this article.</footnote>
(vi) An array **R** of shape $N_{\text{times}} \times 3 \times 3$ containing the rotation-matrices (cf. equation 16) required to convert from equatorial to topocentric coordinates at each time-step. These are computed *outside* the main algorithm loop using any method the user desires (though utilities for computing these are provided in the package). These matrices also *define* the times of the simulation (no other time information is given to the algorithm).

Then, for a particular frequency *v* and time-index *t*, the *matvis* algorithm follows Algorithm 2.

**Algorithm 2** Main algorithm for matvis.

(i) Rotate the pixels/sources above the horizon into topocentric frame: $\mathbf{X}_{\text{hrz}} = \mathbf{R}_{t}^{\top} \mathbf{X}_{\text{eq}}$.

(ii) Interpolate the beams onto source coordinates, i.e. produce the $N_{\text{beam}} N_{\text{feed}} \times N_{\text{ax}} N_{\text{src}}$ complex-valued matrix $\mathbf{A}_{i,j,kl} = A_{i,j,k}(v, \mathbf{X}_{\text{hrz}}^{(l)})$ for unique antenna beam $i$, feed $j$, E-field component $k$ and source $l$.

(iii) Compute the antenna-based exponent: $\tau = -2\pi i v \mathbf{D} \mathbf{X}_{\text{hrz}} / c$, which is an $N_{\text{ant}} \times N_{\text{src}}$ matrix.

(iv) Compute the $N_{\text{ant}} N_{\text{feed}} \times N_{\text{ax}} N_{\text{src}}$ per-antenna complex-valued matrix $\mathbf{Z}_{i,j,kl} = \sqrt{\bar{\mathbf{I}}_i} \mathbf{A}_{d_i j, kl} \exp(\tau_{il})$, the matrix form of equation 14.

(v) Compute the $N_{\text{ant}} N_{\text{feed}} \times N_{\text{ant}} N_{\text{feed}}$ visibility: $\mathbf{V} = \mathbf{ZZ}^{\dagger}$ (equation 11).

There are a few important details to note with respect to this algorithm. First, the present implementation of this algorithm only allows simulating drift-scanning observations, which is the common observing mode of 21 cm cosmology experiments such as HERA. This limitation is due to the fact that the algorithm simply evolves the input source matrix by making a horizontal cut and then computing the visibilities from all sources that are rising for each iteration in time. To enable simulating tracking observations, an additional step that adds time-dependent phase to each antenna would need to be implemented. This would have only a small performance impact, but is currently not in the scope of the *matvis* development team.

Secondly, since steps (iv) and (v) are simply matrix products, highly efficient thread-optimized implementations exist to compute them on CPUs, e.g. in BLAS.<fnref n="8" /> Furthermore, drop-in replacements exist to compute matrix products even more efficiently on GPUs, for example the CUBLAS package. This is highly advantageous, as step (vi) is the computation bottleneck for realistic modern large-*N* arrays.

Lastly, the beam interpolation in step (ii) need only be performed for the number of *unique* antenna beams. In many cases of current interest, all antennas will be assumed to have the same beam, which makes this step negligible in terms of computation. Nevertheless, if the simulation is intended to capture beam irregularities across the antennas, this step can become highly non-negligible. However, GPUs are also extremely well-adapted to spatial interpolation; thus a well-tuned GPU implementation, e.g. by using texture mapping, can significantly reduce this cost.

# 4 THE matvis PACKAGE

The *matvis* package defines an API and an algorithm for computing the RIME, with the intention that different implementations can

be provided as drop-in replacements of each other. Currently, the package provides two such implementations: *matvis.cpu* and *matvis.gpu*, which define CPU-based and GPU-based implementations, respectively. Notably, the API computes the RIME at a single frequency only but includes a loop over time. The majority of the code (all of the CPU-based implementation) are written in pure Python, utilizing the standard *numpy*, *scipy*, and *astropy* libraries. The GPU-based implementation replaces some of the *numpy* operations, particularly steps (ii)–(v) (as described in Section 3.5), with CUDA11 code wrapped by PyCUDA. Each step of the algorithm is implemented as its own CUDA kernel, glued together in Python to maintain code simplicity. Notably, matrix products directly use the highly optimized CUBLAS library.

To maintain feature parity between the CPU-based and GPU-based implementations, the coordinate transformation scheme outlined in Section 3.3 is separately implemented in the function `equatorial_to_eci_coors` in the `conversions` module of the Python package because a GPU-based `astropy` coordinate transformation module does not exist. Users can manually call this function to put sources into the 'perturbed' equatorial coordinates (as described in Section 3.3) before passing them to the main *matvis* routine to improve accuracy of the coordinate rotation in step (i). Although this is not required, i.e. the API has no knowledge of whether this correction has been applied, we highly advise that all simulations use it (see Section 5 on how this influence the simulation outputs).

As discussed, interfacing and interpolation of gridded beams are handled by the UVBeam from the `pyuvdata` package developed and maintained by RASG. For analytic beams, uniform, Gaussian, and Airy beams are readily supported through the *AnalyticBeam* class of the `pyuvsim` package, which is also developed and maintained by RASG, and shares a common beam evaluation interface with the UVBeam interpolation. Other forms of analytic beams are possible by sub-classing the *AnalyticBeam* class, one of which that we have developed is the *PolyBeam* class in the *hera_sim* package, which implements an analytical model of the HERA Phase I primary beam (Choudhuri, Bull & Garsden 2021, see Section 4).

A high-level wrapper is available via the `simulate{vis}r routine, which provides a loop over frequencies, as well as simple implementations to compute the required input equatorial sky coordinates, $\mathbf{X}_{\text{eq}}$ from the more common RA/DEC. This high-level wrapper also provides a simple switch between different back-end implementations (i.e. CPU and GPU). We note that this wrapper does *not* provide in-built capability to perform the coordinate perturbations (though perturbed coordinates can be computed and provided to the wrapper by using the above mentioned function).

A separate wrapper is also available via the * visiblities* module in the *hera_sim* package, a simulation package for HERA-like redundant interferometric arrays, also developed by the HERA collaboration. This has a few key advantages over the in-built wrapper in the *matvis* package, including an option to apply coordinate correction given a reference time, a unified simulation interface that can be used to wrap other simulator codes (with *pyuvsim* wrapper already implemented), and integration with systematics models available in the *hera_sim* package. The main script of this wrapper, *hera-sim-vis.py*, uses *pyuvsim*-style configuration files,<fnref n="9" /> making it much easier to configure and execute a large set of simulations on an HPC cluster. The *hera_sim* package also implements an analytical model of the Phase I HERA dipole beam (Choudhuri et al. 2021) in its *PolyBeam* class for use with

<footnote n="8">⁸Basic Linear Algebra Subroutines: https://www.netlib.org/blas/.</footnote>

<footnote n="9">⁹See https://pyuvsim.readthedocs.io/en/latest/parameter_files.html</footnote>
matvis (and pyuvsim) although we note that the pyuvsim-style
configuration files do not currently support this beam (or any analytic
beam models implemented outside of the AnalyticBeam class of
the pyuvsim package), and thus users must manually build matvis
simulations via the hera_sim unified interface to use this beam.

The matvis package<fnref n="10" /> and all of its associated software that we
have mentioned in this paper, including hera_sim, <fnref n="11" /> pyuvsim, <fnref n="12" />
and pyuvdata, <fnref n="13" /> are publicly available via open-source licenses
on GitHub. PyPI distributions are also available, allowing easy
installation via the pip install command. matvis has been
under active development since May 2021. The development cycle
employs git and GitHub and strictly follows an open-source devel-
opment guideline laid out by the HERA collaboration. <fnref n="14" /> Particularly,
we adopt continuous development and comprehensive code review
processes, thorough unit-testing, including for pyuvsim and pyu-
vdata, and maintain detailed and user-friendly documentation, <fnref n="15" />
including both API-level ‘docstrings’ and a set of tutorials. We
use continuous integration (CI) to run tests and report coverage.
Although the primary contributions to the matvis code have thus
far come from members of the HERA collaboration, we welcome
and appreciate contributions from outside of the collaboration.

5 VALIDATION

In this section, we perform several tests to demonstrate suitability
of matvis simulations for validation of 21 cm power spectrum
analyses. There are two important metrics when quantifying the
accuracy of mock visibilities, and thus the visibility simulators, for
a 21 cm power spectrum analysis. The first is that the dynamic range
of the delay spectrum of foreground visibilities must be > 10¹⁰ to
allow discerning the foreground and 21 cm signals at high delay.
Secondly, the simulated 21 cm visibilities must yield estimated power
spectra that accurately match the input 21 cm power spectrum. We
will show that matvis meets these two criteria in the following
two subsections. Throughout these tests, we will also compare the
accuracy of matvis against a reference implementation of the
RIME in the pyuvsim simulator, which we discussed in Section 2.

In addition to testing the simulator outputs on the targeted pipeline
validation use case, we elaborate on the impact of the point-source
approximate used by matvis (as well as pyuvsim and some
other simulators), particularly with regard to the discretization of
the diffuse sky emission. We will provide a demonstration in the
context of the delay spectrum analysis. We also touch on the imaging
of simulated visibilities from matvis, with a focus on the impact
of matvis' 'fast' coordinate transformation scheme on the point
source position in the image.

With these tests, we aim to demonstrate the accuracy that matvis
can achieve on the type of experiments that it is designed to simulate
and for the targeted scientific measurements that those experiments
aim to measure. However, we also note that the results presented here
should not be treated as an absolute measure of its accuracy. The
numbers will change depending on the simulation parameters and
sky models, and whether they are sufficiently accurate will depend
on specific use cases.

The following software versions are used throughout these tests:
matvis==1.2.1, hera_sim==4.2.0, pyuvsim==1.2.6,
and pyuvdata==2.4.0. For all simulations, we use the CST-
simulated HERA beam from Fagnoni et al. (2021) with bi-linear
and bi-cubic interpolations in spatial and frequency dimensions,
respectively.

5.1 Delay spectrum dynamic range

One of the key concerns in 21 cm power spectrum analysis is the
dynamic range of the foreground signal. This is especially important
for the foreground-avoidance delay-spectrum technique used by
HERA (Parsons et al. 2012). The basis of this measurement is that
the clean and calibrated visibility is Fourier transformed in frequency
with a chosen taper function before forming a power spectrum
over the delay modes. The foreground power will be smooth and
concentrate around zero delay, forming a peak, and taper off at high
delay, if there are no other systematics. The 21 cm delay spectrum,
on the other hand, will be relatively flat, unsmooth, and lies ~ 4 – 5
orders of magnitude lower than the foreground power at zero delay.
Thus, the dynamic range between the peak and the floor of the
foreground delay spectrum is an important factor to allow discerning
of the foreground and 21 cm signals.

To demonstrate that matvis is able to produce simulations
that yield delay spectra of sufficient dynamic range, we perform
simulations with the same sky models used in the validation of the
recent 21 cm power spectrum upper limits from HERA Collaboration
(2023), and derive delay spectra from the simulation outputs. These
sky models consists of a diffuse Galactic synchrotron emission model
taken from the Global Sky Model (de Oliveira-Costa et al. 2008)
represented as a HEALPix map with Nside = 256, a full-sky mock
point source catalogue derived from source count distribution of the
GLEAM catalogue (Franzen et al. 2019), and the extremely bright
‘A-team’ sources extracted from the GLEAM catalogue (Hurley-
Walker et al. 2017). We perform the simulations at two times, LST
= 22.75 h (Field A) and LST = 5.13 h (Field C), each with 5-s
integration over two frequency bands, 117.09–132.62 MHz (Band 1)
and 150.29–167.77 MHz (Band 2). These fields and frequency bands
correspond to those used for measurements in HERA Collaboration
(2023). We simulate two different baseline lengths in East–West and
North–South direction, which correspond to the shortest and longest
baselines in those directions in the current HERA array. We calculate
the delay spectrum as discussed above with a Blackman-Harris
taper function. Source position correction is applied as discussed
in Section 3.3, for which we use the simulated times as the reference
times for coordinate transformations. This is essentially equivalent
to no approximation on the source positions.

Fig. 1 shows the results from this test. The top row shows delay
spectra derived from both simulators at different baseline lengths
and directions within the core of the HERA array. The dashed
line shows delay spectra of the Blackman-Harris taper function,
providing a baseline of how well one can possibly do in terms
of dynamic range. We also peak-normalized all power spectra in
the top row to make it easier to compare their dynamic ranges.
The bottom row shows the fractional residuals between matvis
and pyuvsim, $|V|_{\text{residual}}^2 = |V|_{\text{matvis}}^2/|V|_{\text{pyuvsim}}^2 - 1$. matvis
produces delay spectra with dynamic ranges of more than 10 orders
of magnitude in all cases; the cases with the lowest dynamic range
are the longest baselines, which likely suffer some aliasing of the
fringe from the diffuse sky model (indicating a deficiency in the
input model, rather than the simulator itself). This interpretation of a
sky-model deficiency is supported by the exceptional agreement of

<footnote n="10">¹⁰https://github.com/HERA-Team/matvis</footnote>

<footnote n="11">¹¹https://github.com/HERA-Team/hera_sim</footnote>

<footnote n="12">¹²https://github.com/RadioAstronomySoftwareGroup/pyuvsim</footnote>

<footnote n="13">¹³https://github.com/RadioAstronomySoftwareGroup/pyuvdata</footnote>

<footnote n="14">¹⁴http://reionization.org/wp-content/uploads/2013/03/HERA067_HERA_Software.Community_Guide_6May2019.pdf</footnote>

<footnote n="15">¹⁵https://matvis.readthedocs.io</footnote>
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Delay spectra from *matvis* and *pyuvsim* simulations with mock foreground models at two different frequency bands and times at different baseline lengths and directions (top), and the fractional residuals between the results from the two simulators (bottom). The dashed line shows the delay spectrum of the Blackman-Harris taper function, which indicates the maximum possible dynamic range. All power spectra are peak-normalized to unity to emphasize their dynamic range. *matvis* simulations yield delay spectra with dynamic ranges of > 10¹⁰ in all cases, and agree with *pyuvsim* to better than 10⁻¹⁰ per cent.</figcaption></figure>

*matvis* and *pyuvsim* – at better than a millionth of a percent – for all times, frequencies and baseline lengths. Given the extensive testing of *pyuvsim* (including against analytic solutions; Lanman et al. 2022), this suggests that – beyond simply supporting power-spectrum validation efforts – *matvis* could also be used for detailed calibration purposes, so long as the input sky models are sufficiently detailed, since the rule-of-thumb precision required for such purposes is about 10⁻¹⁰ in power-spectrum units (Barry et al. 2016).

Agreements with *pyuvsim* will degrade if the reference time used for the coordinate transformation differs from the simulated time, for reasons that we discussed in Section 3.3. Fig. 2 demonstrates this point. Here, the top row shows a heat map of delay spectra from *matvis* simulations of field A in frequency band 1. In these maps, the vertical axis captures the offset between the reference and simulation times, and we simulate every 15 min. The bottom row shows their fractional residuals against *pyuvsim*. As expected, *matvis* agrees best with *pyuvsim* when simulating directly at the reference time of the coordinate transformation. When the simulated time differs from the reference time, discrepancies between the two simulators increase rather quickly, with $|V|_{\text{residual}}^2$ doubling when the reference time is ≈ ±1 h from the simulated time. However, the dynamic ranges of the power spectra do not degrade as the simulations move away from the reference time, as illustrated by the vertical constancy of the maps in the top panel. Here, the dashed white lines indicate the horizon delay (Parsons et al. 2012), and we note that at all times, the delay spectra achieve the required 10–12 orders of magnitude of dynamic range at or around the horizon. Thus, we consider the simulations sufficient for the purposes of validating power-spectrum analyses (for which the dynamic range and overall magnitude are the foremost concerns), but too discrepant at large offsets from the reference time for the purposes of calibration. For the latter, we suggest splitting the simulation into chunks of 15 min or less, and performing the coordinate correction (c.f. Section 3.3) at a central reference time for each chunk.

## 5.2 21 cm power spectrum

The ability to obtain an accurate power spectrum estimate from the visibility simulated from a known 21 cm model is also an important factor that goes hand-in-hand with the delay spectrum dynamic range for a visibility simulator to be feasible for 21 cm power spectrum analysis. We test this by simulating from a 21-cm EoR analogue model similar to the one used in HERA Collaboration (2023). The model is a Gaussian random field with a power spectrum slope, $P(k) \propto k^{-2.7}$, created using the *gaussian_random_fields*<fnref n="16" /> code. The EoR Gaussian field is drawn at pixel locations of a set of HEALPIX maps with NSIDE=1024, each map corresponding to a different redshift and simulation frequency. We then discretize this set of maps into an ensemble of point sources as required by *matvis* by multiplying the intensity of each pixel (in Jy sr⁻¹) by the pixel area, converting each pixel into a 'point source' with flux density in Jy.

In general, we expect realistic sky models to be strictly positive-valued, because foregrounds are the dominant component. However, the 21 cm signal is a differential quantity, and can be measured in absorption or emission. When simulating the 21 cm signal in isolation from the foregrounds, we thus expect to encounter negative values, which are a problem for the *matvis* simulator as it is currently implemented. To deal with the negative values for this test, we add a constant value to the entire 21 cm sky model (all pixels at all frequencies) such that it becomes strictly positive. This is equivalent to adding a monopole and will result in a bump in the power spectrum at low *k*, but should not modify the estimated power at larger *k*. This monopole bump can be removed by simulating just the monopole component separately and subtracting the resulting visibilities from the monopole-added visibilities.

For the purposes of this test, we run *matvis* simulations using a small HERA-like configuration: 19 antennas (hexagonal

<footnote n="16">¹⁶https://github.com/zacharymartinot/redshifted_gaussian_fields</footnote>
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Foreground delay spectra from *matvis* simulations with different coordinate transformation reference time (top) and fractional residuals against *pyuvsim* (bottom). It is evident from the residuals that setting the reference time further away from the simulated observation time will result in errors, the larger the differences, the higher the errors. If accuracy is required, e.g. using *matvis* simulation for calibration time, we recommend splitting the simulation into < 2-h time chunk, or simulate one time-step at a time for best results. These errors, however, do not affect the delay spectrum dynamic range as also evident in the top row.</figcaption></figure>

arrangement with 3 antennas on each side), 1 snapshot integration at LST = 0 h, and 131 frequency channels spanning the Band 1 frequencies (117.09–132.62 MHz).

We form the power spectra with the *hera_psepc*<fnref n="17" /> code. This estimate is a complicated matter, and we suggest interested readers to consult, e.g., HERA Collaboration (2023) or Gorce et al. (2023) for details. The top panel of Fig. 3 shows the resulting power spectrum estimates. The input power spectrum, corrected for expected aliasing and instrumental window-function effects (Aguirre et al. 2022), is shown in blue. The red crosses are the baseline-averaged power spectrum from the *matvis* visibility without monopole correction, which exhibit a prominent bump at *k* < 0.20. Although this bump seems problematic, it falls beneath the foreground power in those *k* modes, and thus will not be an issue if the visibility is used for validating foreground avoidance analysis used by HERA. Regardless, it is simple to correct for the monopole (as discussed above), and the estimates with this monopole correction are shown as black circles. These follow the expected power spectrum well for all *k*. Individual-baseline estimates with monopole corrections are shown as grey circles for completeness.

Discrepancies between the estimates and the input power spectrum are expected to arise simply due to sample (sometimes also called ‘cosmic’) variance. The expected sample variance is a decreasing function both of the total integration time and the number of unique baselines at each *k* (Lanman & Pober 2019). We use the sample variance prediction from Lanman & Pober (2019) as integration

time approaches zero, calculated for the HERA beam (their fig. 5), to estimate the sample variance we expect. For a single baseline, the expected RMS fluctuations from sample variance are 100 per cent, and we divide by $\sqrt{N_{bl}} \approx 5.5$, yielding 18 per cent RMS (red dotted line in bottom panel of Fig. 3). The residuals of the estimated power spectra to the input are generally within the expected sample variance.

## 5.3 Simulating from diffuse sources

As discussed extensively already, *matvis* assumes a point-source approximation for a sky model. To simulate visibilities from a resolved radio galaxy with extended jets or diffuse emission, which include the 21 cm signal, the user must first discretize the sky intensity into point sources, essentially a collection of delta functions that are assigned the integrated flux density of surrounding regions. *matvis* is agnostic to the specific choice of this discretization, i.e. the location of each point source and the shape of the region surrounding it, within the dual limitations that (i) the sum of all discretized sources should yield the integrated flux density of the full sky, i.e. the discretization should be a complete tessellation of the sky, and (ii) the fidelity of the simulation will, in general, be higher when the pixels are more compact (e.g. long ‘stripes’ as pixels will be a poorer discretization – for the same number of sources – than pixels with small ‘circular’ surrounding regions). For diffuse emission, one choice of explicit discretization that is consistent with these assumptions is the HEALPIX coordinate grid, by which each of the HEALPIX pixel is treated as a point source. However, the calculation of the RIME is expected to break down on long baselines for a lower resolution HEALPIX grid

<footnote n="17">¹⁷https://github.com/HERA-Team/hera_pspec</footnote>
Redshift $z = 10.37$

<figure><img src="image_5.png" /><figcaption>**Figure 3.** 21 cm power spectrum estimates from *matvis* simulations. In the top panel, the solid blue lines show the expected power spectrum of input model used in the simulation. This model has some negative values, and thus the a monopole has been added by shifting values to be all positive to conform with *matvis* formalism. The power spectrum estimates from the resulting visibility shown in red crosses exhibit a monopole bump at low *k*. This bump is corrected by separately simulating the monopole and subtracting it before estimating the power spectrum. The estimates with this correction are shown in black circles, which follows the expected power spectrum at all *k*. Residuals between the expected and estimated power spectra are primarily due to cosmic variance, for which we show the level predicted by Lanman & Pober (2019) for HERA as the red dotted lines in the lower panel.</figcaption></figure>

because the fringe terms (the exponent in the RIME), which depends on the source positions and baseline lengths, would not be evaluated enough to sufficiently sample the fine fringe patterns of the long baselines.

To demonstrate, we compare the delay spectra of simulation outputs of solely the diffuse component for Band 1 (117.09–132.62 MHz) and Field C from Section 5.1, which was discretized on a HEALPIX grid with NSIDE=256, to those derived from simulations from the same diffuse models discretized on lower resolution HEALPIX grids. Fig. 4 shows these results, for which the top row is the derived delay spectra, and the bottom row is the fractional residuals between the spectra derived from a lower resolution sky model to that of the NSIDE=256 model. The NSIDE parameters of 256, 128, and 64 roughly correspond to angular resolutions of approximately<fnref n="18" /> 13.74, 27.48, and 54.96 arcmin.

First, we note that at 132.63 MHz, the fringe spacing is ≈ 34.5 arcmin for the 226-m EW baseline, the longest baseline in our simulation, whereas the approximate resolution of the HEALPIX grid

at NSIDE=256 is 13.74 arcmin. This means that the fringe patterns are sampled at a slightly higher rate than the Nyquist rate, and thus the discretized diffuse model should capture most of the information in the original GSM map, which should translate into the most accurate delay spectrum for the baselines that we simulate. We have run a test with a diffuse model of higher NSIDE and found it to be in numerical agreement with the NSIDE=256 case.

As we lower the resolution of the diffuse model, the discretized source density reduces, as well as the sampling of the fringes, resulting in aliasing and leading to errors in the power spectrum. This is apparent in the case of NSIDE 64 and 128 for 226-m EW and 143-m NS baselines. The latter corresponds to a fringe spacing of ≈ 54.6 arcmin, and thus provides a close to Nyquist sampling, but, as we mentioned, the HEALPIX pixel resolution is not exact, and this leads to small errors in the delay spectrum as shown.

For the 14-m baseline, the fringe spacing is ≈ 5°, which should be coarse enough to be well sampled by all HEALPIX resolutions in this test, and the resulting delay spectrum matches visually at all NSIDEs. However, their fractional residuals are still ~ 10⁻² or 1 per cent at all delays. This may or may not be accurate enough depending on the use case. In general, a good rule of thumb is to sample the sky model beyond the Nyquist rate.

<footnote n="18">¹⁸Since the shape of the HEALPIX pixels are not all the same, this number is approximated.</footnote>
<figure><img src="image_6.png" /><figcaption>**Figure 4.** Demonstration of the effects of insufficient resolution in the discretization of a diffuse sky model on the delay spectrum. We simulate a diffuse sky model discretized onto a HEALPIX grid at different resolutions defined by the NSIDE parameter, and then form the delay spectra from the simulated output. At a low sky model resolution, the fringe term, which depends on the source positions, may not be evaluated sufficiently for long baselines (e.g. second and fourth panels) with fringe spacing shorter than the approximate angular resolution of the discretization. This aliasing results in errors in the visibility that propagate to the delay spectrum.</figcaption></figure>

## 5.4 Imaging matvis visibilities

Although *matvis* was primarily motivated by the need to validate 21 cm drift-scan-style delay-spectrum experiments, and currently does not support simulating tracking observations, imaging the output visibilities is certainly possible. In this case, the main concern is how errors from our fast coordinate rotation scheme manifest in the resulting images.

We simulate multiple 5-s snapshots of a single, 100-Jy, point source at RA= 12 h and Dec= 33.6° transiting near the zenith of the full, 350-antenna, HERA array (maximum baseline ≈ 800 m). In each case, the snapshot time remains the same, but we modify the coordinate transform ‘reference time’. We form dirty images with the tclean task in CASA with 5 arcmin pixel size. The top row of Fig. 5 shows the resulting images, which are visually indistinguishable. The central column of the figure corresponds to the simulation in which the simulation time matches the reference time.

However, the bottom row shows each image as a residual to the central column – i.e. the difference between assuming an offset reference time compared to a precise coordinate transformation. Here, it is apparent that the point source is slightly shifted as the reference time diverges from the simulation time. Similar to how this affects the delay spectrum in Fig. 2, the further the difference between the observation time and the reference time, the larger the error.

This is expected as source position errors will lead to phase errors and will thus result in shifting or smearing of the source in the image. This error will likely be more acute when simulating an array with longer baselines. We leave this demonstration as a precaution for the readers, and obviously, the best practice is to simulate one snapshot at a time and set the reference time to the observation time.

## 6 COMPUTATIONAL PERFORMANCE AND OPTIMIZATION

An important feature of *matvis* is its computational performance. Indeed, the consideration of performance is the driving motivation behind the matrix-based algorithm presented in Section 3.5, and the structuring of the Python implementation in a way that GPUs can be used as a drop-in replacement for CPUs. In this section, we will

look at the relative contribution to overall run-time from each of the algorithm steps, determining the bottleneck(s). We will also explore how this run-time scales with different simulation parameters, numbers of cores, and how much acceleration is afforded by the GPU. We will also consider the memory footprint of the algorithm and associated concerns, and finally we will take a deeper look at the largest bottleneck – the matrix–matrix multiplication of step (v) – considering its theoretical complexity and potential optimizations.

### 6.1 Runtime performance

Fig. 6 summarizes the runtime performance of *matvis*. Specifically, it shows the timing of a number of simulations with different numbers of source $N_{sources}$, antennas $N_{antennas}$, and times/snapshots $N_{times}$, using the script *matvis profile* included in the package. The script generates simulation configurations that take advantage of the fact that the performance cannot depend on the properties or locations of sources, the locations of antennas, or the precise frequencies or times of observation: only their total number. In addition to these dimensions, the performance will depend on the number of unique beams, $N_{beams}$, and whether they are analytic or simulated (and if simulated, the size of the grid onto which they are simulated), as well as the number of feeds per antenna, $N_{feed}$ and whether the beams are *E*-field beams or not.

In Fig. 6 we consider only the case where $N_{beams} = 1$, i.e. all antennas have the same beam (which is a common assumption). This beam is a complex-valued *E*-field beam evaluated on a 1-deg regular mesh for two feeds. To be representative of common simulation setups, the sources are distributed uniformly across the sky, which means that at any particular integration, half of them are below the horizon. In the current implementation of *matvis*, these sources are ignored (for that given integration), which potentially impacts performance.

Fig. 6 shows the impact of three parameters. Our ‘default’ value for each is $N_{sources} = 2^{19} \approx 525$ k, $N_{antennas} = 256$ and $N_{times} = 16$. Keeping two parameters fixed at their default, we run the third parameter over a logarithmic grid: $N_{sources} \in 2^{[10,...,19]}$, $N_{antennas} \in 2^{[3,...,8]}$, and $N_{times} \in 2^{[10,...,7]}$. Run-times are shown for the CPU (orange) and GPU (blue), with the default run shown in black. We also show the run-time for the default set using *pyuvsim* as a black
<figure><img src="image_7.png" /><figcaption>**Figure 5.** Demonstration of the effects of errors in matvis 'fast' coordinate transformation on the dirty images of its output visibilities. Here, a snapshot of a single 100-Jy point source transiting near the zenith of the full HERA array was simulated, and then converted to a dirty image with the tclean task on CASA. The top row shows the imaging results, for which each column is derived from a simulation with a different matvis coordinate transformation reference time relative to the simulated time. The bottom row shows the residual differences of each column to the middle column, for which the start and reference times are the same, thus providing the most accurate source position. As expected, the source position error due to matvis coordinate transformation scheme leads to a position shift of the point source in the image.</figcaption></figure>

square. Different marker symbols indicate steps of the algorithm (cf. Section 3.5), which we will discuss in more detail momentarily. These profiling simulations were performed on the Ilifu cloud-computing cluster<fnref n="19" /> in Cape Town, South Africa. The CPU profiles are computed on a 2.6 GHz Xeon Processor using a single core for each run. The GPU profiles are computed on NVIDIA V100 GPU. Note that the results shown are for a single frequency; there is almost perfect parallelism between frequencies, and virtually no overhead from computing a single frequency in a single job. There are two caveats to this: first, the computation of the array **R** for all time-steps can be done once for all frequencies, and secondly, writing out the simulated data to disc can be more efficient when frequencies are grouped together. The computation of **R** is included in the total time shown in Fig. 6, contributing to the small overhead. I/O overheads from reading the beam are included, but not from writing the results to disc. The total run-time (for both CPU and GPU) is calculated by simply measuring the wall-time of the entire calculation. The run-time of individual algorithm steps is measured on the CPU by using the line profiler Python package. On the GPU, since each step corresponds to a particular CUDA kernel, we use the nvprof utility which specifically measures performance of kernels at runtime. We note that the sum of the run-times of individual steps should be strictly less than the total time, which also takes into account overheads outside the main loop of the algorithm as well as internal overheads (e.g. Python for-loop overhead).

The first thing we notice is that the scaling with $N_{\text{times}}$ is almost perfectly linear for all steps of the calculation. This is exactly as we expect, since each step is run once per loop and the overhead outside the loop is very minimal compared to the run-time of the

default parameters. The scaling with $N_{\text{antennas}}$ and $N_{\text{sources}}$ is far more interesting, and varies amongst the different algorithm steps.

A second obvious takeaway is that the GPU accelerates the computation by a factor of ~100 for the larger simulations tested here. This is clearly dependent on the properties of both the CPU and GPU used, but should be a good rule-of-thumb. Even the slower CPU implementation is ~40× faster than the pyuvsim simulator in its current form.

For both GPU and CPU, the major component in terms of run-time is the matrix–matrix multiplication, step (v). For the majority of problem sizes on the GPU, this step is a factor of ~5 more expensive than the next-most-expensive step, which is generally the calculation of the per-antenna matrix **Z**, i.e. step (iv). Interestingly, on the CPU this is only true for the largest $N_{\text{antennas}}$ considered here: for smaller $N_{\text{antennas}}$, step (iv) is dominant. The scaling of these two steps with $N_{\text{antennas}}$ is expected to be different, with step (v) naively scaling more strongly than step (iv), which explains this cross-over for large $N_{\text{antennas}}$. Another caveat to this picture is that on the GPU, when $N_{\text{sources}} \lesssim 10^5$ or $N_{\text{antennas}} \lesssim 128$, the dominant contribution to run-time is in fact overheads, which for this setup appear to be about one second. While some of this overhead is from outside the time-loop, inspection of the profile reveals that much of the overhead is from within the time-loop (between GPU operations), and therefore will scale with $N_{\text{times}}$. This will present an issue for simulations with small $N_{\text{ant}}$ with large $N_{\text{times}}$, but can perhaps be overcome with further optimizations to the Python code. A final caveat is that while the beam interpolation is a negligible component here, this is partially because we have used $N_{\text{beam}} = 1$, whereas in principle it could scale with $N_{\text{antennas}}$ if all antennas are simulated with varying beams. If this is the case, it would become the dominant step on the CPU up to very large $N_{\text{antennas}}$ (where step (v) would overtake due to its hyperlinear scaling), and on the GPU it would be similar to step (iv) in its impact.

<footnote n="19">¹⁹https://www.ilifu.ac.za/</footnote>
<figure><img src="image_8.png" /><figcaption>**Figure 6.** Scaling test results of the CPU and GPU implementation of *matvis*. Orange and blue colours indicate CPU and GPU results, while the different markers indicate steps in the calculation as described in Section 3.5. The wall time for the default simulation parameters ($N_{\text{antennas}} = 256$, $N_{\text{sources}} = 2^{19}$, $N_{\text{times}} = 16$) are shown in a black circle markers for *matvis* and a black square marker for *pyuvsim* for comparison. Note that scaling tests of *pyuvsim* indicate that it scales similarly to *matvis* except at very small simulation sizes where it has some higher overheads.</figcaption></figure>

Since the algorithm’s performance is indifferent to whether baseline vectors are unique or not, the scalings shown in Fig. 6 hold both for highly redundant (e.g. HERA) and less redundant arrays (e.g. MWA, SKA). The performance of *matvis* for next-generation SKA-LOW with 512 antennas can thus be predicted to be approximately 3 s per channel per time-stamp for a sky model with $\sim 5 \times 10^5$ sources.

## 6.2 Multi-core CPU performance

Modern linear algebra and interpolation routines, such as those in NumPy and SciPy, can take advantage of multiple CPU cores through multi-threading. To measure this performance, we run the profiling script with the default parameter set giving 2, 4, 8, 16, and 32 number of CPUs per job. We compile NumPy with the Intel Math Kernel Library, and set the MKL_NUM_THREADS environment variable to the number of CPUs used in each job to specify the number of threads. Fig. 7 shows the scaling test results. Different markers correspond to different steps in *matvis* as in Fig. 6. The grey dashed line indicate linear scaling (i.e. -1 slop). Notice that only step (v), which forms the per-baseline visibilities through the matrix product of the per-antenna visibilities, benefits from multi-threading. The scaling is linear until ≈ 4 CPUs, and then it flatten off. While further CPU cores do continue to speed up step (v), step (iv) becomes the bottleneck by that point since it is simply performing complex array arithmetic, and it does not benefit from multi-threading. Thus, we recommend that each job is given 2–4 cores when simulating with the CPU implementation.

## 6.3 Optimization

Given the above observations, the crucial components to optimize for the majority of current use cases are steps (iv) and (v) – and particularly (v) as the number of antennas grow larger in modern experiments (for example, SKA Phase II will have 512 antennas, double the maximum tested in this section).

Step (iv) naively has complexity $O(N_{\text{ant}}N_{\text{src}})$, and this is borne out in Figs 6 and 7. It requires calculation of a complex-exponential for each of the terms, which can be quite computationally expensive. However, it is also naively parallelism over either the source or antenna axes. Which is more convenient depends on the memory layout of the input and output arrays.

Step (v) is a more interesting problem. Naively, its complexity is $O(N_{\text{ant}}^2 N_{\text{src}})$. This scaling applies inescapably to traditional simulators that compute the RIME as a set of $N_{\text{baseline}}$ tasks, of each which computes equation (3), since $N_{\text{baseline}} \propto N_{\text{ant}}^2$. However, our matrix-based approach makes use of highly optimized algorithms for one of the most common operations in all computing: the matrix–matrix product. The matrix $Z$ has shape $N \times M$, where $N = N_{\text{ant}}$ and $M = N_{\text{src}}$ ($N_{\text{src}} \sim 10^{2-3}$ and $M = N_{\text{ax}} N_{\text{src}} \sim 10^{5-6}$). For square matrix multiplication ($M = N$), the optimal complexity is known to be $O(N^\omega)$, with $\omega \approx 2.37$ or better (Vassilevska Williams et al. 2023). For more generic rectangular matrix multiplication, Le Gall & Urrutia (2018) places $\omega(k) < 2 + k$ as an upper bound for multiplying an $N \times N^k$ matrix by an $N^k \times N$ matrix for $k \le 5$ (e.g. $\omega \le 6.16$ for $k = 5$). Thus, an optimized matrix multiplication algorithm should perform better than $O(N_{\text{ant}}^2)$ even in the realistic situation where $N_{\text{src}} \gg N_{\text{ant}}$.
<figure><img src="image_9.png" /><figcaption>**Figure 7.** Scaling test results for multi-cores CPU performance of matvis. The grey dashed line indicates a linear scaling.</figcaption></figure>

We have checked several commonly used Python routines, and found the complexities in the range of $O(N_{ant}^{1.3-2})$, reaffirming this claim. Thus, the algorithm presented here has every chance of being unbeatable for simulating current and future large-$N_{ant}$ experiments.

In detail, some choices need to be made about the implementation of the matrix–matrix product. For both the CPU and GPU the current implementation uses the ZGEMM routine from BLAS to perform the multiplication. Specifically, on the CPU, the call is made via the `numpy.dot` function, which internally uses a backend implementation of BLAS that is configurable by the user. In Figs 6 and 7, we used the MKL backend. On the GPU, we use the CuBLAS implementation. Comparison to other implementations (e.g. `ArrayFire`<fnref n="20" /> and `Magma`<fnref n="21" />) is ongoing.

Besides the choice of implementation package, there are further optimizations worth considering. For example, the product $V = ZZ^†$ contains Hermitian symmetry that in principle can be exploited by the ZHERK routine in BLAS, reducing the computations by a factor of two. However, we have found this approach to be less performant than the full ZGEMM for a range of realistic problem sizes. This is most likely due to the highly irregular size of the matrices involved (i.e. $M \gg N$).

Another set of symmetries that this approach does not exploit are baseline redundancies: two baselines: $ij$ and $kl$ have equal visibilities if their beams are the same ($A_i \equiv A_k$ and $A_j \equiv A_l$) and their baseline vector is the same: $b_{ij} = b_{kl}$, regardless of the absolute positions of the antennas. For compact arrays like HERA – under the common assumption of ideal beams – the fraction of baselines which are redundant is very large ($\gtrsim 95$ per cent), which means that ZGEMM is naively performing a factor of $> 20 - 40$ more calculations than necessary. We consider this problem in a little more detail in Appendix A, noting here only that for very large $N_{ants}$, the optimal scaling of the matrix-multiplication will still win, even with the redundant calculations. In particular, matvis should be highly performant for upcoming large-$N$ non-redundant arrays such as the SKA.

## 6.4 Memory considerations

The way that the matvis algorithm is structured means that all steps inside the algorithm loop need only be stored in memory for a single snapshot, and only the final result needs to be kept for all snapshots. This final output is a complex array of size $N_{ant}^2 N_{feed}^2 N_{time}$. For the common case of $N_{feed} = 2$ and double-precision calculations, this is an array of size

$$ size(V) \approx \left( \frac{N_{ant}}{128} \right)^2 \frac{N_{time}}{1000} \text{ GB.} \qquad (18) $$

Inside the loop, the largest array is $Z$, with

$$ size(Z) \approx \frac{N_{ant}}{128} \frac{N_{src}}{130\,000} \text{ GB.} \qquad (19) $$

For the GPU implementation, this latter size is very important, as it is the dominant source of device-based memory. If $N_{src}$ is very large, the peak memory (and device memory) can be easily reduced by splitting the sources into chunks, for minimal overhead (since the matrix–matrix calculation scales linearly with $N_{src}$), and matvis does this automatically.

Likewise, the total peak memory can be controlled by chunking the number of times, which incurs minimal overhead (according to Fig. 6) as long as $N_{ant}$ and $N_{src}$ are reasonably large.

## 7 CONCLUSIONS

In this paper, we introduced the matvis visibility simulator and its matrix-based algorithm for calculating the RIME. By first calculating per-antenna responses and then expanding them to visibilities through a simple outer product, matvis is able to utilize highly efficient matrix operations, especially those implemented for GPUs, to accelerate the calculation of the RIME. The algorithms assume a point-source approximation for the sky model, treating each of the components in the model as an unresolved point source with a defined total flux density at a given frequency. To simulate diffuse sky emission or extended sources, they must first be discretized into an ensemble of point sources. To improve the performance of matvis, equatorial-to-horizontal coordinate transformations are performed as simple rotations, which introduces a slight degradation of accuracy if simulating a series of snapshots over a long interval. To mitigate this degradation, we introduce a computational trick that uses accurate astropy-based transformations to partially correct the errors in the fast geometric rotations, providing a reference time for the correction. This ‘fast’ coordinate transformation scheme yields the same results as astropy, and thus the most accurate source positions, when simulating from a single time snapshot.<fnref n="22" />

The primary use case of matvis, which is also the motivation for its development, is to generate mock visibilities from many-element 21 cm arrays, particularly HERA, and for validating the power spectrum analysis pipelines of these experiments. However, the presented algorithm and the matvis Python package itself are general and can be applied to other use cases that require mock visibilities. We have also developed the `hera_sim` package that provides HERA-specific systematic models and a unified interface for visibility simulators with already built-in wrappers for matvis and `pyuvsim`, a reference simulator that prioritizes

<footnote n="22">²²We note again that this approximate scheme is no longer required in versions of matvis greater than 1.3.0, and fast-and-accurate schemes can be used instead.</footnote>

<footnote n="20">²⁰https://arrayfire.org</footnote>

<footnote n="21">²¹https://icl.utk.edu/magma/</footnote>
accuracy and has been tested against analytical calculations. The `hera_sim` wrapper is our recommended way to get started using `matvis`.

We validated `matvis` for its 21 cm analysis use cases, particularly by looking at the dynamic range of the delay spectrum and power spectrum estimates from the simulations. We also compared `matvis` against `pyuvsim` to test its accuracy when simulating single time snapshots. For the delay spectrum, we simulated a single time snapshot from a mock full-sky point source catalogue derived from the source count distribution of the GLEAM (Franzen et al. 2019) catalogue, and a GSM diffuse sky model (Oliveira-Costa et al. 2008) with both simulators, and show that delay spectrum derived from these foreground simulations are well-matched to `pyuvsim` to ~$10^{-10}$ in fractional residuals. However, when simulating long-duration observations for the purposes of calibration, we recommend that the simulation be split into ~15 min time chunks to alleviate the errors from the approximate coordinate transforms. We also show that the obtained foreground delay spectrum maintains a peak-to-floor dynamic range $\gtrsim 10^{12}$ regardless of the reference time used for `matvis`' fast' coordinate transforms. We also show that simulated visibilities drawn from a 21 cm-only sky model, with known input power spectrum, produce power spectra that closely match the input power, within sample variance. These tests support `matvis` as a suitable simulator for the validation of 21 cm delay-spectrum pipelines. Furthermore, since residuals against `pyuvsim` are $\lesssim 10^{-10}$ when simulating a single time snapshot, `matvis` simulations performed in this way (one snapshot at a time) are of sufficient accuracy for the purposes of calibration of wide-field 21 cm data (given a suitable sky model).

We also demonstrate the fundamental limitation of `matvis`'s point-source approximation (and other simulators that make this assumption) when simulating from discretized diffuse emission at long baselines. In short, the resulting delay spectrum (or any measurement) derived from the simulated visibility will contain errors due to aliasing in the fringe term of the RIME. The rule of thumb is to make sure that the diffuse sky is discretized on a coordinate grid, e.g. HEALPIX, at least twice the angular resolution of the longest baseline in the simulation, i.e. following the Nyquist-Shannon sampling theorem.

There are currently a few limitations to `matvis`: (i) it can only simulate drift-scan observations; (ii) it can only simulate unpolarized (Stokes I) sky models, although it accounts for polarized beam effects; and (iii) it cannot simulate just a subset of the pairs of a given set of antennas, for example just the non-redundant baselines in HERA. One way to enable 'tracking observations' to address (i) is to perform a post-processing step that rephases the output visibilities to a chosen phase centre. As for (ii), the algorithm presented in equations (10) and (11), which is a general case of the current `matvis` algorithm, is already capable of handling a fully polarized sky. However, implementing it would require an additional tensor summation step, which must be tested and will affect the performance. Solving (iii) would require changing (or adding) how the visibility is computed in equation (3) to support a partial antenna matrix. Lastly, although not an essential limitation, the fast coordinate transformation approximation is a rough edge that can currently only be overcome by a slow outer-loop chunking simulations over snapshot times. This limitation is currently being removed so that arbitrarily high-precision long-duration simulations will be possible with future versions of `matvis`. We leave these as features for future versions of `matvis` and welcome contributions from the community.

## ACKNOWLEDGEMENTS

This material is based upon work supported by the National Science Foundation under grants #1636646 and #1836019 and institutional support from the HERA collaboration partners. This research is funded in part by the Gordon and Betty Moore Foundation through Grant GBMF5212 to the Massachusetts Institute of Technology. HERA is hosted by the South African Radio Astronomy Observatory, which is a facility of the National Research Foundation, an agency of the Department of Science and Innovation.

PK acknowledges the financial assistance of the South African Radio Astronomy Observatory (SARAO; www.sarao.ac.za). This result is part of a project that has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (Grant agreement No. 948764; PB). PB and HG acknowledge support from UK Science and Technology Facilities Council (STFC) Grants ST/T000341/1 and ST/X002624/1. SGM has received funding from the European Union's Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie grant agreement No. 101067043. CC is supported by the Beus Center for Cosmic Foundations. This work used Bridges-2 at the Pittsburgh Supercomputing Centre (PSC) through allocation PHY210042 from the Advanced Cyberinfrastructure Coordination Ecosystem: Services & Support (ACCESS) program, which is supported by National Science Foundation (NSF) grants #2138259, #2138286, #2138307, #2137603, and #2138296. We acknowledge the use of the ilifu cloud computing facility – www.ilifu.ac.za, a partnership between the University of Cape Town, the University of the Western Cape, Stellenbosch University, Sol Plaatje University, the Cape Peninsula University of Technology, and the South African Radio Astronomy Observatory. The ilifu facility is supported by contributions from the Inter-University Institute for Data Intensive Astronomy (IDIA – a partnership between the University of Cape Town, the University of Pretoria and the University of the Western Cape), the Computational Biology division at UCT and the Data Intensive Research Initiative of South Africa (DIRISA).

## DATA AVAILABILITY

`matvis` and `hera_sim` have been developed or contributed by members of the HERA collaboration. `pyuvsim` and `pyuvdata` are developed and maintained by the Radio Astronomy Software Group (RASG). All packages mentioned are publicly available on GitHub (see Section 4 for details). Visibility files from Section 5 can be made available upon request to the corresponding author.

## REFERENCES

* Aguirre J. E. et al., 2022, ApJ, 924, 85

* Astropy Collaboration, 2013, A&A, 558, A33

* Astropy Collaboration, 2018, AJ, 156, 123

* Astropy Collaboration, 2022, ApJ, 935, 167

* Barry N., Chokshi A., 2022, ApJ, 929, 64

* Barry N., Hazelton B., Sullivan I., Morales M. F., Pober J. C., 2016, MNRAS, 461, 3135

* Barry N., Beardsley A. P., Byrne R., Hazelton B., Morales M. F., Pober J. C., Sullivan I., 2019, Publ. Astron. Soc. Aust., 36, e026

* Bassett N., Rapetti D., Tauscher K., Nhan B. D., Bordenave D. D., Hibbard J. J., Burns J. O., 2021, ApJ, 923, 33

* Berkhout L. M. et al., 2024, PASP, 136, 045002

* Cheng C. et al., 2018, ApJ, 868, 26
Choudhuri S., Bull P., Garsden H., 2021, MNRAS, 506, 2066

Cohen A., Fialkov A., Barkana R., 2018, MNRAS, 478, 2193

DeBoer D. R. et al., 2017, PASP, 129, 045001

Fagnoni N., de Lera Acedo E., Drought N., DeBoer D. R., Riley D., Razavi-Ghods N., Carey S., Parsons A. R., 2021, IEEE Trans. Antennas Prop., 69, 8143

Franzen T. M. O., Vernstrom T., Jackson C. A., Hurley-Walker N., Ekers R. D., Heald G., Seymour N., White S. V., 2019, Publ. Astron. Soc. Aust., 36, e004

Gehlot B. K. et al., 2021, MNRAS, 506, 4578

Gorce A. et al., 2023, MNRAS, 520, 375

Górski K. M., Hivon E., Banday A. J., Wandelt B. D., Hansen F. K., Reinecke M., Bartelmann M., 2005, ApJ, 622, 759

HERA Collaboration, 2023, ApJ, 945, 124

Hamaker J. P., 2000, A&AS, 143, 515

Hamaker J. P., 2006, A&A, 456, 395

Hamaker J. P., Bregman J. D., Sault R. J., 1996, A&AS, 117, 137

Hazelton B. J., Jacobs D. C., Pober J. C., Beardsley A. P., 2017, J. Open Source Softw., 2, 140

Hurley-Walker N. et al., 2017, MNRAS, 464, 1146

Jacobs D. C. et al., 2017, PASP, 129, 035002

Josaitis A. T., Ewall-Wice A., Fagnoni N., de Lera Acedo E., 2022, MNRAS, 514, 1804

Kern N. S., Parsons A. R., Dillon J. S., Lanman A. E., Fagnoni N., de Lera Acedo E., 2019, ApJ, 884, 105

Kern N. S. et al., 2020, ApJ, 888, 70

Kohn S. A. et al., 2019, ApJ, 882, 58

Kolopanis M. et al., 2019, ApJ, 883, 133

Lanman A. E., Kern N., 2019, Astrophysics Source Code Library, record-ascl:1907.002

Lanman A. E., Pober J. C., 2019, MNRAS, 487, 5840

Lanman A., Hazelton B., Jacobs D., Kolopanis M., Pober J., Aguirre J., Thyagarajan N., 2019, J. Open Source Softw., 4, 1234

Lanman A. E., Murray S. G., Jacobs D. C., 2022, ApJS, 259, 1, 22

Le Gall F., Urrutia F., 2018, Proc. Annual ACM-SIAM Symposium on Discrete Algorithms (SODA). p. 1029

Line J. L. B. et al., 2020, Publ. Astron. Soc. Aust., 37, e027

Martinot Z. E., 2022, PhD thesis, University of Pennsylvania. Available at: https://repository.upenn.edu/dissertations/AAI29169424

Mort B. J., Dulwich F., Salvini S., Adami K. Z., Jones M. E., 2010, in IEEE International Symposium on Phased Array Systems and Technology. IEEE, Waltham, Massachusetts, p. 690

Noordam J. E., Smirnov O. M., 2010, A&A, 524, A61

de Oliveira-Costa A., Tegmark M., Gaensler B. M., Jonas J., Landecker T. L., Reich P., 2008, MNRAS, 388, 247

Pagano M. et al., 2023, MNRAS, 520, 5552

Parsons A. R., Pober J. C., Aguirre J. E., Carilli C. L., Jacobs D. C., Moore D. F., 2012, ApJ, 756, 165

Perkins S. J., Marais P. C., Zwart J. T. L., Natarajan I., Tasse C., Smirnov O., 2015, Astron. Comput., 12, 73

Santos M. G., Cooray A., Knox L., 2005, ApJ, 625, 575

Sault R. J., Teuben P. J., Wright M. C. H., 1995, in Shaw R. A., Payne H. E., Hayes J. J. E., eds, ASP Conf. Ser. Vol. 77, Astronomical Data Analysis Software and Systems IV. Astron. Soc. Pac., San Francisco, p. 433

Shaw J. R., Sigurdson K., Pen U.-L., Stebbins A., Sitwell M., 2014, ApJ, 781, 57

Shaw J. R., Sigurdson K., Sitwell M., Stebbins A., Pen U.-L., 2015, Phys. Rev. D, 91, 083514

Shaw J. R., Foreman S., Nitsche R., Masui K., Zuo S., Berger P., 2023, Astrophysics Source Code Library, record ascl:2308.013

Smirnov O. M., 2011a, A&A, 527, A106

Smirnov O. M., 2011b, A&A, 527, A107

Smirnov O. M., 2011c, A&A, 527, A108

Smirnov O. M., 2011d, A&A, 531, A159

Spreeuw H., Yatawatta S., van Werkhoven B., Diblen F., 2020a, preprint (arXiv:2010.15661)

Spreeuw H., Yatawatta S., Van Werkhoven B., Diblen F., 2020b, preprint (arXiv:2010.16209)

Spreeuw H., van Werkhoven B., Yattawatta S., Diblen F., 2020c, in Pizzo R., Deul E. R., Mol J. D., de Plaa J., Verkouter H., eds, ASP Conf. Ser. Vol. 527, Astronomical Data Analysis Software and Systems XXIX. Astron. Soc. Pac., San Francisco, p. 747

Sullivan I. S. et al., 2012, ApJ, 759, 17

The CASA Team, 2022, PASP, 134, 114501

Vassilevska Williams V., Xu Y., Xu Z., Zhou R., 2023, preprint (arXiv:2307.07970)

Wilensky M. J. et al., 2023, ApJ, 957, 78

## APPENDIX A: UTILIZING BASELINE REDUNDANCIES

Let the matrix $Z$ have shape $N \times M$, where $N = N_{\text{feed}}N_{\text{ant}}$ and $M = N_{\text{ax}}N_{\text{src}}$. Then the matrix $V = ZZ^\dagger$ has shape $N \times N$. As discussed in Section 6.3, one of the disadvantages of the current matvis implementation is that all $N^2$ antenna-feed pairs must be simulated, including those that are 'redundant'.

To be clear, this redundancy occurs when two baselines, $ij$ and $kl$, have equivalent beams and baseline vectors ($b_{ij} = b_{kl}$). We can split the baselines of any array into a number $N'_u$ of sets of baselines, where each set contains all $n_i$ baselines that are redundant with each other, such that $\sum_i n_i = N_{\text{ant}}(N_{\text{ant}} + 1)/2$ (the sets do not contain reverse-pairs since their visibilities are simply the conjugate of the reverse). Letting $N_u = N_{\text{feed}}^2 N'_u$ be the number of non-redundant elements of $V$, we can write the 'non-redundancy fraction',

$$ \xi = N_u / N^2, \qquad (A1) $$

which is very low for highly redundant compact arrays like HERA. Thus, calculating the full $N^2$ matrix is performing far more calculations than strictly necessary.

In this appendix, we consider a few alternative approaches aimed at improving performance in this situation. Far from being a fully fledged examination, we will limit our focus to a few ideas, and test their performance only on 'HERA-like' arrays, by which we mean that the antennas are arranged in a hexagonal core, split into three parallelogram shards (reducing redundancy) with $N_{\text{hex}}$ antennas along a side. In this setup, $\xi$ varies between ~1 and 3 per cent, and $N'_u$ varies between 500 and 6000 (see Fig. A1).

<figure><img src="image_10.png" /></figure>

**Figure A1.** Characteristics of the 'HERA-like' arrays used in this appendix. The x-axis denotes the size of the array, defined by $N_{\text{hex}}$ (the number of antennas along a side of the hexagon) and $N_{\text{ants}}$. The orange and green curves (left spines) show characteristics of the arrays: the number of unique pairs (orange) and the non-redundancy fraction, $\xi$. The blue and red curves (right spines) show indicate properties of the sub-matrix algorithm outlined in this appendix: the total size of the sub-matrices (in this case, three sub-matrices were used) in blue, and this size as a fraction of the full matrix. Unique pairs grow almost linearly with $N_{\text{ants}}$.
<figure><img src="image_11.png" /><figcaption>**Figure A2.** Performance characteristics of four different approaches to calculating **V**. The LHS panel shows performance as a function of $N_{ant}$, where the array is a 'HERA-like' array with characteristics shown in Fig. A1. Here $N_{side} = 256$ is used. On the RHS, $N_{ant}$ is fixed at 320 and the performance as a function of $N_{src}$ is shown. The full matrix–matrix method adopted by *matvis* is in blue, the matrix–vector approach in orange, the vector–vector dot product approach in green, and the more general sub-matrix approach in red and purple (for different numbers of sub-matrices). Black dashed line indicates linear scaling, and black dotted indicates quadratic scaling.</figcaption></figure>

Fig. A2 shows in blue the performance of the default implementation of *matvis*: the full matrix–matrix (MM) calculation. On the LHS the performance as a function of $N_{ant}$ is shown, and on the RHS as a function of $N_{src}$. The black-dashed line is an indication of linear scaling with each parameter, whereas the black-dotted line indicates scaling with $N_{ant}^2$. It is clear that the GPU implementation used here scales close to $N_{ant}^2$.

The most straightforward way to use the redundancy of the array is to compute only the $N_u$ vector–vector (VV) dot products required within a for-loop. This has complexity $O(N_u M \equiv \xi N^2 M)$, which is better than the $O(N^2 M)$ that the full matrix–matrix calculation exhibits in Fig. A2. On the other hand, we expect this to have some overhead associated with the for-loop, and this is borne out in Fig. A2, where this method is shown in green. For small $N_{ant}$, the overhead associated with the for-loop is large enough that this method is ~6x slower than the MM method. However, it grows with $N_u$ instead of $N^2$, which in the HERA-like case considered here is close to $N_u \propto N$. For the full HERA array, the VV loop is more performant than the MM calculation. Considering the performance against $N_{src}$, the VV loop has very poor performance for small $N_{src}$, and only starts to hit the theoretical linear growth for $N_{src} \gg 10^6$. This must be due to low utilization of the many GPU cores: each VV dot product is performing $N_{src}$ calculations split between thousands of cores.

A more general alternative approach is to deconstruct **V** into a set of smaller matrices {$V'_x$}, for which each $V'_x$ contains some non-zero number of unique pairs in **V**, and is a smaller yet large enough matrix to take advantage of modern matrix multiplication routines,

$$V'_x = A_x B'_x \qquad (A2)$$

Here, $A_x$ and $B_x$ are both formed from (possibly non-adjacent) full rows of **Z**, having shape $N_A \times M$ and $N_B \times M$ respectively. The rows chosen for $B_x$ correspond to the columns of $V'_x$. We can then attempt to optimise the sets {$A_x$} and {$B_x$} to yield the best possible computation time for {$V'_x$} (and thus for **V**) given a particular set

of unique antenna pairs. This optimization is obviously dependent on the antenna distribution. Note that both the vector dot-product approach described above and the full matrix product $V = ZZ^\dagger$ are special cases of equation A2, where the number of matrices in {$V'_x$} is $N_u$ and 1, respectively.

For the full HERA array, the number of possible choices for {$V'_x$} is extremely large (since $N = 700)^{23}$, which makes equation A2 a daunting optimization problem. However, the structure of the matrix can potentially lead to some useful ideas.

One such intermediate idea is for each $A_x$ to be a $1 \times M$ matrix, and each $B_x$ to include only the antennas for which the single antenna in $A_x$ appears in non-redundant pairs. This computes only the redundant pairs, like the VV method above. However, it is liable to utilize vectorization better than the VV approach as it includes fewer but larger matrix-vector (MV) products. On the other hand, memory management may be a problematic overhead in this approach since for each of the $N$ MV products, the matrix needs to be 'formed' by coalescing memory from non-adjacent rows in the full **Z** matrix. This is borne out in Fig. A2, where this method is shown in orange. While the scaling of the MV method with both $N_{ant}$ and $N_{src}$ is what we expect (and is therefore better than MM for some large $N$ where $\xi$ is very small), it is much slower than all the other methods, due to both the overhead in the for-loop ($N_{ant}$ iterations are required) and in the movement of data on the GPU to form $B_x$.

A more generic idea to optimise the choice of {$A_x$} and {$B_x$} is to place down sub-matrices that cover all unique entries of **V**, but are as small as possible in size. A heuristic algorithm for this is the following:

(i) Form a list of pairs of antennas, $U = [(a_1, b_1), (a_2, b_2), \dots, (a_{N_u}, b_{N_u})]$ where each pair is a representative from one of the

<footnote>23There are $N_u$ choices of number of matrices in the set, and for each such choice there a multitude of ways to configure the shape of each sub-matrix.</footnote>
<figure><img src="image_12.png" /><figcaption>**Figure A3.** A visualization of our sub-matrix procedure for taking advantage of antenna-pair redundancies in HERA. The plot is meant to be read from top-left to bottom-right. The top-left represents **V** for HERA-320, where black pixels indicate unique baselines. These black pixels are less than 3 per cent of the pixels. In the top-right, the rows have been re-ordered such that the top row has the most black dots, and the bottom row the least. On the bottom panel, we show the result of our algorithm for partitioning the rows into different numbers of subsets (1, 2, and 3 from left to right). In each case, after deciding on a vertical extent for the matrix (coloured region), the columns *within* the matrix are shuffled such that columns occupied by any row within the matrix are moved to the right. This is simply to visually indicate the final matrix density. All black pixels must be covered by a sub-matrix. The coloured-shaded regions show the matrices {$V'_x$} that can be calculated faster than the full matrix product. See the text for details.</figcaption></figure>

‘redundant groups’. Since the choice of representative is arbitrary, there are many ways to do this, but one simple way is to use the first pair of the group in lexicographic order.

(ii) Let $\alpha = \{a_1, a_2, \dots\}$ be the set of all ‘first’ antennas in $\mathcal{U}$, and $\beta$ the corresponding set of ‘second’ antennas. Let $N_\alpha$ be the set of multiplicities of $a_x$ in $\mathcal{U}$.

(iii) Obtain $\alpha'$, which is a list of the elements of $\alpha$, sorted in descending order of their multiplicity. That is, sort antennas in decreasing order of how many times they appear as first-antenna in a unique pair.

(iv) Decide on a number $N_{\text{sub}}$ of sub-matrices into which **V** should be partitioned: our partition will be limited to an ordered partition of $\alpha'$, i.e. the placement of $N_{\text{sub}} - 1$ cuts in $\alpha'$ such that we have a new set of sub-lists, $\{\alpha'_x\}$, for $0 \le x < N_{\text{sub}}$. Let the size of each subset be $n_x$, so that $\sum_x n_x = |\alpha'|$; then $\alpha'_0$ corresponds to the first $n_0$ values in $\alpha'$, $\alpha'_1$ the next $n_1$ values, etc. This defines the set $\{\mathbf{A}_x\}$ such that each element is the matrix formed by taking the rows of **Z** corresponding to the indices in $\alpha'_x$. Moreover, it uniquely defines the corresponding set $\{\beta_x\}$ which must be the set of all $b_i$ that exist in a pair in $\mathcal{U}$ where the first element is in $\alpha'_x$. The set $\{\mathbf{B}_x\}$ is the

analogue of $\{\mathbf{A}_x\}$ for $\{\beta_x\}$. Letting $V_x$ be the outer-product of pairs from $\alpha'_x$ and $\beta'_x$, $V_x = \{(i, j)\}\forall i \in \alpha'_x, \forall j \in \beta_x$, we denote the ‘size’ of a particular sub-matrix as $S_x = |V_x|$, and its ‘density’ $D_x$ is the ratio of the number of unique pairs in $\mathcal{U}$ that appear in $V_x$ to $S_x$, i.e. $D_x = |\mathcal{U} \cap V_x|/|V_x|$.

(v) To decide where the $N_{\text{sub}} - 1$ partitions should be placed, use an iterative procedure. First, notice that increasing the size of a sub-matrix $\mathbf{A}_x$ by including one extra index from $\alpha'$ is almost always going to reduce the density $D_x$, since $\alpha'$ is ordered by multiplicity. Using this, we simply do the following: for each partition point, defined by $n_x$, try all values from $n_x = 1, 2, \dots$, each time computing $D_x$. Once $D_x$ reduces below a threshold $\tau$, move to the next cut, $n_{x+1}$ and do the same. After placing all cuts, compute $D_{N_{\text{sub}}}$ for the last sub-matrix. If it is greater than $\tau$, then we’re done. If not, reduce $\tau$ by some pre-defined small amount, and redo the loop over $x$ with the new threshold.

The output of this algorithm is the ordered list $\alpha'$ and the set of cuts, $\{n_1, n_2, \dots, n_{N_{\text{sub}}}\}$, such that the densities $D_x$ are as high as possible while all remaining similar. Since computation time will
scale roughly as the total size, $S = \sum_x S_x$, we note that $S = N_u/D$, where by construction $D < \min\{D_x\}$. This process is illustrated in Fig. A3. This algorithm is by no means guaranteed to return the best possible decomposition of $V$: it is constrained to partitions only in $\alpha'$, not the joint space of $\alpha'$ and $\beta$ (in the visual representation of Fig. A3, different sub-matrices never appear on the same row), and it does not attempt to use other pairs from the redundant groups instead of the arbitrarily chosen representative pair. Nevertheless, it obviously does reasonably well in the case of HERA, with a compression factor of more than 15 for 3 sub-matrices compared to the full calculation. A trivial adaptation to the method is to reverse the pairs (i.e. pass $(b, a)$ into the algorithm instead of $(a, b)$). Our tests show a compression factor of ~25 in this case.

Fig. A2 shows the results of this sub-matrix method in red and purple (triangles), for 2- and 3-sub-matrices, respectively. We first note that these two decompositions yield very similar performance. Overall they scale as we would expect: linearly with $N_{\text{ant}}$ like the other methods that utilize redundancy, and linearly with $N_{\text{src}}$. Since each of the sub-matrices is reasonably large, we do not see issues of underutilization of the GPU parallelization for small number of sources, like we did in the VV method. The reason that both 2- and 3-sub-matrix cases perform similarly well seems to stem from the fact that they have similar total sizes (11.4 per cent and 6.6 per cent, respectively, c.f. Fig. A3), and that the modest decrease in raw computations afforded by three sub-matrices is counterbalanced by the overhead of an extra independent matrix multiplication. It is unlikely, at least in the case of HERA, that further sub-matrices would be effective, for this reason.

The strongest downside to this method is that, similar to the MV method, it needs to copy memory on the GPU to form the $A_x$ and $B_x$. Our tests show that this memory copying constitutes about 1/3 of the total runtime in our non-optimized implementation. Better memory management could make this method even more competitive. This is seen also in the $N_{\text{side}} = 512$ case, where there is a sudden jump in runtime. This is due to the extra memory consumption of this method; to get around this, a loop over chunks of $N_{\text{src}}$ is required, with the extra memory copies associated. The sub-matrix method may also not work when baselines are not perfectly redundant since matrix compression utilize redundancies. We plan to improve on the method and conduct thorough tests before adding it to matvis in the near future.

## APPENDIX B: THE SQUARE ROOT OF THE COHERENCY MATRIX

The form of the matrix $M$, where

$$MM^{\dagger} = C = \begin{pmatrix} I+Q & U+iV \\ U-iV & I-Q \end{pmatrix}, \qquad (B1)$$

is

$$M = \frac{1}{\sqrt{2}} \begin{pmatrix} \frac{(Q-T)\sqrt{\frac{I^2+I-T-Q}{(U-iV)}}}{\sqrt{(Q/T+1)(I-T)}} & \frac{(T+Q)\sqrt{-\frac{I^2+I+T-Q}{(U-iV)}}}{\sqrt{(1-\frac{Q}{T})(I+T)}} \\ \frac{\sqrt{\frac{I^2+I-T-Q}{(U-iV)}}}{\sqrt{(Q/T+1)(I-T)}} & \frac{\sqrt{-\frac{I^2+I+T-Q}{(U-iV)}}}{\sqrt{(1-\frac{Q}{T})(I+T)}} \end{pmatrix}, \qquad (B2)$$

where $T^2 = Q^2 + U^2 + V^2$.

<footnote n="1">¹Department of Physics and Astronomy, University of Western Cape, Cape Town 7535, South Africa</footnote>

<footnote n="2">²South African Radio Astronomy Observatory, Black River Park, 2 Fir Street, Observatory, Cape Town 7925, South Africa</footnote>

<footnote n="3">³Scuola Normale Superiore, I-56126 Pisa, PI, Italy</footnote>

<footnote n="4">⁴School of Earth and Space Exploration, Arizona State University, Tempe, AZ 85287, USA</footnote>

<footnote n="5">⁵Queen Mary University London, Department of Physics and Astronomy, London E1 4NS, UK</footnote>

<footnote n="6">⁶Jodrell Bank Centre for Astrophysics, University of Manchester, Manchester M13 9PL, UK</footnote>

<footnote n="7">⁷Department of Astronomy, University of California, Berkeley, CA 94720, USA</footnote>

<footnote n="8">⁸Department of Physics and Astronomy, University of Pennsylvania, Philadelphia, PA 19104, USA</footnote>

<footnote n="9">⁹Department of Physics, Winona State University, Winona, MN 55987, USA</footnote>

<footnote n="10">¹⁰Department of Physics and Trottier Space Institute, McGill University, 3600 University Street, Montreal, QC H3A 2T8, Canada</footnote>

<footnote n="11">¹¹INAF – Istituto di Radioastronomia, via Gobetti 101, I-40129 Bologna, Italy</footnote>

<footnote n="12">¹²Department of Physics and Electronics, Rhodes University, PO Box 94, Grahamstown 6140, South Africa</footnote>

<footnote n="13">¹³National Radio Astronomy Observatory, Charlottesville, VA 22903, USA</footnote>

<footnote n="14">¹⁴Cavendish Astrophysics, University of Cambridge, Cambridge CB3 0HE, UK</footnote>

<footnote n="15">¹⁵National Radio Astronomy Observatory, Socorro, NM 87801, USA</footnote>

<footnote n="16">¹⁶MIT Kavli Institute, Massachusetts Institute of Technology, Cambridge, MA 02139, USA</footnote>

<footnote n="17">¹⁷Department of Physics, Massachusetts Institute of Technology, Cambridge, MA 02139, USA</footnote>

<footnote n="18">¹⁸Centre for Strings, Gravitation and Cosmology, Department of Physics, Indian Institute of Technology Madras, Chennai 600036, India</footnote>

<footnote n="19">¹⁹Radio Astronomy Lab, University of California, Berkeley, CA 94720, USA</footnote>

<footnote n="20">²⁰Department of Physics, University of California, Berkeley, CA 94720, USA</footnote>

<footnote n="21">²¹Department of Physics and Astronomy, University of California, Los Angeles, CA 90095, USA</footnote>

<footnote n="22">²²Institut d'Astrophysique Spatiale, CNRS, Université Paris-Saclay, F-91405 Orsay, France</footnote>

<footnote n="23">²³School of Physics, University of Melbourne, Parkville, VIC 3010, Australia</footnote>

<footnote n="24">²⁴Department of Physics, University of Washington, Seattle, WA 98195, USA</footnote>

<footnote n="25">²⁵eScience Institute, University of Washington, Seattle, WA 98195, USA</footnote>

<footnote n="26">²⁶Department of Physics, Brown University, Providence, RI 02912, USA</footnote>

<footnote n="27">²⁷Department of Computer Science, University of Nevada, Las Vegas, Las Vegas, NV 89154, USA</footnote>

<footnote n="28">²⁸Department of Physics and McGill Space Institute, McGill University, 3600 University Street, Montreal, QC H3A 2T8, Canada</footnote>

<footnote n="29">²⁹Department of Physics, Stellenbosch University, Matieland, Western Cape 7602, South Africa</footnote>

<footnote n="30">³⁰School of Chemistry and Physics, University of KwaZulu-Natal, Westville Campus Private Bag X54001, Durban 4000, South Africa</footnote>

<footnote n="31">³¹International Centre for Radio Astronomy Research, Curtin University, Bentley, WA 6102, Australia</footnote>

<footnote n="32">³²ARC Centre of Excellence for All Sky Astrophysics in 3 Dimensions (ASTRO 3D), Bentley, WA 6102, Australia</footnote>

<footnote n="33">³³School of Physics, University of Melbourne, Parkville, VIC 3010, Australia</footnote>

<footnote n="34">³⁴Raman Research Institute, Sadashivanagar, Bengaluru 560 080, India</footnote>

<footnote n="35">³⁵Commonwealth Scientific and Industrial Research Organisation (CSIRO), Space and Astronomy, P.O. Box 1130, Bentley, WA 6102, Australia</footnote>

<footnote>This paper has been typeset from a TeX/LA TeX file prepared by the author.</footnote>