# Lossy Compression of Large-Scale Radio Interferometric Data

M Atemkeng,1 S Perkins,2 E Seck,3 S Makhathini,2,4,5 O Smirnov,2,4 L
Bester,2,4 B Hugo,2,4
1Department of Mathematics, Rhodes University, Grahamstown 6139, South Africa
2South African Radio Astronomy Oberservatory, Black River Park, 2 Fir Street,
Observatory, Cape Town, 7925, South Africa
3African Institute for Mathematical Sciences (AIMS), Limbe, Cameroon
4 Department of Physics and Electronics, Rhodes University, PO Box 94,
Grahamstown, 6140, South Africa
5School of Physics, University of the Witwatersrand, 1 Jan Smuts Avenue
Johannesburg, 2001, South Africa E-mail<EMAIL_ADDRESS>(MA)

(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

Radio telescopes produce vast amounts of data and the data volume is set to
increase as upcoming telescopes (e.g. SKA, ngVLA) come online. The vast
amounts of data are an important issue to deal with in the context of
calibration, deep wide-field imaging, storage and archiving for sub-sequence
processing. This work proposes to reduce visibility data volume using a
baseline-dependent lossy compression technique that preserves smearing at the
edges of the field-of-view. We exploit the relation of the rank of a matrix
and the fact that a low-rank approximation can describe the raw visibility
data as a sum of basic components where each basic component corresponds to a
specific Fourier component of the sky distribution. As such, the entire
visibility data is represented as a collection of data matrices from
baselines, instead of a single tensor. This has the benefit of parallel
computation and allows the investigation of baseline-dependent rank$-r$
decomposition. The proposed methods are formulated as follows: provided a
large dataset of the entire visibility data; the first algorithm, named simple
SVD projects the data into a regular sampling space of rank$-r$ data matrices.
In this space, the data for all the baselines has the same rank, which makes
the compression factor equal across all baselines. The second algorithm, named
BDSVD projects the data into an irregular sampling space of rank$-r_{pq}$ data
matrices. The subscript $pq$ indicates that the rank of the data matrix varies
across baselines $pq$, which makes the compression factor baseline-dependent.
MeerKAT and the European Very Long Baseline Interferometry Network are used as
reference telescopes to evaluate and compare the performance of the proposed
methods against traditional methods, such as traditional averaging and
baseline-dependent averaging (BDA). For the same spatial resolution threshold,
both simple SVD and BDSVD show effective compression by two-orders of
magnitude higher than traditional averaging and BDA. At the same space-saving
rate, there is no decrease in spatial resolution and there is a reduction in
the noise variance in the data which improves the S/N to over $1.5$ dB at the
edges of the field-of-view. The proposed compression methods offer superior
compression but requires processing data at full data resolution, for reduced
implementation complexity.

###### keywords:

Instrumentation: interferometers, Methods: data analysis, Methods: numerical,
Techniques: interferometric

††pubyear: 2015††pagerange: Lossy Compression of Large-Scale Radio
Interferometric Data–A

##  1 Introduction

Radio interferometric arrays consist of an assembly of radio antennas that are
correlated in pairs to produce complex data values, known as visibilities or
sampled spatial frequencies. The data volume grows quadratically with the
number of antennas, very large sky surveys, high spectral and temporal
resolutions. Processing and storing this volume of visibility data has become
a challenge, for example calibrating and translating the data to the image
space via resampling and fast Fourier transform operations. The large volume
of visibility data is an important problem to deal with in the context of
calibration and deep wide-field imaging with current radio interferometric
arrays such as MeerKAT (Jonas, 2009), ASKAP (Johnston et al., 2008), NenuFAR
(Zarka et al., 2015), LOFAR (van Haarlem et al., 2013) and future radio
interferometers, including the Square Kilometre Array (SKA)(Dewdney et al.,
2009). To resolve the morphological structure of compact sources, the SKA will
sample very high spatial frequencies with long baselines. This high spatial
frequency sampling will take the SKA to an unprecedented visibility data
volume era; this will require more computation and improved data compression
strategies. Long-term data archiving of calibrated products is necessary, and
will benefit from the implementation of new data compression algorithms.

Data compression is critical to reduce the costs involved in storing,
processing and archiving data. Many data compression strategies exist,
however, choosing a relevant strategy depends on the science case and how the
information related to the science case is encoded in the data. Therefore, a
high data compression rate requires that we understand the information in the
data and the science case. For a given dataset it is possible that only a few
data points contain the majority of information relevant to the science, while
other points contain noise or information irrelevant to the case (Bobin et
al., 2008).

In visibility data, for example, the noise from system electronics, signal
from unwanted areas of the sky and spatial frequencies from redundant
baselines can be removed from the data. There are mathematical tools that can
break down this type of dataset into a new form that makes it easier to
understand how the information in the dataset is represented based on some
degree of importance. The properties of each of these mathematical tools for
breaking down the dataset are different and their suitability depends on the
science case. Additional factors must also be considered when choosing a
compression strategy, to find the best tradeoff between the compression ratio,
signal loss and science.

In practice, the visibility data are integrated and averaged over finite time
and frequency intervals which introduces decorrelation effects. To limit
decorrelation, the finite time and frequency intervals are kept very small so
that the phase term on long baselines is preserved. The decorrelation
attenuates and changes the morphology of sources at the edges of the field-of-
view (FoV). If the time and frequency intervals are scaled up to some limit,
for a given FoV, then traditional averaging can be used to reduce the data
volume. Deconvolution routines must correctly take into account the smearing
effects to avoid limiting the dynamic range of the image; this is discussed in
details by several authors (Cotton, 2009; Atemkeng et al., 2016; Tasse et al.,
2018; Atemkeng et al., 2020).

Since high-resolution sampling of the spatial frequency is only required on
long baselines, baseline-dependent averaging (BDA) methods have recently
gained popularity and can be used to compress the visibility data and suppress
sources out of the FoV while maintaining spectral and temporal resolutions on
the long baselines. While BDA can potentially offer compression capabilities
on short baselines and maintain the high spectral and temporal resolution
required on long baselines, it still requires further investigation. The
Measurement Set v2.0 (MS)(Kemball & Wieringa, 2000) format does not support
storing variable frequency bins as imposed by BDA. The Measurement Set v2.0
format can, in principle, store variable frequency bins required by BDA, thus
this approach requires the creation of a spectral window for each
decomposition of the frequency domain. Current calibration algorithms expect
the data to be regularly sampled. However, to make use of BDA compression
capabilities, calibration algorithms will need to adapt, especially for the
different spectral and temporal resolutions along baselines. It is also
important to choose appropriate solution intervals for calibration when
choosing BDA parameters. Well-defined interpolation windows and weighting
schemes are needed for imaging. This is not a problem for imaging algorithms
as they place unstructured samples on a regular grid.

Visibility data is frequently archived for later use; traditional averaging
and/or BDA is applied before archiving. However, many criteria must be
considered when choosing compression parameters, as high spectral and temporal
resolution is imperative for certain science cases (e.g. VLBI science,
spectral line studies) and for calibration routines that take into account the
variation of the primary beam (PB) and the effect of the ionosphere, for
example. Traditional averaging and/or BDA is accompanied by a decrease in the
spectral and temporal resolution at least on shorter baselines for BDA and
must therefore be applied with care. There is, therefore, a need to
investigate data compression techniques that can archive visibility data with
little to no decrease in spectral and temporal resolution.

In this work, we use low-rank approximation for visibility data compression. A
low-rank approximation minimises a loss function that measures the best-
approximating data matrix with a reduced-rank relative to the original data
matrix i.e. a low-rank approximation represents a higher-rank data matrix by a
reduced-rank data matrix with little to no loss of information. A data matrix,
$\boldsymbol{V}$ of size $M\times N$ has a rank$-r$ approximation given by:

$\boldsymbol{V}\simeq\boldsymbol{A}\boldsymbol{C},$ (1)

where the data matrices $\boldsymbol{A}$ and $\boldsymbol{C}$ are of size
$M\times r$ and $r\times N$, respectively. The number of entries in the
rank$-r$ approximation is $r(M+N)$, which can now be stored in the place of
$\boldsymbol{V}$ to save memory and/or computation since $r(M+N)$ is
considerably smaller than $MN$ if $r$ is relatively small and $M,N$ are
relatively large. Thus, $1\leq r<MN/(M+N)$. $\mathcal{O}(r(M+N))$ disk space
is now required to save the rank$-r$ approximation of $\boldsymbol{V}$ rather
than $\mathcal{O}(MN)$.

As of now, in the big data regime, many fields applying image processing
project a high-rank dataset to a lower-rank by dimensionality reduction. That
is a process of feature (dimension) extraction that preserves the most
relevant information (reduction) in the data. Dimensionality reduction
primarily functions as a pre-processing step to reduce the need for high
computing resources, and is classified into two different classes; linear and
nonlinear. The most well-known low-rank approximation that belongs to the
linear class of dimensionality reduction algorithms is the singular value
decomposition (SVD). The SVD decomposes the data into separate sets of
relevant features, noisy and redundant components. The properties of the SVD
are fully explored in other fields. It is popular in signal processing due to
its orthogonal matrix decomposition. Several authors use SVD to tackle
different tasks with interferometric data; for example Offringa et al. (2010)
discussed using SVD to separate strong radio frequency interference from weak
astronomical signals. To speed up imaging and deconvolution, Kartik et al.
(2017) used SVD to compress interpolated and gridded visibility data. Recent
work proposes using holographic measurement of the PB and SVD to deal with
spatial and spectral variation of the PB; e.g. Iheanetu et al. (2019)
discussed the possibility of determining a sparse representation of the PB
with a few features and Ansah-Narh et al. (2020) discussed mitigating noise in
the PB by choosing the strongest components of the decomposition.

The interest in exploiting the SVD to compress the visibility data is not
accidental; high-sensitivity radio interferometric arrays are dominated by
short baselines: SVD exploits redundancies in short baselines and can isolate
noise in the data. A rank$-r$ SVD decomposition of a data matrix of size
$M\times N$ requires $\mathcal{O}(rMN)$ operations: $\mathcal{O}(N^{3})$ in
the case of a square data matrix (Gu & Eisenstat, 1996; Demmel et al., 2007;
Shishkin et al., 2019). Computation of this magnitude is impractical in a big
data regime. However, the SVD methods described in this paper should not be
confused with a method that treats the entire visibility data as a single
tensor and where there is an assumption of similarity between different
baselines. We exploit the relation of the SVD to the rank of a matrix and the
fact that the SVD describes the raw visibility data as a sum of components,
where each component corresponds to a specific Fourier component of the sky
distribution. As such, the entire visibility data is represented as a
collection of matrices derived for each baseline, instead of a single tensor.
This facilitates baseline-dependent rank$-r$ decomposition that can be
performed independently for each baseline, allowing for parallel computing.
The proposed methods are formulated below. Given a dataset of visibility data:

  * •

The first algorithm projects the data into a regular sampling space of
rank$-r$ data matrices. In this space, the data for all the baselines has the
same rank, which makes the compression factor equal across all baselines. We
referred to this algorithm as simple SVD. The simple SVD is shown to be
effective in compressing the visibility data by two-orders of magnitude higher
than traditional averaging with a reduction in the noise variance in the data.

  * •

The second algorithm projects the data into an irregular sampling space of
rank$-r_{pq}$ data matrices. The extra subscript $pq$ indicates that the rank
of the data matrix varies across baselines, which makes the compression factor
baseline-dependent. We referred to this method as baseline-dependent SVD
(BDSVD).

It should be noted that simple SVD and BDSVD are used differently compared to
traditional averaging and BDA. The compressed visibility data from traditional
averaging and BDA are used directly as input to other tasks such as imaging or
calibration, with a lower computation cost as opposed to the uncompressed
data. This leads to a decrease in spectral and temporal resolution. The simple
SVD and BDSVD are used to archive the visibility data and the data must be
decompressed before subsequent use. This has the advantage of maintaining the
spectral and temporal resolution of the original data. The use cases are
summarised in Table 1. In this paper, we implement the simple SVD and BDSVD
serially and provide detailed justification for parallel implementation.
Whereafter we suggest a parallel algorithm that can reduce the computation
from $\mathcal{O}\big{(}N_{b}MN^{2})$ to
$\mathcal{O}\big{(}N_{b}MN^{2}/N_{p})$; where $N_{b}$ and $N_{p}$ are the
number of baselines and the number of computing nodes, respectively. The
implementation of this algorithm is not part of this work, we leave it for
future work.

The rest of the paper is organised as follows: Section 2 gives an overview of
radio interferometric measurements and the theory behind SVD. Section 3
describes the mathematical formalism behind the algorithms. Section 4 analyses
the effects of the proposed algorithms on the image and the noise. Simulation
results are discussed in Section 5 and Section 6 draws conclusions and
suggests future works. Our mathematical notations are summarised in Table 2.

Methods |  Compression |  Preserved spectral and temporal resolution within the FoV? |  Speedup processing?
---|---|---|---
Traditional averaging |  yes |  no |  yes
BDA |  yes |  no |  yes
Simple SVD |  yes |  yes |  no
BDSVD |  yes |  yes |  no
Table 1: The use cases of traditional averaging and BDA vs. simple SVD and
BDSVD

##  2 Background

In this section, we discuss the mathematical framework of the visibility data
measurements provided by the radio interferometric measurement equation (RIME)
formalism (Hamaker et al., 1996; Smirnov, 2011). The framework is used to
understand, describe and measure the performance of simple SVD and BDSVD.

###  2.1 Full-sky RIME and traditional averaging

Under the RIME formalism and for a single point source, Smirnov (2011)
establishes the relationship between a measured visibility at time $t$ and
frequency $\nu$ and the averaged visibility over some time and frequency
integration interval $[t_{0},t_{1}]\times[\nu_{0},\nu_{1}]$. We extend this
relationship for a full-sky RIME formalism:

$\displaystyle\mathcal{V}(\boldsymbol{u}_{pq}(t,\nu))$
$\displaystyle=\sum_{\boldsymbol{l}}\big{(}\mathcal{G}_{p\boldsymbol{l}t\nu}\mathcal{I}_{\boldsymbol{l}}\mathcal{G}_{q\boldsymbol{l}t\nu}^{\dagger}\big{)}\mathcal{K}^{\boldsymbol{l}}_{pqt\nu}.$
(2)

The scalar term
$\mathcal{K}^{\boldsymbol{l}}_{pqt\nu}=\text{e}^{-2i\pi(\boldsymbol{l}-\boldsymbol{l}_{0})\boldsymbol{u}_{pqt\nu}}$
describes the effects of the position of antennas $p,q$ and the rotation of
the baseline vector $\boldsymbol{u}_{pq}(t,\nu)\equiv\boldsymbol{u}_{pqt\nu}$
which tracks a source located in the direction of the unit vector
$\boldsymbol{l}$. A compensating delay $\boldsymbol{l}_{0}$ is introduced by
the correlator to enforce that $\mathcal{K}^{\boldsymbol{l}}_{pqt\nu}\equiv 1$
at the phase centre of the observation:

$\displaystyle\boldsymbol{u}_{pqt\nu}=\frac{\nu}{c}\begin{bmatrix}u_{pqt}\\\
v_{pqt}\\\ w_{pqt}\end{bmatrix},\boldsymbol{l}=\begin{bmatrix}l\\\ m\\\
n\end{bmatrix},\boldsymbol{l}_{0}=\begin{bmatrix}0\\\ 0\\\ 1\end{bmatrix}.$

In the above equations, † is the complex transpose operator, $c$ stands for
the speed of light, $\mathcal{I}_{\boldsymbol{l}}$ is the sky brightness,
$\mathcal{G}_{p\boldsymbol{l}t\nu}$ and
$\mathcal{G}_{q\boldsymbol{l}t\nu}^{\dagger}$ group the direction-dependent
Jones matrices of antenna $p$ and $q$, respectively.

Below is an approximate of the average of Eq. 2 over
$[t_{0},t_{1}]\times[\nu_{0},\nu_{1}]$:

$\displaystyle V_{pq(t\nu)_{i}}$
$\displaystyle\simeq\sum_{\boldsymbol{l}}\text{sinc}\frac{\Delta\Phi}{2}\text{sinc}\frac{\Delta\psi}{2}\big{(}\mathcal{G}_{p\boldsymbol{l}(t\nu)_{i}}\mathcal{I}_{\boldsymbol{l}}\mathcal{G}_{q\boldsymbol{l}(t\nu)_{i}}^{\dagger}\big{)}\mathcal{K}^{\boldsymbol{l}}_{pq(t\nu)_{i}},$
(3)

where $(t\nu)_{i}=t_{i}\nu_{i}$ with $t_{i}=(t_{1}+t_{0})/2$ and
$\nu_{i}=(\nu_{1}+\nu_{0})/2$. The phase differences, $\Delta\Phi$ and
$\Delta\psi$ are defined as:

$\displaystyle\Delta\Phi$
$\displaystyle=\text{arg}\mathcal{G}_{p\boldsymbol{l}t_{1}\nu_{i}}+\text{arg}\mathcal{G}_{q\boldsymbol{l}t_{1}\nu_{i}}^{\dagger}+\text{arg}\mathcal{K}^{\boldsymbol{l}}_{pqt_{1}\nu_{i}}$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}-\Big{(}\text{arg}\mathcal{G}_{p\boldsymbol{l}t_{0}\nu_{i}}+\text{arg}\mathcal{G}_{q\boldsymbol{l}t_{0}\nu_{i}}^{\dagger}+\text{arg}\mathcal{K}^{\boldsymbol{l}}_{pqt_{0}\nu_{i}}\Big{)}$
$\displaystyle\Delta\psi$
$\displaystyle=\text{arg}\mathcal{G}_{p\boldsymbol{l}t_{i}\nu_{1}}+\text{arg}\mathcal{G}_{q\boldsymbol{l}t_{i}\nu_{1}}^{\dagger}+\text{arg}\mathcal{K}^{\boldsymbol{l}}_{pqt_{i}\nu_{1}}$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}-\Big{(}\text{arg}\mathcal{G}_{p\boldsymbol{l}t_{i}\nu_{0}}+\text{arg}\mathcal{G}_{q\boldsymbol{l}t_{i}\nu_{0}}^{\dagger}+\text{arg}\mathcal{K}^{\boldsymbol{l}}_{pqt_{i}\nu_{0}}\Big{)}.$

The average in Eq. 3 always results in a net loss of amplitude for an off-
phase centre source; in practice the conditioning $t_{1}-t_{0}\rightarrow 0$
and $\nu_{1}-\nu_{0}\rightarrow 0$ cannot be satisfied because the phase terms
in $\mathcal{K}^{\boldsymbol{l}}_{pqt\nu},\mathcal{G}_{p\boldsymbol{l}t\nu}$
and $\mathcal{G}_{q\boldsymbol{l}t\nu}^{\dagger}$ vary over time and
frequency. This loss of amplitude is known as time and frequency decorrelation
(or smearing) when it is caused only by the scalar
$\mathcal{K}^{\boldsymbol{l}}_{pqt\nu}$. The general term decoherence is used
when the amplitude loss is caused by the cumulative phase variation of
$\mathcal{K}^{\boldsymbol{l}}_{pqt\nu},\mathcal{G}_{p\boldsymbol{l}t\nu},\mathcal{G}_{q\boldsymbol{l}t\nu}^{\dagger}$
since $\mathcal{G}_{p\boldsymbol{l}t\nu}$ and
$\mathcal{G}_{q\boldsymbol{l}t\nu}^{\dagger}$ also hold complex phases which
vary in time and frequency. Eq. 3 shows that for a full sky RIME, the
decoherence can be measured separately as the sum of the contribution of each
individual source. When considering the decoherence introduced by all the
complex terms, an approximation to the average of each individual source that
is part of the sum in Eq. 3 can be expressed as the phase differences. In the
decoherence formulation discussed above, we have neglected the case of
decoherence caused by $\mathcal{I}_{\boldsymbol{l}}$ because this work focuses
on the decoherence of point sources rather than extended sources for which the
phase term in $\mathcal{I}_{\boldsymbol{l}}$ varies considerably over time and
frequency.

Notation |  Description
---|---
$\boldsymbol{V}$ |  The data matrix for $N_{b}$ baselines before compression
$\boldsymbol{V}_{n}$ |  The data matrix for $N_{b}$ baselines after compression
$\mathcal{V}$ |  The unsampled visibility data
$V_{pqt\nu}$ |  The sampled visibility data for $pq$ at $t,\nu$
$\boldsymbol{u}_{pq}$ |  The baseline $pq$ vector
$\mathcal{I}_{\boldsymbol{l}}$ |  The sky brightness in the direction of $\boldsymbol{l}$
$\boldsymbol{V}_{pq}$ |  The rank$-r$ data matrix for baseline $pq$
$\boldsymbol{V}_{pq,n}$ |  The rank$-n$ approximation of $\boldsymbol{V}_{pq}$
$\eta_{pqk}$ |  The singular value
$CF$ |  The compression factor
$N_{b}$ |  The number of baselines
$M,N$ |  The size of matrices
$\epsilon$ |  The maximum threshold error
$\varepsilon$ |  The minimum percentage signal to preserve
$\boldsymbol{I}^{\mathrm{d}}_{\mathrm{pq}}$ |  The dirty image for baseline $pq$
$\widetilde{\boldsymbol{I}}^{\mathrm{d}}_{\mathrm{pq}}$ |  The compressed dirty image for baseline $pq$
$\boldsymbol{I}^{\mathrm{d}}$ |  The dirty image matrix
$\widetilde{\boldsymbol{I}}^{\mathrm{d}}$ |  The compressed dirty image matrix
$\boldsymbol{I}_{\mathrm{loss}}^{\mathrm{d}}$ |  The loss image matrix
$\mathcal{F}$ |  The Fourier transform operator
$S/R$ |  The signal to noise
Table 2: Main mathematical notations: lowercase bold letters are vectors and
uppercase bold letters are matrices. Calligraphic capital letters are used to
designate functional forms. Everything else is a constant.

###  2.2 BDA and implementation

To aggressively compress visibility data and mitigate amplitude lost caused by
traditional averaging, several authors have discussed the potential of using
BDA (Cotton, 2009; Atemkeng et al., 2018; Wijnholds et al., 2018). As shown in
Eq. 3 the product of the three complex terms
$\mathcal{K}^{\boldsymbol{l}}_{pqt\nu},\mathcal{G}_{p\boldsymbol{l}t\nu}$ and
$\mathcal{G}_{q\boldsymbol{l}t\nu}^{\dagger}$ is attenuated by $sinc$
functions. The degree of attenuation at the edges of the FoV is determined by
the width of these $sinc$ functions which depend on the spatial frequencies
that each baseline samples over $[t_{0},t_{1}]\times[\nu_{0},\nu_{1}]$. On
long spacing, traditional averaging with a wide
$[t_{0},t_{1}]\times[\nu_{0},\nu_{1}]$ will result in a narrow $sinc$ and
significant drop in source amplitude. In order to keep amplitude loss equal
across all baselines for a given FoV, the widths of the $sinc$ functions must
remain constant, so that the time-frequency interval over which the data is
averaged varies. This method is substantial for dense-core interferometric
arrays where data is aggressively averaged on shorter baselines. A small time-
frequency interval is required to reach the width limit of the $sinc$
functions on longer baselines, resulting in less compression.

It is advantageous for a BDA implementation to integrate with existing
software, formats and specifications to benefit from compression, and for
testing with real observational data. The contemporary specification for radio
astronomy data is the MS v2, based on the CASA Table Data System (CTDS)
format. We have, therefore, developed a BDA implementation targeting this
specification and format.

MS rows are grouped by baseline, sorted by time and aggregated into bins whose
$\text{sinc}(\Delta\Phi/2)$ does not exceed the acceptable decorrelation
tolerance defined by $\mathcal{K}^{\boldsymbol{l}}_{pqt\nu}$. Each bin’s
$\text{sinc}(\Delta\Phi/2)$ is inverted, firstly to obtain
$\text{sinc}(\Delta\psi/2)$ and secondly, the bin frequency integration
interval $[\nu_{0},\nu_{1}]$. Practically, each bin can have a different
frequency interval and this channelisation does not conveniently fit into the
CTDS format which has a fixed number of channels per spectral window.

Therefore, to fit BDA data into this format, we discretise the channelisation
by subdividing the original spectral window by the integral factors of its
channels to produce new spectral windows representing each discretisation.
Then, each BDA bin is output to a single row in the output MS and each row is
associated with a different spectral window, thereby trading a small factor of
compression for compatibility with the MS.

###  2.3 Low-rank approximation: SVD

The SVD break down a given complex data matrix $\boldsymbol{V}$ having $M$
rows and $N$ columns into three independent matrices:

$\displaystyle\boldsymbol{V}=\boldsymbol{A}\boldsymbol{\Lambda}\boldsymbol{C}^{\dagger},$
(4)

where $\boldsymbol{\Lambda}$ is a diagonal matrix of size $M\times N$,
$\boldsymbol{A}$ of size $M\times M$ and $\boldsymbol{C}$ of size $N\times N$
are unitary matrices. These matrices are defined as:

$\displaystyle\boldsymbol{\Lambda}=\begin{bmatrix}\boldsymbol{\Lambda}_{r}&\boldsymbol{0}\\\
\boldsymbol{0}&\boldsymbol{0}\end{bmatrix},~{}~{}\boldsymbol{A}=[\boldsymbol{A}_{r},\boldsymbol{a}_{r+1},\cdots,\boldsymbol{a}_{M}],~{}~{}\boldsymbol{C}=\begin{bmatrix}\boldsymbol{C}_{r}\\\
\boldsymbol{c}_{r+1}\\\ \vdots\\\ \boldsymbol{c}_{N}\end{bmatrix}.$ (5)

The decomposition does not require $\boldsymbol{V}$ to be a square matrix.
Matrix
$\boldsymbol{\Lambda}_{r}=\mathrm{diag}(\eta_{1},\eta_{2},\cdots,\eta_{r})$ of
size $r\times r$ is diagonal with $r=\min(M,N)$ and $\eta_{k}$ are the
singular values of $\boldsymbol{V}$. As discussed in Stewart (1998),
$\eta_{k}=\sqrt{\lambda_{k}}$ where $\lambda_{k}$ are the eigenvalues of
$\boldsymbol{V}^{\dagger}\boldsymbol{V}$. We note that all the entries in the
diagonal of $\boldsymbol{\Lambda}_{r}$ are all non-zeros and are in decreasing
order. Here,
$\boldsymbol{A}_{r}=[\boldsymbol{a}_{1},\boldsymbol{a}_{2},\cdots,\boldsymbol{a}_{r}]$
and
$\boldsymbol{C}_{r}=[\boldsymbol{c}_{1},\boldsymbol{c}_{2},\cdots,\boldsymbol{c}_{r}]^{\mathrm{T}}$,
where T is the transpose operator. The vectors, $\boldsymbol{a}_{k}$ of size
$M$ and $\boldsymbol{c}_{k}$ of size $N$ are the left and right singular
vectors of $\boldsymbol{V}$, respectively. $\boldsymbol{a}_{k}$ and
$\boldsymbol{c}_{k}$ are eigenvectors of
$\boldsymbol{V}\boldsymbol{V}^{\dagger}$ and
$\boldsymbol{V}^{\dagger}\boldsymbol{V}$, respectively. Eq. 4 remains valid if
rewritten as:

$\displaystyle\boldsymbol{V}$
$\displaystyle=\boldsymbol{A}_{r}\boldsymbol{\Lambda}_{r}\boldsymbol{C}_{r}^{\dagger}$
(6)
$\displaystyle=\sum_{k=1}^{r}\eta_{k}\boldsymbol{a}_{k}\boldsymbol{c}_{k}^{\dagger}.$
(7)

Each $\eta_{k}$ quantifies how best the corresponding component
$\eta_{k}\boldsymbol{a}_{k}\boldsymbol{c}_{k}^{\dagger}$ contributes to the
relevant features in Eq. 7; such as in the reconstruction of $\boldsymbol{V}$.
The larger $\eta_{k}$ is, the more
$\eta_{k}\boldsymbol{a}_{k}\boldsymbol{c}_{k}^{\dagger}$ effectively
contributes to the reconstruction of $\boldsymbol{V}$. Note that
$\boldsymbol{a}_{k}$ and $\boldsymbol{c}_{k}^{\dagger}$ determine the geometry
of the features contained in
$\eta_{k}\boldsymbol{a}_{k}\boldsymbol{c}_{k}^{\dagger}$.

##  3 Proposed lossy compression methods

The low-rank approximation methods described in this section are applied to
the raw visibility data (continuous data), and at each baseline separately. We
note that the methods can also be applied to gridded visibility data. However,
the gridded data lie on a regular grid where data for all baselines have been
interpolated together, making it difficult to gauge acceptable thresholding of
the nonzero singular values. Also, since each baseline observes a different
portion of the sky, the signal distortion and attenuation are baseline-
dependent and the noise variance is different per-visibility data. These
effects cannot be taken into account in gridded visibility data. It is
insufficient to store the images (therefore gridded visibility data) because
typically calibration problems dominate and data has, out of necessity, be
recalibrated most of the time to create science ready products for specific
continuum, transient and line subdomains. This cannot be done well with all
the baselines being convolved onto grid spacing’s.

Following the above recommendations, we formally express the visibility data
(for a given baseline) in a compact and robust matrices formalism. Assume a
non-polarized sky, a single channel timeslot visibility measured by a baseline
$pq$ is the complex Stokes value:

$\displaystyle V_{pqt_{i}\nu_{j}}$
$\displaystyle=\delta_{pqij}\big{(}\mathcal{V}+\mathcal{E}\big{)},$ (8)

where $\delta_{pqij}=\delta(\boldsymbol{u}-\boldsymbol{u}_{pqt_{i}\nu_{j}})$
is a delta-function shifted to the $uv$-point being sampled, $\mathcal{V}$ is
discussed in Eq. 2 and $\mathcal{E}$ is the contaminating random noise with
zero mean and r.m.s. $\Sigma$. Assuming an observation of $M$ timeslots and
$N$ channels, we can package the visibility data for single baseline $pq$ into
a single data matrix, $\boldsymbol{V}_{pq}$ of size $M\times N$:

$\boldsymbol{V}_{pq}=\begin{bmatrix}V_{pqt_{1}\nu_{1}}&\cdots&V_{pqt_{1}\nu_{N}}\\\
\vdots&\cdots&\vdots\\\ V_{pqt_{M}\nu_{1}}&\cdots&V_{pqt_{M}\nu_{N}}\\\
\end{bmatrix}.$ (9)

The formalism can be extended to data structure with orthogonal bases capable
of fully describing the sky. We discuss the compression of
$\boldsymbol{V}_{pq}$ and then, in Section 4, we provide a detailed analysis
of the compression effect on $\Sigma$ and the sky image. Note that it is
possible to remodel $\boldsymbol{V}_{pq}$ into sub-matrices and independently
find the low-rank of each sub-matrix. The latter should be considered in cases
where computation is a bottleneck (i.e. $MN\rightarrow\infty$) and where
approximating local low-rank matrices is an advantage (Lee et al., 2016).

###  3.1 Method 1: Simple SVD

Finding a very small low-rank data matrix of $\boldsymbol{V}_{pq}$ without
information loss is challenging because the entries in $\boldsymbol{V}_{pq}$
are strongly correlated. Uncorrelated $\boldsymbol{V}_{pq}$ is rarely possible
in a real observation due to system electronics, and the fact that each
visibility data point comes from a pairwise correlation between antenna
voltages, etc. This implies that there exists no low-rank approximation,
$\boldsymbol{V}_{pq,n}$ of rank$-n$ that can perfectly reconstruct
$\boldsymbol{V}_{pq}$ of rank$-r$:

$\displaystyle\boldsymbol{V}_{pq,n}$
$\displaystyle\simeq\boldsymbol{V}_{pq},~{}1\leq n<r.$ (10)
$\displaystyle\boldsymbol{V}_{pq,n}$
$\displaystyle=\mathbf{A}_{pq,n}\mathbf{\Lambda}_{pq,n}\mathbf{C}^{\dagger}_{pq,n}$
(11)
$\displaystyle=\sum_{k=1}^{n}\eta_{pqk}\boldsymbol{a}_{pqk}\boldsymbol{c}_{pqk}^{\dagger}.$
(12)

For simplicity, Eq. 11 can also be written as a Kronecker product, $\otimes$
after vectorisation:

$\displaystyle\mbox{vec}(\boldsymbol{V}_{pq,n})$
$\displaystyle=(\mathbf{C}^{\dagger}_{pq,n}\otimes\mathbf{A}_{pq,n})\boldsymbol{\eta}_{pq},$
(13)

where
$\boldsymbol{\eta}_{pq}=[\eta_{pq1},\eta_{pq2},\cdots,\eta_{pqn}]^{\mathrm{T}}$.

This simple SVD only allows an equal compression factor across all the
baselines as $n$ is fixed. The compression factor is the ratio between the
number of data points of the uncompressed data and the number of data points
of the compressed data. The entries of $\mathbf{A}_{pq,n}$ and
$\mathbf{C}^{\dagger}_{pq,n}$ are complex numbers, while the entries of
$\mathbf{\Lambda}_{pq,n}$ are real numbers. In terms of computer storage
requirements, if a complex entry counts as one entry, a real entry should
count as $0.5$. Therefore, if the size of $\boldsymbol{V}_{pq}$ is $M\times N$
then one can show that the number of elements in the sub-matrices needed to
compute $\boldsymbol{V}_{pq,n}$ is $n(M+N+0.5)$ which leads to an overall
compression factor, $CF$ of:

$CF=\frac{MN}{n(M+N+0.5)}.$ (14)

In this setting, one has to choose $CF$ carefully to determine the number of
singular values to retain on all the baselines:

$n=\big{\lceil}\frac{MN}{CF(M+N+0.5)}\big{\rceil},$ (15)

where $\lceil.\rceil$ is the ceiling operator. The compression loss is
computed from the tensor norm:

$\displaystyle||\boldsymbol{V}-\boldsymbol{V}_{n}||$
$\displaystyle\leq\epsilon||\boldsymbol{V}||,$ (16)

where $\epsilon$ is the maximum threshold error provided to control the
divergence between the visibility data tensors $\boldsymbol{V}$ and
$\boldsymbol{V}_{n}$ representing the data for the $N_{b}$ baselines before
and after compression, respectively. In this work, we defined a tensor norm as
the sum of the Frobenius norm, $||.||_{\mathrm{F}}$ per baseline:

$\displaystyle||\boldsymbol{V}-\boldsymbol{V}_{n}||$
$\displaystyle=\sum_{pq}||\boldsymbol{V}_{pq}-\boldsymbol{V}_{pq,n}||_{\mathrm{F}}.$
(17)

Eq. 17 can be computed from the sum of the Euclidean norm of the singular
values not retained at each baseline:

$\displaystyle||\boldsymbol{V}-\boldsymbol{V}_{n}||$
$\displaystyle=\sum_{pq}\sqrt{\sum_{k=n+1}^{r}\eta_{pqk}^{2}}.$ (18)

Note that in Eq 18, $k$ runs from $n+1$ to $r$ i.e. $(r-n)$ singular values
are not retained in the compression. $||\boldsymbol{V}||$ is defined as:

$\displaystyle||\boldsymbol{V}||$
$\displaystyle=\sum_{pq}||\boldsymbol{V}_{pq}||_{\mathrm{F}}$ (19)
$\displaystyle=\sum_{pq}\sqrt{\sum_{k=1}^{r}\eta_{pqk}^{2}}.$ (20)

Instead of using a threshold error, a minimum percentage of the signal to
preserve, $\varepsilon$ can be specified:

$\displaystyle||\boldsymbol{V}_{n}||$
$\displaystyle\geq\varepsilon||\boldsymbol{V}||,$ (21)
$\displaystyle||\boldsymbol{V}_{n}||$
$\displaystyle=\sum_{pq}||\boldsymbol{V}_{pq,n}||_{\mathrm{F}}$ (22)
$\displaystyle=\sum_{pq}\sqrt{\sum_{k=1}^{n}\eta_{pqk}^{2}}.$ (23)

If the constraint $\epsilon$ or $\varepsilon$ is given then it is
computationally cheap to obtain the corresponding value of $n$ that satisfies
the constraint as shown in Algorithm 1. Note that these norms could directly
be obtained from computing the Euclidean norm of
$\mbox{vec}(\boldsymbol{V}_{pq,n})$, $\mbox{vec}(\boldsymbol{V})$ and
$\mbox{vec}(\boldsymbol{V}_{n})$; the vectorisation versions of
$\boldsymbol{V}_{pq,n},\boldsymbol{V}$ and $\boldsymbol{V}_{n}$, respectively.

Data: $\boldsymbol{V};\epsilon~{}\text{or}~{}\varepsilon$

Result: the number of singular values to retain, $n$

1 $n\leftarrow 0;\boldsymbol{V}_{pq,n}\leftarrow\boldsymbol{0}$

2 if  _$\epsilon$ is given_ then

3 do

4 for  _all baseline $pq$_ do

5
$\boldsymbol{V}_{pq,n}\leftarrow\boldsymbol{V}_{pq,n-1}+\eta_{pqn}\boldsymbol{a}_{pqn}\boldsymbol{c}_{pqn}^{\dagger}$

6

7 end for

8 $n\leftarrow n+1$

9 while
_$\frac{\sum_{pq}||\boldsymbol{V}_{pq}-\boldsymbol{V}_{pq,n}||_{\mathrm{F}}}{\sum_{pq}||\boldsymbol{V}_{pq}||_{\mathrm{F}}}
>\epsilon$ and $n\leq r$_;

10else

11 do

12 for  _all baseline $pq$_ do

13
$\boldsymbol{V}_{pq,n}\leftarrow\boldsymbol{V}_{pq,n-1}+\eta_{pqn}\boldsymbol{a}_{pqn}\boldsymbol{c}_{pqn}^{\dagger}$

14

15 end for

16 $n\leftarrow n+1$

17 while
_$\frac{\sum_{pq}||\boldsymbol{V}_{pq,n}||_{\mathrm{F}}}{\sum_{pq}||\boldsymbol{V}_{pq}||_{\mathrm{F}}}\times
100 <\varepsilon$ and $n\leq r$_;

18 end if

Algorithm 1 Finding $n$ for simple SVD. Here, all the $\boldsymbol{V}_{pq}$
are taken from $\boldsymbol{V}$

###  3.2 Source position in the sky and baselines effects on the distribution
of singular values

As discussed in Section 2.1 the baselines observe different portions of the
sky due to the variation of the phase terms in
$\mathcal{K}^{\boldsymbol{l}}_{pqt\nu},\mathcal{G}_{p\boldsymbol{l}t\nu}$ and
$\mathcal{G}_{q\boldsymbol{l}t\nu}^{\dagger}$, as well as the different
Nyquist sampling due to the baselines geometry and distribution. With this
information, we noted that decorrelation is baseline-dependent. Given this
finding, the following questions arise: Is the distribution of the singular
values equal across all baselines? If not, how do the strength and decay of
the singular values differ at each baseline? Can the variation in strength and
decay of the singular values be used to aggressively compress the data and
what limits the compression?

To answer the above questions, Figure 1 shows the first $20$ singular values
of three baselines; the longest baseline (left-column panels), medium-length
baseline (middle-column panels) and shortest (right-column panels) of a 1 Jy
source at the phase centre (top-row panels), 5 deg (middle-row panels) and 10
deg (bottom-row panels) away from the phase centre. The singular values are
obtained by simulating the MeerKAT telescope at 1.4 GHz. The simulation
corresponds to $10~{}000$ timeslots per second with a total bandwidth of 0.8
MHz divided into 10 channels each of width 80 kHz. This figure shows some
important behaviors of the distribution of the singular values across
baselines, therefore, it deserves a detailed explanation. Regardless of the
position of the source, the first component in each baseline captures the most
features in the original data. In addition, the strength and decay of the
singular values are different at each baseline and source position, with
shorter baselines decaying rapidly compared to longer baselines where the
decay is slow. At the phase centre of the observation, the singular values
decay faster; the degree of decay becomes slower as the source moves away from
the phase centre. In other words:

  * •

For a source at the phase centre, the first component in all the baselines
captures the entire features in the original data, in other words $\eta_{pq1}$
is very large and $\forall k>1,\eta_{pqk}\sim 0$. This result aligns with what
is known; for all baselines, decorrelation is negligible at the phase centre
of the observation.

  * •

The spread of the singular values is baseline-dependent; the singular values
are spread to several components on the longer baselines while for shorter
baselines only a few of the first components capture features; for example,
say $pq$ and $\alpha\beta$ are two baselines with length
$||\textbf{u}_{pq}||_{2}<||\textbf{u}_{\alpha\beta}||_{2};\exists i$ such that
$\forall k>i$:

$\displaystyle\eta_{pqi}$ $\displaystyle\gg\eta_{\alpha\beta
i}~{}\text{and}~{}\eta_{\alpha\beta k}\gg\eta_{pqk},$ (24)

where $||.||_{2}$ is the Euclidean norm. This result aligns with the fact that
longer baselines observed severe amplitude loss compared to short baselines.

  * •

The spread of the singular values is position-dependent; the singular values
are spread in many components for far-field sources compared to nearby field
sources. This is easy to explain - far-field sources experience severe
amplitude loss compared to nearby field sources.

The above results show that if we were to consider the same number of $n$
components to retain (i.e. $r-n$ components to discard) across all the
baselines as described in Section 3.1, this will have a direct impact on
baseline compression errors,
$||\boldsymbol{V}_{pq}-\boldsymbol{V}_{pq,n}||_{\mathrm{F}}$, which will be
small on short baselines compared to long baselines. Figure 2 shows these
compression errors from the simulation in Figure 1, where the first eigh
components of each of the baselines are retained. In this figure, the
compression errors are plotted against baseline lengths; left-, middle-, and
right panels are the compression errors when the source is at the phase
centre, $5$ deg and $10$ deg, respectively. The following is observed:

  * •

At the phase centre,
$||\boldsymbol{V}_{pq}-\boldsymbol{V}_{pq,8}||_{\mathrm{F}}\sim 0$ for all
baselines $pq$.

  * •

On longer baselines,
$||\boldsymbol{V}_{pq}-\boldsymbol{V}_{pq,8}||_{\mathrm{F}}$ increased rapidly
for far-field sources.

  * •

For any two baselines of mismatched length, say $pq$ and $\alpha\beta$ with
$||\textbf{u}_{pq}||_{2}<||\textbf{u}_{\alpha\beta}||_{2}$ we have:

$\displaystyle||\boldsymbol{V}_{pq}-\boldsymbol{V}_{pq,8}||_{\mathrm{F}}<||\boldsymbol{V}_{\alpha\beta}-\boldsymbol{V}_{\alpha\beta,8}||_{\mathrm{F}}.$
(25)

We note that the singular values experience a slow decay for out-of-phase
centre sources and on long baselines. Since the strength of a singular value,
$\eta_{pqk}$ indicates the degree of features that the component,
$\eta_{pqk}\boldsymbol{a}_{pqk}\boldsymbol{c}_{pqk}^{\dagger}$ contributes to
the original data for baseline $pq$, depending on how the singular values are
decaying for the given baseline, the components with smaller singular values
could be discarded. By so doing, we observe the following:

  * •

A significant number of components could be discarded on shorter baselines
with little or no loss of source amplitude. This leads to aggressive data
compression with negligible effect on image fidelity.

  * •

While attempting to preserve the loss of source amplitude and image fidelity,
only a few components could be discarded on longer baselines.

This makes the compression baseline-dependent. However, for a given FoV, a
sophisticated compression approach is to carefully choose $n$ on each of the
baselines so that the errors in Eq. 25 are equal on all the baselines. This
method is described in Section 3.3.

Figure 1: The first 20 singular values of three baselines; the longest
baseline (left-column panels), medium-length baseline (middle-column panels)
and shortest (right-column panels) of a 1 Jy source at the phase centre (top-
row panels), 5 deg (middle-row panels) and 10 deg (bottom-row panels) away
from the phase centre. The singular values are obtained from simulating the
MeerKAT telescope at 1.4 GHz. The data is sampled at 1 s and 80 kHz during
$166$ min $40$ s with 0.8 MHz bandwidth.

Figure 2: Compression errors from the simulation in Figure 1, where the first
eight components of each of the baselines are retained. At the phase centre,
the compression error tends to zero for all baselines. On long baselines, the
compression error increased rapidly for far-field sources. The compression
error is small on short baselines compared to long baselines.

###  3.3 Method 2: Baseline-dependent SVD (BDSVD)

The BDSVD method finds the corresponding baseline-dependent number of
components to retain on each baseline so that the compression errors between
baselines do not vary, for example for any two baselines $pq$ and
$\alpha\beta$ with $||\textbf{u}_{pq}||_{2}<||\textbf{u}_{\alpha\beta}||_{2}$
find the number of components $n_{pq}$ and $n_{\alpha\beta}$, respectively so
that:

$\displaystyle||\boldsymbol{V}_{pq}-\boldsymbol{V}_{pq,n_{pq}}||_{\mathrm{F}}\equiv||\boldsymbol{V}_{\alpha\beta}-\boldsymbol{V}_{\alpha\beta,n_{\alpha\beta}}||_{\mathrm{F}},$
(26)

where

$\displaystyle\boldsymbol{V}_{pq,n_{pq}}$
$\displaystyle=\sum_{k=1}^{n_{pq}}\eta_{pqk}\boldsymbol{a}_{pqk}\boldsymbol{c}_{pqk}^{\dagger}.$
(27)

It should be noted that the maximum threshold error $\epsilon$ (or minimum
percentage of the signal to preserve $\varepsilon$) is now verified on each
individual baseline as opposed to Eq. 16 (or to Eq. 21):

$\displaystyle||\boldsymbol{V}_{pq}-\boldsymbol{V}_{pq,n_{pq}}||_{\mathrm{F}}$
$\displaystyle\leq\epsilon||\boldsymbol{V}_{pq}||_{\mathrm{F}}$ (28)
$\displaystyle||\boldsymbol{V}_{pq,n_{pq}}||_{\mathrm{F}}\geq\varepsilon||\boldsymbol{V}_{pq}||_{\mathrm{F}}.$
(29)

Intuitively, since $||\textbf{u}_{pq}||_{2}<||\textbf{u}_{\alpha\beta}||_{2}$
this means that $n_{pq}<n_{\alpha\beta}$ while the cumulative strength of the
singular values of a compact source are equal on all baselines,

$\displaystyle\sum_{k=1}^{n_{pq}}\eta_{pqk}^{2}=\sum_{k=1}^{n_{\alpha\beta}}\eta_{\alpha\beta
k}^{2}.$ (30)

With this information, Eq. 32 is correct and valid.

$\displaystyle||\boldsymbol{V}_{n}||$
$\displaystyle=\sum_{pq}\sqrt{\sum_{k=1}^{n_{pq}}\eta_{pqk}^{2}}$ (31)
$\displaystyle=N_{b}\sqrt{\sum_{k=1}^{n_{pq}}\eta_{pqk}^{2}}.$ (32)

In the above, the number of components to retain is baseline-dependent, thus
the compression factor also becomes baseline-dependent:

$CF_{pq}=\frac{MN}{n_{pq}(M+N+0.5)},$ (33)

which is larger for shorter baselines compared to longer baselines, for
example $CF_{pq}>CF_{\alpha\beta}$ since $n_{pq}<n_{\alpha\beta}$. As opposed
to Eq. 14, the BDSVD overall $CF$ follows:

$CF=\frac{N_{b}MN}{(M+N+0.5)\sum_{pq}n_{pq}}.$ (34)

Algorithm 2 describes an iterative process to find the baseline-dependent
number of component, $n_{pq}$ given $\epsilon$ or $\varepsilon$.

For compression factor $CF$, the space-saving $SS$ is measured as:

$SS=(1-CF^{-1})\times 100\%.$ (35)

###  3.4 Computation complexity and parallelisation

Low-rank approximation has computational drawbacks. It is expensive to compute
the SVD of a very large matrix. For visibility data compression as discussed
in this paper, for some cases (e.g. to archive the data), we can overlook this
computational drawback since we have to compress the data only once at a time.
In other cases, for example big data radio interferometers such as the SKA,
the amount of visibility data on each of the baselines is very large, even the
once-off compression can become a challenge since finding the exact $CF_{pq}$
involves the computation of singular values at each baseline sequentially,
which is an expensive computational task. For example, if the visibility data
for each baseline has size $M\times N$ with $M\geq N$, then for $N_{b}$
baselines, the sequential full rank$-n$ simple SVD scales as:

$s_{\mathrm{cost}}\sim\mathcal{O}\big{(}N_{b}MN^{2}\big{)}.$ (36)

BDSVD scales closely with simple SVD for equal $CF$. However, the good news is
that visibility data has some natural basis consisting of row data. This
information can be used to speed up the compression via a parallel algorithm.
The row visibility data for a given baseline can be shared with multiple
compute nodes; for example, for each baseline, we can subdivide the row data
into chunks of $\boldsymbol{V}_{pq}^{\mathrm{[i]}}$ and compute the SVD for
all the chunks in parallel. The subscript [i] indicates the $i^{th}$ chunk of
the row visibility data. We briefly discuss this distributed compression
process in Algorithm 3, however, further investigation is needed to assess its
computational efficiency in practice and to determine at what rate of
visibility data size the algorithm should be enforced.

Algorithm 3 shows that we can compute the SVD for each of the baselines in
parallel (line 1); assume $N_{p}$ parallel nodes. Also (as shown in line 5)
for each baseline, to find $\mathbf{A}_{pq},\mathbf{\Lambda}_{pq}$ and
$\mathbf{C}_{pq}$ involve computing all
$\boldsymbol{V}_{pq}^{\mathrm{[i]\dagger}}\boldsymbol{V}_{pq}^{[i]}$ in
parallel. $\boldsymbol{V}_{pq}^{\dagger}\boldsymbol{V}_{pq}$ is then obtained
from summing all
$\boldsymbol{V}_{pq}^{\mathrm{[i]\dagger}}\boldsymbol{V}_{pq}^{[i]}$. The
latter suggests that the complexity of Algorithm 3 scales as:

$\displaystyle p_{\mathrm{cost}}$
$\displaystyle\sim\mathcal{O}\bigg{(}N_{b}MN^{2}/N_{p}\bigg{)}.$ (37)

Ideally, we would hope to demonstrate strong scaling over Eq. 36, but this
remains outside the scope of this work.

##  4 Effects on the image and noise: analytical quantification

The SVD is not perfectly orthogonal in practice, both noise and signal are
present in each of the components of the decomposition. A trade-off between
the amount of signal and noise to remove is necessary. Simple SVD and BDSVD
also result in a small loss of signal and retain a small amount of noise. Our
goal now is to provide the mathematical models to quantify the signal loss and
the noise removed as a form of contribution at all baselines.

###  4.1 Effect on the image

As the Fourier transform $\mathcal{F}$ is linear, the compressed dirty image
$\widetilde{\boldsymbol{I}}^{\mathrm{d}}_{\mathrm{pq}}$ of a single baseline
is the sum:

$\displaystyle\widetilde{\boldsymbol{I}}^{\mathrm{d}}_{\mathrm{pq}}$
$\displaystyle=\mathcal{F}\boldsymbol{V}_{pq,n_{pq}}$ (38)
$\displaystyle=\sum_{k=1}^{n_{pq}}\eta_{pqk}\boldsymbol{I}_{pqk}^{\mathrm{d}},$
(39)

where $\eta_{pqk}\boldsymbol{I}_{pqk}^{\mathrm{d}}$ is the $k^{th}$ linearly
independent component of the sky image seen by the baseline $pq$ and
$\mathcal{F}$ is a unitary linear operator. In relation to Parseval’s theorem,
we have:

$\displaystyle\boldsymbol{I}_{pqk}^{\mathrm{d}}=\mathbf{a}_{pqk}^{\prime}\mathbf{c}^{\prime\dagger}_{pqk},$
(40)

where $\mathbf{a}_{pqk}^{\prime}$ and $\mathbf{c}^{\prime}_{pqk}$ are vectors.
The singular value, $\eta_{pqk}$ indicates how best
$\eta_{pqk}\boldsymbol{I}_{pqk}^{\mathrm{d}}$ contributes to the quality and
fidelity of $\widetilde{\boldsymbol{I}}^{\mathrm{d}}_{\mathrm{pq}}$. A closer
look at Eq. 39 shows that the singular values in both the visibility and image
domains are equal, thanks to the linear and unitary properties of the Fourier
transform. This means that the choice of the domain where the data should be
compressed does not matter; the number of components to be retained in the
visibility domain would eventually be the same in the image domain.

Data: $\boldsymbol{V};\epsilon~{}\text{or}~{}\varepsilon$

Result: the list of the number of singular values to retain at each baseline
$\big{\\{}n_{pq}\big{\\}},\forall pq$

1 $n_{pq}\leftarrow 0;\boldsymbol{V}_{pq,n_{pq}}\leftarrow\boldsymbol{0}$

2 if  _$\epsilon$ is given_ then

3 for  _all baseline $pq$_ do

4

5 do

6 $\boldsymbol{V}_{pq,n_{pq}}\leftarrow\boldsymbol{V}_{pq,n_{pq}-1}$

7
$~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}+\eta_{pqn_{pq}}\boldsymbol{a}_{pqn_{pq}}\boldsymbol{c}_{pqn_{pq}}^{\dagger}$

8 while
_$\frac{||\boldsymbol{V}_{pq}-\boldsymbol{V}_{pq,n_{pq}}||_{\mathrm{F}}}{||\boldsymbol{V}_{pq}||_{\mathrm{F}}}
>\epsilon$ and $n_{pq}\leq r$_;

9 $n_{pq}\leftarrow n_{pq}+1$

10

11 end for

12

13else

14 for  _all baseline $pq$_ do

15 do

16 $\boldsymbol{V}_{pq,n_{pq}}\leftarrow\boldsymbol{V}_{pq,n_{pq}-1}$

17
$~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}+\eta_{pqn_{pq}}\boldsymbol{a}_{pqn_{pq}}\boldsymbol{c}_{pqn_{pq}}^{\dagger}$

18 while
_$\frac{||\boldsymbol{V}_{pq,n_{pq}}||_{\mathrm{F}}}{||\boldsymbol{V}_{pq}||_{\mathrm{F}}}\times
100 <\varepsilon$ and $n_{pq}\leq r$_;

19 $n_{pq}\leftarrow n_{pq}+1$

20

21 end for

22

23 end if

Algorithm 2 Finding $\big{\\{}n_{pq}\big{\\}},\forall pq$ using BDSVD. Here,
all the $\boldsymbol{V}_{pq}$ are taken from $\boldsymbol{V}$.

Alternatively, the dirty image is derived from summing Eq. 39 across all the
baselines:

$\displaystyle\widetilde{\boldsymbol{I}}^{\mathrm{d}}$
$\displaystyle=\sum_{pq}\widetilde{\boldsymbol{I}}^{\mathrm{d}}_{\mathrm{pq}}$
(41)
$\displaystyle=\sum_{pq}\bigg{(}\sum_{k=1}^{n_{pq}}\eta_{pqk}\boldsymbol{I}_{pqk}^{\mathrm{d}}\bigg{)}.$
(42)

If
$\boldsymbol{I}_{\mathrm{}}^{\mathrm{d}}=\sum_{pq}\mathcal{F}\boldsymbol{V}_{pq}$
is the Fourier transform of the uncompressed data (i.e. the uncompressed image
of the sky), then to quantify the net loss in signal per-pixel,
$\boldsymbol{I}_{\mathrm{loss}}^{\mathrm{d}}$, the following difference is
adopted as a standard fidelity metric:

$\displaystyle\boldsymbol{I}_{\mathrm{loss}}^{\mathrm{d}}$
$\displaystyle=\boldsymbol{I}^{\mathrm{d}}-\widetilde{\boldsymbol{I}}^{\mathrm{d}}$
(43)
$\displaystyle=\sum_{pq}\bigg{(}\sum_{k=1}^{r}\eta_{pqk}\boldsymbol{I}_{pqk}^{\mathrm{d}}-\sum_{k=1}^{n_{pq}}\eta_{pqk}\boldsymbol{I}_{pqk}^{\mathrm{d}}\bigg{)}.$
(44)

For $n_{pq}<r$, we have

$\displaystyle\boldsymbol{I}_{\mathrm{loss}}^{\mathrm{d}}$
$\displaystyle=\sum_{pq}\sum_{k=n_{pq}+1}^{r}\eta_{pqk}\boldsymbol{I}_{pqk}^{\mathrm{d}},$
(45)

where the entries of $\boldsymbol{I}_{\mathrm{loss}}^{\mathrm{d}}$ are
different from $0$.

Data: $\boldsymbol{V}$; $N_{p}$

Result: The SVD of $\boldsymbol{V}$ decomposed per baseline:
$\big{\\{}\mathbf{A}_{pq},\mathbf{\Lambda}_{pq},\mathbf{C}_{pq}\big{\\}},~{}\forall
pq$

1 for  _all baseline $pq$, in parallel_ do

2 for  _all $i$, in parallel_ do

3
$\boldsymbol{V}_{pq}^{\dagger}\boldsymbol{V}_{pq}+=\boldsymbol{V}_{pq}^{\mathrm{[i]\dagger}}\boldsymbol{V}_{pq}^{\mathrm{[i]}}$

4 end for

5

6
$\mathbf{C}_{pq}\mathbf{\Lambda}_{pq}^{2}\mathbf{C}_{pq}^{\dagger}=\boldsymbol{V}_{pq}^{\dagger}\boldsymbol{V}_{pq}$

7 $\mathbf{\Lambda}_{pq}=\sqrt{\mathbf{\Lambda}_{pq}^{2}}$

8 for  _all $i$, in parallel_ do

9
$\mathbf{A}_{pq}^{\mathrm{[i]}}=\boldsymbol{V}_{pq}^{\mathrm{[i]}}\mathbf{C}_{pq}\mathbf{\Lambda}_{pq}^{-1}$

10 end for

11
$\mathbf{A}_{pq}=\big{\\{}\mathbf{A}_{pq}^{\mathrm{[i]}}\big{\\}},~{}\forall
i$

12 end for

Algorithm 3 Finding
$\big{\\{}\mathbf{A}_{pq},\mathbf{\Lambda}_{pq},\mathbf{C}_{pq}\big{\\}},\forall
pq$ using parallel computing. Here, all the $\boldsymbol{V}_{pq}$ are taken
from $\boldsymbol{V}$, and $N_{p}$ is the number of compute nodes.

###  4.2 Noise filtering and S/N

Three important questions arise surrounding noise filtering when simple SVD
and BDSVD are used i) how does the noise behave for simple SVD where the same
number of components are retained on all the baselines? ii) how does BDSVD and
its varying compression factor affect the noise at each baseline and how does
the filtered noise differ from the simple SVD and traditional averaging? and
iii) on what type of baseline is noise heavily filtered when using simple SVD
and BDSVD? In this section, we address these questions through well-posed
conditioning and discussion.

The entries of the compact visibility data matrix $\boldsymbol{V}_{pq}$ come
from sampling the sum of $\mathcal{V}$ and noise $\mathcal{E}$ as shown in Eq.
8. Thus, each $\eta_{pqk}$ reflects the strength of the signal contribution
from a component which is sampled from $\mathcal{V}+\mathcal{E}$. Small
$\eta_{pqk}$ corresponds to components of $\boldsymbol{V}_{pq}$ that are
heavily corrupted by noise, which means that if one retains components with
the bigger $\eta_{pqk}$ this should be equivalent to removing noise and
keeping useful signal. This means that with simple SVD and BDSVD, the noise in
the compressed data is reduced for $n<r$ and $n_{pq}<r$, respectively. Below,
we provide details on the noise filtering capability of simple SVD and BDSVD
compared to traditional averaging at the same $CF$.

For the same compression factor, the analytical visibility noise penalty
estimates of simple SVD and BDSVD is the relative decrease in noise over
traditional averaging:

$\displaystyle\Xi_{X}$ $\displaystyle=\frac{\Sigma_{X}}{\Sigma_{avg}},$ (46)

where $\Sigma_{X}$ and $\Sigma_{avg}$ are compressed noise using simple SVD
(or BDSVD) and traditional averaging, respectively. The goal is to show that
$\Xi_{X}<0$ which means $\Sigma_{X}<\Sigma_{avg}$.

Assuming that for all baselines and samples, the noise term has constant
r.m.s. $\Sigma$, when denoising a signal by averaging $n_{avg}$ samples, the
reduction in noise is well understood to be
$\Sigma_{avg}=\Sigma/\sqrt{n_{avg}}$ if the noise is not correlated between
samples; which confirms that $\Sigma_{avg}<\Sigma$. As discussed, the SVD
removed some noise in the data, therefore, $\Sigma_{X}<\Sigma$. However, at
the same compression factor, it is not trivial to see that
$\Sigma_{X}<\Sigma_{avg}$ analytically (we refer the reader to Appendix A for
a detailed mathematical discussion). Empirical measurements are used in
Section 5 to show that at the same compression factor
$\Sigma_{X}<\Sigma_{avg}$ and therefore $\Xi_{X}<0$.

In the image domain, the noise penalty estimate of the centre pixel is given
by

$\displaystyle\Xi_{X}^{W}$
$\displaystyle=\frac{\sum_{pqt_{i}\nu_{j}}W_{pqt_{i}\nu_{j}}^{2}\Xi_{X}^{2}}{(\sum_{pqt_{i}\nu_{j}}W_{pqt_{i}\nu_{j}})^{2}},$
(47)

where $W_{pqt_{i}\nu_{j}}$ is the imaging weight per visibility.

As discussed above, components with smaller $\eta_{pqk}$ are strongly noisy
and are candidates for removal when noise is of concern. In addition,
components with larger $\eta_{pqk}$ contain the most signal to retain. Simple
SVD and BDSVD only retain components with larger $\eta_{pqk}$ for given
thresholding and components that do not meet this threshold are rejected.
BDSVD differs from simple SVD in that the number of components, $n_{pq}$ to
retain is baseline-dependent. As mentioned in Section 3.2, the first
$\eta_{pqk}$ are large for shorter baselines and decrease faster compared to
longer baselines, where $\eta_{pqk}$ decreases slowly (Figure 1 clearly shows
this behaviour of $\eta_{pqk}$). This means that if the same number of
components is retained on all the baselines as in simple SVD, then some of the
components that contain the signal will be discarded on longer baselines,
hence distorting the compressed data. On the other hand, on the shorter
baselines, several noisy components are retained in addition to all the
components with a strong signal. BDSVD takes advantage of this drawback and
constructs a baseline-dependent number of components to retain such that only
components with strong signal strength are retained on each baseline. This
heavily filters out noise and maintains signal fidelity compared to simple SVD
at the same $CF$. With this in mind, then at the same compression factor, we
can write:

$\displaystyle\Sigma_{\mathrm{svd}}$
$\displaystyle>\Sigma_{\mathrm{bdsvd}}\implies\Xi_{\mathrm{svd}}^{W}>\Xi_{\mathrm{bdsvd}}^{W},$
(48)

where $\Sigma_{\mathrm{svd}}$ and $\Sigma_{\mathrm{bdsvd}}$ (respectively
$\Xi_{\mathrm{svd}}^{W}$ and $\Xi_{\mathrm{bdsvd}}^{W}$) are the noise
(respectively the noise penalty) using simple SVD and BDSVD, respectively.
Analytically, Eq. 48 clearly shows that BDSVD reduces noise compared to simple
SVD. Empirical measurements are used in Section 5 to confirm the analytical
result in Eq. 48.

The metric we used to measure the signal to noise, $S/R$ in decibel in each
pixel of the compressed image is:

$\displaystyle S/R$
$\displaystyle=10\log_{10}\frac{\widetilde{\mathcal{I}}_{\boldsymbol{l}}}{\sigma_{\mathrm{pix}}+c_{\mathrm{noise}}},$
(49)

where $\widetilde{\mathcal{I}}_{\boldsymbol{l}}$ is the compressed version of
$\mathcal{I}_{\boldsymbol{l}}$; the sky without any electronics corruption and
effects that can disrupt the signal in the path towards the instrument (see
Eq. 2). In this formulation, $c_{\mathrm{noise}}$ represents the signal coming
from sources that are outside the FoV and $\sigma_{\mathrm{pix}}$ is the per-
pixel noise in the dirty image:

$\displaystyle\sigma_{\mathrm{pix}}^{2}$
$\displaystyle=\frac{\sum_{pqt_{i}\nu_{j}}W_{pqt_{i}\nu_{j}}^{2}\Sigma_{X}^{2}}{(\sum_{pqt_{i}\nu_{j}}W_{pqt_{i}\nu_{j}})^{2}}.$
(50)

##  5 Simulations

The MeerKAT and the The European Very Long Baseline Interferometry Network
(EVN) are used as reference telescopes in this section to evaluate and compare
the performance of each method; traditional averaging, BDA, simple SVD and
BDSVD. Three different metrics are evaluated; i) source amplitude loss is
measured vs. baseline lengths, baseline-dependent compression factor vs.
baseline lengths and east-west baseline lengths, ii) amplitude loss is
measured relative to the source’s position in the sky; iii) to measure the
spectral and temporal resolution, the Point Spread Function (PSF) shape is
measured relative to the source’s position in the sky and (iv) the S/R is also
measured relative to the source’s position in the sky.

Note that for a rank$-n$ SVD decomposition of two-dimensional visibility data,
$\boldsymbol{V}_{pq,n}$ as defined in Eq. 11; $M$ and $N$ are the number of
time steps and channels, respectively. In the space where
$\boldsymbol{V}_{pq,n}$ is defined, the row and column data are the time and
frequency observations, respectively. When $\boldsymbol{V}_{pq,n}$ is
projected into the rank$-n$ SVD space, the time and frequency observations are
mixed, which makes it difficult to completely determine the direction in which
the compression of rank$-n$ is performed.

To compare the proposed methods with traditional averaging and BDA, we need to
evaluate the above metrics in a precise compression direction, such as in time
or frequency. We adopt the following strategy, for example, to evaluate the
metrics in time: $\boldsymbol{V}_{pq,n}$ is scaled from a data matrix of size
$M\times N$ to size $M_{1}\times M_{2}\times N$, where rank$-n$ SVD is only
applied in the time direction with chunk of size $M=M_{1}\times M_{2}$. The
notation $CF=CF_{t}\times CF_{\nu}$ adopted in this section means that the
data is compressed by a factor of $CF_{t}$ in time and $CF_{\nu}$ in
frequency. The discussion in this section does not include the compression in
the frequency direction ($CF_{\nu}=1$); however, equivalent performance is
observed in the frequency direction if $\boldsymbol{V}_{pq,n}$ is scaled to
$M\times N_{1}\times N_{2}$ and the rank$-n$ SVD is applied only to entries
$N=N_{1}\times N_{2}$ of $\boldsymbol{V}_{pq,n}$.

###  5.1 MeerKAT telescope

For the MeerKAT telescope at 1.4 GHz, a point source is simulated during $166$
min $40$ s which gives $10~{}000$ time steps of 1 s each. The simulation has a
total bandwidth of 0.8 MHz divided into 10 channels each 80 kHz wide. The 1.4
GHz MeerKAT telescope observes a full FoV with radius $\sim 2.25$ deg; the
second null of the PB without time and frequency smearing effects with a time
step of 1 s and a channel width of 80 kHz. This high-resolution dataset is
scaled to $100\times 100\times 10$, where simple SVD and BDSVD are applied
separately at each frequency chunk of size $100\times 100$. The compression
factors we adopt are $CF=25\times 1$, $CF=35\times 1$ and $CF=50\times 1$,
which translates to a space-saving of $96.0\%$ , $97.14\%$ and $98.0\%$ in the
time direction, respectively. With this high-resolution dataset, $CF=50\times
1$ is the maximum time compression factor that can be achieved with the
proposed methods. This is because $CF_{t}=50$ for rank-$1$ SVD; in other words
only one component is kept across all baselines and the compression factor as
described in Eq. 33 becomes $100\times 100/(100+100+0.5)\sim 50$.

The strategies adopted for traditional averaging and BDA are:

  * •

Two low-resolution MSs with bins size $25s\times 5kHz$ and $50s\times 5kHz$
are created to receive the resampled visibilities for traditional averaging
with $CF=25\times 1$ and $CF=50\times 1$, respectively.

  * •

Using the method described in Section 2.2, a third MS is created to receive
the resampled visibility for BDA with $CF=25\times 1$ and $CF=50\times 1$.

  * •

A fourth MS of the same copy as the high-resolution MS is created to receive
resampled visibility for simple SVD and BDSVD. No backup policy is required to
save the decompressed visibility data for simple SVD and BDSVD. The adopted
compression factors are $CF=25\times 1$ and $CF=50\times 1$ for simple SVD,
while for BDSVD $CF=25\times 1$ and $CF=35\times 1$ are adopted. Since
$CF=50\times 1$ is the maximum compression factor BDSVD with $CF=50\times 1$
will be equivalent to simple SVD with $CF=50\times 1$. Due to the latter, we
do not run BDSVD with $CF=50\times 1$.

####  5.1.1 Amplitude vs. east-west baseline lengths

This first test aims to quantify the decorrelation of a single point source of
1 Jy amplitude placed at 2.25 deg; the second null of the MeerKAT at 1.4 GHz
PB. The results are shown in Figure 3.

For $\varepsilon=99\%$, Figure 3 (top panel) shows the baseline-dependent
compression factor $CF_{pq}$ in logarithmic scale against increasing east-west
baseline lengths. The rank$-n$ SVD-related methods, $CF_{pq}$ is computed from
Eq.33 after calculating all $n_{pq}$ using Algorithm 2, while for BDA
$CF_{pq}$ is computed following the discussion in Atemkeng et al. (2018). This
is an important result of this study, so it deserves a detailed explanation.

Although we understand the $CF_{pq}$ of BDA as discussed in Atemkeng et al.
(2018), the $CF_{pq}$ of BDSVD shows a different pattern and strength for
east-west projected baselines. It can be observed that for BDA the $CF_{pq}$
are strictly decreasing for increasing east-west projected baseline length,
whereas for BDSVD, east-west projections of equal lengths can have different
$CF_{pq}$ where some $CF_{pq}$ are significantly larger than others. However,
this behaviour still shows that more data is compressed on the smaller east-
west baseline lengths compared to the longer; although the pattern and
strength are not similar to that of BDA. Figure 3 (middle panel) shows the
$CF_{pq}$ as a function of increasing baseline length. We observe that
baselines of equal lengths can have different $CF_{pq}$ when BDA or BDSVD are
in effect. Whereas $CF_{pq}$ shows more decreasing strengths with BDSVD and
randomly distributed with BDA, which confirms that decorrelation is not a
degree of baseline length but rather the east-west baseline length. At the
same overall compression factor $CF$, BDA and BDSVD see strictly different
$CF_{pq}$.

Figure 3 (bottom panel) shows the source amplitude of 1 Jy as a function of
increasing east-west projection lengths. It is clear from this result that
BDSVD outperformed BDA. The source amplitude is attenuated with BDA compared
to BDSVD which retains more than $\varepsilon=99\%$ of the source amplitude.
Additionally, we observe with BDSVD that the amplitude of the decorrelated
source is not constant for some equal east-west baseline lengths. Indeed, as
observed in Figure 3 (top panel); $CF_{pq}$ varies for certain equal east-west
baseline lengths, which can result in different degrees of decorrelation.
Similar behaviour can be observed with BDA, where decorrelation is revealed
faster for east-west baseline lengths belonging to the same averaging range.
For example, all east-west baselines with a length between $2.2~{}km$ and
$3.5~{}km$ fall within the same averaging range where $\sim 2$ bins are
averaged together. Since the length of the east-west baseline in the same
averaging range is different; explaining the different degrees of
decorrelation.

Figure 3: The top panel shows the baseline-dependent compression factor
$CF_{pq}$ in the logarithmic scale against increasing east-west baseline
lengths while the middle panel shows $CF_{pq}$ in logarithmic scale as a
function of baseline length. The bottom panel shows the amplitude of the 1 Jy
point source against increasing east-west baseline lengths. BDA sees strictly
different $CF_{pq}$ than BDSVD while BDSVD outperformed BDA in preserving the
1 Jy source amplitude.

####  5.1.2 Amplitude vs. source position

This section examines the result of the amplitude-decorrelation with respect
to the source position in the sky. The noise penalty is also measured and
compared for each of the compression methods. To measure the thermal noise,
the high-resolution dataset as described above is used to populate an empty
sky with 1 Jy thermal noise, where the different compression methods are
applied.

Figure 4: MeerKAT telescope observing at 1.4 GHz during $166$ min $40$ s with
a total bandwidth of 0.8 MHz; demonstrating the degree of a 1 Jy source
amplitude loss at various sky positions and the associated noise penalty.
Smearing against source distance from the phase centre, for traditional
averaging and BDA with $CF=25\times 1$ and $CF=50\times 1$ and for simple SVD,
BDSVD with $CF=25\times 1$, $CF=35\times 1$ and $CF=50\times 1$. The space-
saving $SS$ and the noise penalty $\Xi$ are given relative to the traditional
averaging bins.

Figure 4 shows the results of time compression. There are important points to
note about these results. A compression factor of $CF=25\times 1$ is needed to
retain at least $\varepsilon=99\%$ of the source amplitude. At this
compression regime, simple SVD and BDSVD keep the source amplitude almost flat
with negligible attenuation of $1\%$ that starts around a radius of $\sim 4$
deg upwards, while the $1\%$ attenuation starts around a radius of $\sim 1.32$
deg and $\sim 0.8$ deg upward for the BDA and traditional averaging,
respectively. A compression factor of $CF=50\times 1$ provides up to a FoV
with radius $\sim 2.9$ deg for $\varepsilon=95\%$ whereas at this compression
factor, traditional averaging and BDA can only provide a FoV with radius up to
$\sim 1.32$ deg and $\sim 2$ deg, respectively.

It can also be noted that traditional averaging only provides a FoV with a
radius of $\sim 1.32$ deg for $CF=25\times 1$ whereas to obtain the same FoV
with a radius of $\sim 1.32$ deg using simple SVD, we can compress the data by
two-order of magnitude higher than traditional averaging. Traditional
averaging or BDA cannot compress data by a factor of $CF=25\times 1$ and
achieve $\sim 2.25$ deg radius with $1\%$ attenuation. However, at this $1\%$
attenuation rate, BDSVD achieves $\sim 2.25$ deg radius for $CF=35\times 1$.
BDSVD draws its potential from the advantages of both BDA and simple SVD.

The values of $\Xi\sim 1$ for BDA, this result has also been observed in
Atemkeng et al. (2018) and at the same compression factor $CF$, the reduction
in noise of both BDA and traditional averaging varies with $\sqrt{CF}$. Values
of $\Xi<1$ for simple SVD, which demonstrates the common use of SVD to denoise
signals. The values of $\Xi$ for BDSVD are lower than that of simple SVD which
confirms the theoretical result discussed in Section 4.2; the noise
performance of BDSVD exceeds that of simple SVD at the same compression
factor.

####  5.1.3 PSF distortion vs. source position

Figure 5: MeerKAT telescope observing at 1.4 GHz during $166$ min $40$ s with
a total bandwidth of 0.8 MHz. The PSF of a source at $2.25$ deg is shown when
the high-resolution dataset is compressed for $CF=25\times 1$ with BDSVD (top
left), simple SVD (top right), BDA (middle left), traditional averaging
(middle right), the vertical (bottom left) and horizontal (bottom right)
directions of the normalised PSFs. Figure 6: MeerKAT telescope observing at
1.4 GHz during $166$ min $40$ s with a total bandwidth of 0.8 MHz. The average
of the radial and tangential PSF resolution measured at the FWHM are shown at
various sky positions when traditional averaging and BDA with $CF=25\times 1$
and $CF=50\times 1$ and for simple SVD, BDSVD with $CF=25\times 1$,
$CF=35\times 1$ and $CF=50\times 1$ are applied. The PSF is not distorted at
the edges of the FoV when the SVD-related methods applied. Figure 7: MeerKAT
telescope observing at 1.4 GHz during $166$ min $40$ s with a total bandwidth
of 0.8 MHz. The S/N as defined in Eq 49 is measured at various sky positions
for compression factors $CF=25\times 1$, $CF=35\times 1$ and $CF=50\times 1$
with traditional averaging, BDA, simple SVD and BDSVD. The SVD-related methods
offer superior S/N.

As discussed in Section 1, a negative aspect of averaging the visibility data
is the distortion of the PSF for which the longer baselines are the major
contributors. The amplitude of the PSF at a given position in the sky provides
a measure of the signal loss while its width (say at the FWHM) describes how
widely the source is spread out at that position in the sky. In this section,
we evaluate the PSF distortion by measuring the width of the PSF at the FWHM
when each of the discussed compression methods are applied in time. The
simulated high-resolution dataset of $1s\times 80kHz$ bins discussed in
Section 5.1 is reused where the three compression factors $CF=25\times 1$,
$CF=35\times 1$ and $CF=50\times 1$ are compared for all the different
methods.

Results are shown in Figures 5 and 6. Figure 5 depicts the PSF of a source at
$2.25$ deg when the high-resolution dataset is compressed for $CF=25\times 1$
with BDSVD (top left), simple SVD (top right), BDA (middle left) and
traditional averaging (middle right). Figure 5 displays the vertical (bottom
left) and horizontal (bottom right) directions of the normalised PSFs as well.
It is clear from this visual inspection that the width of the PSF is different
for all the compression methods, with traditional averaging having a wider PSF
compared to the other methods. Although BDA is a potential compression method
that retains amplitude loss, the resulting PSF shape differs completely from
that of traditional averaging and the SVD-related methods.

In order to quantify the width of the PSFs, Figure 6 displays the average of
the radial and tangential PSF resolutions measured at the FWHM against
distance from the phase centre of the observation. BDSVD shows excellent
capabilities that maintain the PSF without any distortion up to a radius
greater than $3$ deg with $CF=25\times 1$ and $CF=35\times 1$. A similar
result is observed for simple SVD with $CF=25\times 1$. However, for simple
SVD the PSF distortion starts at $\sim 1.2$ deg radius for $CF=50\times 1$.
When traditional averaging or BDA is in action, the PSF quickly begins to
observe distortion starting at around $\sim 1.2$ deg for each of the
compression factors. These results confirm that BDSVD and simple SVD can
maintain the PSF without any distortion up to above a radius of $2.25$ deg;
the maximum radius that the MeerKAT telescope at 1.4 GHz is capable to
achieve.

####  5.1.4 Relative S/N

The S/N as discussed theoretically in Section 4.2 is shown in Figure 7 as a
function of source position in the sky. The simulation in Section 5.1.2 is
used to measure the source amplitude when the compression methods are applied
with the different compression factors.

Two simulations are included in this section to measure the thermal noise and
noise from sources outside the FoV. The high-resolution dataset as described
above is used to populate i) an empty sky with 1 Jy thermal noise and ii) a 10
Jy source is simulated at 20 deg away from the phase centre, and the different
compression methods are applied separately to each simulation. The far-field
contamination is measured from images of size $2^{10}\times 2^{10}$. Note that
apart from the signal from the far-field source, this image should be empty.
BDSVD and simple SVD benefit from a strong improvement of S/N to about $\sim
1.5$ dB at $2.5$ deg compared to traditional averaging. However, BDA does not
improve the S/N compared to traditional averaging; this is understood since
thermal noise reduction scales are the same for BDA and traditional averaging,
while BDA does not suppress sources that are outside the FoV compared to
traditional averaging.

Figure 8: EVN telescope observing at 1.6 GHz during $15$ min with a total
bandwidth of $50$ kHz; demonstrating the degree of smearing for a 1 Jy source
amplitude loss at various sky positions and the associated noise penalty.
Smearing against source distance from the phase centre, for traditional
averaging, BDA, simple SVD and BDSVD with $CF=10\times 1$, $CF=25\times 1$,
$CF=50\times 1$, $CF=95\times 1$ and $CF=150\times 1$. The noise penalty $\Xi$
is given relative to the traditional averaging bins.

###  5.2 EVN

We investigate the application of simple SVD and BDSVD in VLBI to keep
decorrelation down to a certain level while significantly compressing the
data. We simulate a 15 min full EVN (i.e. Badary, Effelsberg, Hartebeesthoek,
Jodrell Bank, Medicina, Noto, Onsala, Shanghai, Svetloe, Torun, Westerbork,
Zelenchukskaya) at 1.6 GHz observation. With a total bandwidth of $50$ kHz
channelised into 10 channels each of width $5$ kHz, the 15 min observation is
sampled at each $0.01$ s which gives $90~{}000$ time steps. To apply the
simple SVD and BDSVD, the $90~{}000$ time steps are scaled to $300\times 300$,
which allows us to investigate compression factors of $CF=10\times 1$,
$CF=95\times 1$ and $CF=150\times 1$. A single point source with 1 Jy
amplitude is simulated and the amplitude of the source is measured and
compared to that of traditional averaging and BDA. The results are shown in
Figure 8. It is observed that for a smearing factor of $1\%$ and $CF=10\times
1$, traditional averaging and BDA would result in a FoV of 6 arcmin and 18
arcmin, respectively. Whereas with $CF=95\times 1$, BDSVD can image a FoV of
above 30 arcmin with a superior noise reduction capability. Although these
results give a good taste in the VLBI regime, the implementation of the
different compression methods on real data is imperative.

##  6 Conclusions

One of the major contributors to the large data volumes is the long baselines
of an array as they influence the degree of data sampling required to avoid
decorrelation of the astrophysical signal. With BDA, the sampling rate is
baseline-dependent as short baselines can be sampled far more coarsely leading
to a high compression rate. BDA is an established technique for compressing
radio interferometric data. However, this technique results in irregularly
sampled data which, while technically supported by the MS format via variable
frequency bins supported by many spectral windows, requires the data to be
restructured in ways that reduce performance when processing. This work shows
an approach that uses a low-rank matrix approximation to achieve greater
compression rates while significantly minimising smearing compared to BDA or
traditional averaging. What is even more exciting is that when combined the
low-rank approximation and baseline-dependent formalism, BDSVD, effectively
eliminates smearing within the FoV for the MeerKAT telescope and the EVN.
However, this method has three caveats.

Firstly, SVD is computationally expensive and may be impractical for large
datasets when performed on every baseline independently. Fortunately, we show
that the computations can be distributed across multiple computer nodes as
discussed in Algorithms 3.

Secondly, simple SVD and BDSVD do not affect current calibration algorithms or
software. The possibility of calibrating visibility data remains since the
compressed data is decompressed before it is calibrated. This implies that the
temporal and frequency resolutions are recovered and consequently, do not
modify the solution intervals required for calibration. We also note that the
capability of current imaging software is maintained. A very important result
to note is that linear transformation between the compressed data, the
compressed noise and the source is maintained, thus showing that it is
possible to image the sky in the compressed domain if
$\boldsymbol{I}_{pqk}^{\mathrm{d}}$ is well understood. This opens a potential
research investigation that could possibly see the data be imaged in the SVD
space.

Thirdly, the implementation discussed in Section 2.2 describes a method which
compresses BDA data to an MS-compatible format on both disk and in memory. By
contrast the SVD compression schemes described here produce singular values
that must be expanded "on-the-fly" to full resolution in memory. There is,
therefore, an interesting contrast between the two approaches as the first
decreases both the amount of data and computation (FLOPs) performed by
algorithms on that data. This is achieved at the cost of implementation
complexity, especially with regards to calibration techniques which must
reason about solution intervals that lie across multiple data points on the
sparse domain in which BDA data lies. Simple SVD and BDSVD offer superior
compression to BDA and traditional averaging, but requires processing data at
full data resolution, for reduced implementation complexity. Decompressing
data "on-the-fly" is already offered in the PyDATA ecosystem by packages such
as BLOSC (Haenel, 2014), which compresses data into chunks and decompresses
them exactly into the L1 cache of a CPU. This reduces the memory intensity of
an algorithm, which is beneficial since modern CPUs are starved for data by
their memory system (Alted, 2010) and radio interferometric data is very
large. In such a regime, it may be possible to transmit the singular values
into the L1 cache of a CPU and expand the data to full resolution for
consumption by algorithms. Additional possible future work would investigate
an online compression-decompression strategy using the full archiving
potential of SVD-related methods.

Real data that is contaminated with Radio Frequency Interference (RFI), power
is distributed across all visibilities when decompressing the SVD. There are
two ways to avoid this situation: either the RFIs are removed from the data
before applying the SVD or the SVD is used to remove the RFI during
compression. The latter requires intense investigation as future work, since
removing the RFI is equivalent to removing the higher singular value
components. Removing higher singular value components is problematic because
it is unclear where weak radio signals will begin to be suppressed in the
process.

Real data usually has a few entries flagged for some reason. Interesting
future work would be to investigate how these flagged entries affect singular
values. We expect that the situation will fall somewhere between the two below
scenarios:

  * •

If the flagged real data entries belong to the same neighbourhood and are
assigned the same value; this increases the similarity in the neighbourhood
and hence SVD will result in a lower rank (consequently, the compression
factor will be bigger) compared to no flagging; and

  * •

If the flagged real data entries are not from the same neighbourhood (e.g.
randomly flagged data), the rank of the SVD will be higher (consequently, the
compression factor will be smaller) compared to no flagging.

To address this problem, a potential solution would be to perform in-painting
on the model data where flags are present and then use this in-painted data to
fill in the flagged entries in the uncompressed data prior to applying the
SVD.

## Acknowledgements

This work is based upon research supported by the South African Research
Chairs Initiative of the Department of Science and Technology and National
Research Foundation. The MeerKAT telescope is operated 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. MA
acknowledges support from Rhodes University.

## Data Availability

The data underlying this article will be shared on reasonable request to the
corresponding author

## References

  * Alted (2010) Alted F., 2010, Computing in Science & Engineering, 12, 68
  * Ansah-Narh et al. (2020) Ansah-Narh T., Akoto-Danso A., Proven-Adzri E., Takyi B., 2020, Journal of Instrumentation, 15, P04020
  * Atemkeng et al. (2016) Atemkeng M., Smirnov O., Tasse C., Foster G., Jonas J., 2016, Monthly Notices of the Royal Astronomical Society, 462, 2542
  * Atemkeng et al. (2018) Atemkeng M., Smirnov O., Tasse C., Foster G., Keimpema A., Paragi Z., Jonas J., 2018, Monthly Notices of the Royal Astronomical Society, 477, 4511
  * Atemkeng et al. (2020) Atemkeng M., Smirnov O., Tasse C., Foster G., Makhathini S., 2020, Monthly Notices of the Royal Astronomical Society, 499, 292
  * Bobin et al. (2008) Bobin J., Starck J.-L., Ottensamer R., 2008, IEEE Journal of Selected Topics in Signal Processing, 2, 718
  * Cotton (2009) Cotton W., 2009, OBIT Development Memo Series
  * Demmel et al. (2007) Demmel J., Dumitriu I., Holtz O., 2007, Numerische Mathematik, 108, 59
  * Dewdney et al. (2009) Dewdney P. E., Hall P. J., Schilizzi R. T., Lazio T. J. L., 2009, Proceedings of the IEEE, 97, 1482
  * Gu & Eisenstat (1996) Gu M., Eisenstat S. C., 1996, SIAM Journal on Scientific Computing, 17, 848
  * Haenel (2014) Haenel V., 2014, CoRR, abs/1404.6383
  * Hamaker et al. (1996) Hamaker J., Bregman J., Sault R., 1996, Astronomy and Astrophysics Supplement Series, 117, 137
  * Iheanetu et al. (2019) Iheanetu K., Girard J., Smirnov O., Asad K., de Villiers M., Thorat K., Makhathini S., Perley R., 2019, Monthly Notices of the Royal Astronomical Society, 485, 4107
  * Johnston et al. (2008) Johnston S., et al., 2008, Experimental astronomy, 22, 151
  * Jonas (2009) Jonas J. L., 2009, Proceedings of the IEEE, 97, 1522
  * Kartik et al. (2017) Kartik V. S., Carrillo R. E., Thiran J.-P., Wiaux Y., 2017, Monthly Notices of the Royal Astronomical Society, 468, 2382
  * Kemball & Wieringa (2000) Kemball A., Wieringa M., 2000, URL: http://casa. nrao. edu/Memos/229. html
  * Lee et al. (2016) Lee J., Kim S., Lebanon G., Singer Y., Bengio S., 2016, Proceedings of Machine Learning Research
  * Offringa et al. (2010) Offringa A., De Bruyn A., Biehl M., Zaroubi S., Bernardi G., Pandey V., 2010, Monthly Notices of the Royal Astronomical Society, 405, 155
  * Shishkin et al. (2019) Shishkin S. L., Shalaginov A., Bopardikar S. D., 2019, Numerical Linear Algebra with Applications, 26, e2246
  * Smirnov (2011) Smirnov O. M., 2011, Astronomy & Astrophysics, 527, A106
  * Stewart (1998) Stewart G. W., 1998, Matrix algorithms: volume 1: basic decompositions. SIAM
  * Tasse et al. (2018) Tasse C., et al., 2018, Astronomy & Astrophysics, 611, A87
  * Wijnholds et al. (2018) Wijnholds S., Willis A., Salvini S., 2018, Monthly Notices of the Royal Astronomical Society, 476, 2029
  * Zarka et al. (2015) Zarka P., et al., 2015, in 2015 International Conference on Antenna Theory and Techniques (ICATT). pp 1–6
  * van Haarlem et al. (2013) van Haarlem M. P., et al., 2013, Astronomy & astrophysics, 556, A2

##  Appendix A Keeping track of the noise

Given an observation (assuming $M<N$)

$\underset{M\times N}{Y}=\bar{Y}+\epsilon=\underset{M\times
M}{U}\underset{M\times M}{\Lambda}\underset{M\times N}{V},$ (51)

with

$\mbox{vec}(\epsilon)\quad\sim\quad N(0,\Sigma).$ (52)

Can we keep track of the noise properties after thresholding the singular
values? We can use the decomposition to keep track of the covariance of the
singular values. To do this we note that, since $U$ and $V$ are unitary
matrices, we have

$\Lambda=U^{\dagger}(\bar{Y}+\epsilon)V^{\dagger}=\bar{\Lambda}+\epsilon_{\Lambda},$
(53)

with $\bar{\Lambda}=U^{\dagger}\bar{Y}V^{\dagger}$ and
$\epsilon_{\Lambda}=U^{\dagger}\epsilon V^{\dagger}$. Now the covariance in
SVD space is given by

$\displaystyle\Sigma_{\Lambda}$ $\displaystyle=$
$\displaystyle\mathbb{E}[\mbox{vec}(\epsilon_{\Lambda})\mbox{vec}(\epsilon_{\Lambda}^{\dagger})]$
(54) $\displaystyle=$ $\displaystyle\mathbb{E}[(V\otimes
U^{\dagger})\mbox{vec}(\epsilon)\mbox{vec}(\epsilon)^{\dagger}(V\otimes
U^{\dagger})^{\dagger}]$ (55) $\displaystyle=$ $\displaystyle(V\otimes
U^{\dagger})\mathbb{E}[\mbox{vec}(\epsilon)\mbox{vec}(\epsilon)^{\dagger}](V\otimes
U^{\dagger})^{\dagger}$ (56) $\displaystyle=$ $\displaystyle(V\otimes
U^{\dagger})\Sigma(V\otimes U^{\dagger})^{\dagger},$ (57)

where we have used the identity $\mbox{vec}(ABC)=(C^{\dagger}\otimes
A)\mbox{vec}(B)$ and the fact that

$\mathbb{E}[\mbox{vec}(\epsilon)\mbox{vec}(\epsilon)^{\dagger}]=\Sigma.$ (58)

Notice that, from the mixed product property of the Kronecker product, the
noise in SVD space will be i.i.d. whenever $\epsilon$ is i.i.d (i.e. $\Sigma$
is a multiple of the identity). Keeping track of the uncertainty after
thresholding is trivial when this is the case because the singular values are
independent. Things are a bit more complicated when there are off diagonal
entries in $\Sigma_{\Lambda}$. It is likely that the $\Sigma_{\Lambda}$ will
be diagonally dominant when $\Sigma$ is diagonal but this has to be checked.
If it is we can probably safely just treat it as a diagonal matrix.