{'fraction_non_alphanumeric': 0.0550354609929078, 'fraction_numerical': 0.01475177304964539, 'mean_word_length': 3.96164672765658, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Performance Analysis of Distributed Radio Interferometric Calibration

Distributed calibration based on consensus optimization is a computationally efficient method to calibrate large radio interferometers such as LOFAR and SKA. Calibrating along multiple directions in the sky and removing the bright foreground signal is a crucial step in many science cases in radio interferometry. The residual data contain weak signals of huge scientific interest and of particular concern is the effect of incomplete sky models used in calibration on the residual. In order to study this, we consider the mapping between the input uncalibrated data and the output residual data. We derive an analytical relationship between the input and output probability density functions which can be used to study the performance of calibration.

Abstract-Distributed calibration based on consensus optimization is a computationally efficient method to calibrate large radio interferometers such as LOFAR and SKA. Calibrating along multiple directions in the sky and removing the bright foreground signal is a crucial step in many science cases in radio interferometry. The residual data contain weak signals of huge scientific interest and of particular concern is the effect of incomplete sky models used in calibration on the residual. In order to study this, we consider the mapping between the input uncalibrated data and the output residual data. We derive an analytical relationship between the input and output probability density functions which can be used to study the performance of calibration.

I. INTRODUCTION
Most challenging science cases in modern radio astronomy are after weak signals that are hidden under noise and bright foregrounds (see, e.g., [1], [2]). The main goal of calibration is the correction for systematic errors in the data and the removal of contaminating foregrounds from this data to reveal such weak signals. Consensus optimization [3] has proved to be a computationally efficient solution for calibration [4]- [8] as well as for imaging [9]- [12] massive amounts of radio interferometric data. Calibration is always imperfect due to the errors in the input sky model as well as the consensus polynomials being used. Cramer-Rao lower bounds [13]- [16] have been used to study the asymptotic variance of estimation error of calibration parameters. Translating this bound to the error in the residual is however, cumbersome. In order to overcome this, calibration is considered as a nonlinear regression and Jacobian leverage [17]- [19] is proposed in [20], [21] to get limits on the variance of the residuals.
In this paper, we consider the mapping between the input uncalibrated data and the output residual data, where the residual is obtained after calibration and removal of the bright foreground signals. We derive an analytic relationship between the probability density functions (PDFs) of the input and output. In order to do that we use developments in bi-level optimization [22], [23] and matrix differentiation [24] to find derivatives of argmin function used in calibration. This work is supported by Netherlands eScience Center (project DIRAC, grant 27016G05).
The rest of the paper is organized is as follows. In section II, we give an overview of distributed calibration in radio interferometry. In section III, we derive analytic relationships for the performance of distributed calibration. Next, in section IV we derive a relationship between the input and output PDFs and finally, we draw our conclusions in section V.
Notation: Lower case bold letters refer to column vectors (e.g. y). Upper case bold letters refer to matrices (e.g. C). Unless otherwise stated, all parameters are complex numbers. The matrix inverse, transpose, Hermitian transpose, and conjugation are referred to as (.) −1 , (.) T , (.) H , (.) ⋆ , respectively. The matrix Kronecker product is given by ⊗. The vectorized representation of a matrix is given by vec(.). The identity matrix of size N × N is given by I N . Estimated parameters are denoted by a hat, (.). All logarithms are to the base e, unless stated otherwise. The Frobenius norm is given by . .

II. RADIO INTERFEROMETRIC CALIBRATION
In this section we give a brief overview of the data model used in radio interferometric calibration [25], [26]. We consider the radio frequency sky that is part of the sky model to be composed of discrete sources, far away from the earth such that the approaching radiation from each one of them appears to be plane waves. There are N receiving elements with dual polarized feeds in the array and at the p-th station, this plane wave causes an induced voltage, which is dependent on the beam attenuation as well as the radio frequency receiver chain attenuation. Consider the correlation of signals at the pth receiver and the q-th receiver, at frequency f , with proper signal delay. After correlation, the correlated signal of the p-th station and the q-th station (named as the visibilities), V pqf (∈ C 2×2 ) is given by In (1), J pf and J qf (∈ C 2×2 ) are the Jones matrices describing systematic errors at frequency f , at stations p and q, respectively. These matrices represent the effects of the propagation medium, the beam shape and the receiver. The noise matrix is given as N pqf (∈ C 2×2 ). The intrinsic signal on baseline pq is given by the coherency matrix C pqf (∈ C 2×2 ). For a linearly polarized source in the sky, with Stokes parameters I pqf , Q pqf , U pqf , V pqf , we have where φ pqf is the Fourier phase component that depends on the direction in the sky as well as the separation of stations p and q. For baseline coordinates u pq , v pq , w pq (in wavelengths) and direction cosines l, m, n, we have φ pqf = −2π(u pq l + v pq m + w pq (n − 1)). The noise matrix N pqf is assumed to have elements with zero mean, complex Gaussian entries with equal variance in real and imaginary parts but the statistics will vary because of the unmodelled structure [27]- [29]. The cost function that is minimized is given as where the systematic errors for all N stations are grouped as we can select the systematic errors for the station p as A p J f . Note that in (3), the summation is taken over all the baselines pq that have data, within a small bandwidth and time interval within which the systematic errors are assumed to be fixed. Consensus optimization problem is formulated as follows. First, we create the augmented Lagrangian as where the subscript (.) f denotes data (and parameters) at frequency f . In (6), g f (J f ) is the original cost function as in (3). The Lagrange multiplier is given by representing the basis functions evaluated at frequency f . The regularization parameter is given by ρ (∈ R + ). The alternating direction method of multipliers (ADMM) iterations n = 1, 2, . . . for solving (6) are given as where we use the superscript (.) n to denote the n-th iteration where (7) to (10) are executed in order. The steps (7), (9) and (10) are done for each f in parallel, at each compute (slave) node. The slave nodes are distributed across a network of computers. The update of the global variable in (8) is done in closed form at the fusion center. The extension of this data model to a multi-source scenario can be found in e.g., [8].
III. PERFORMANCE ANALYSIS At convergence, the closed form update of the global variable Z is where B fi corresponds to the consensus polynomial terms evaluated at frequency f i . We separate one frequency f i = f from the other P − 1 frequencies to get where and Note that R (∈ C 2F N ×2 ) in (14) has no dependence on the variables at frequency f , i.e., J f and Y f . Substituting (12) to (6), we get The gradients of (15) with respect to J f , Y f are given as and 2N ×2N ) and r 1 (R 1 ),r 2 (R 2 ) are the remainder terms that are independent of J f and Y f . The proof is given in appendix I. The gradient of the original cost function (3) is the derivation and be found in [4]. At a local minimum, we have Consider x p ′ q ′ r ∈ R to be one data point out of many that completes a full observation V pq , p, q ∈ [1, N ],p = q. This data point belongs to p = p ′ , q = q ′ correlation pair and r ∈ [1,8]. We select the value of r to represent one real or imaginary value of V p ′ q ′ (∈ C 2×2 ). Note that each complex number is considered as two data points. For instance, if r = 1, we represent the real part of V p ′ q ′ (1, 1). If r = 2, the imaginary part of V p ′ q ′ (1, 1) is selected, and so on. In order to find ∂J ∂x p ′ q ′ r , we take the derivative of both sides of (19) and (20) as in [22], [23] and we get The proof can be found in appendix II. Note that ∂V p ′ q ′ ∂x p ′ q ′ r in (21) will give a matrix (∈ C 2×2 ) with all zeros, except one real or imaginary value equal to 1, depending on the value of r.
The residual is calculated by subtracting the calibrated model from the data as We take the derivative of the residual with respect to x p ′ q ′ r and using [24], we get Note that vec (24) is zero except when p = p ′ , q = q ′ . Using (21) and (24), we can study the behavior of the residual with respect to small changes in input data. In section IV, we develop this further to consider the relationship between the input data and output residual PDFs.

IV. PROBABILITY DENSITY FUNCTIONS
We reformulate (3) as a vector optimization problem, for the sake of simplicity. The vectorized form of (1), v pq = vec(V pq ) can be written as where n pqf = vec(N pqf ). Depending on the time and frequency interval within which calibration solutions are obtained, we can stack up all cross correlations within that interval as x is a vector of size D × 1 of real data points. For a single time sample, D = 8N (N − 1)/2 because each (unique) cross correlation produces 8 real data points. One element out of this vector is x p ′ q ′ r (considered in section III), where p ′ , q ′ denote the pair of receivers forming the correlation and r is one data point out of the 8 produced by each correlation. We have the data model where θ is the real parameter vector (size M × 1) that is estimated by calibration. The parameters θ are the elements of J pf -s, with real and imaginary parts considered separately. The maximum likelihood (ML) estimate of θ under zero mean, white Gaussian noise is obtained by minimizing the least squares cost where The residual using calibration solution θ is obtained as The CRLB [13]- [16] is used to find a lower bound to the variance of θ. However, relating this lower bound to the residual y is not simple. Using Jacobian leverage [17]- [19], it is possible to obtain limits for the variance of y [20], [21] but we are after a simpler approach.
For m = [1, . . . , D], consider x m to be one element of x, and this is the same x p ′ q ′ r considered in section III, except we use m as the index instead of p ′ , q ′ , r. Taking the derivative of the residual with respect to x m , we have Using the chain rule (at θ = θ) where ∂s(θ) ∂θ T ∈ R D×M and ∂θ ∂xm ∈ R M×1 . At the solution, the gradient of the cost function is zero, i.e., ∂f (x,θ) ∂θ | θ= θ = 0. Following [22], [23], we have 1 where Consider the mapping from x to y, where T(·) is a composite of calibration and model subtraction to get the residual. Let the joint probability density functions of x and y be p X (x) and p Y (y), respectively. We can find p Y (y) by looking at the statistics of the residual, but scientific interest is in finding p X (x), so we use where J ∈ R D×D is the Jacobian of the mapping T(·), We use (31) to evaluate each column of J . Using (31), we can rewrite (38) as where θ) . . . f XD θ (x, θ)]| θ= θ (40) The determinant of J can be given using eigenvalues of A, i.e. λ j (A) as [30] |J | = exp The evaluation of (40) and moreover its eigenvalues is an expensive task. However, we can use the results of section III to simplify this. The closed form expressions (21) and (24) can be used to evaluate elements of (40) in closed form. The only requirement is the careful mapping of index m to indices p ′ , q ′ , r and vice versa. Moreover, the inversion of (f θθ (x, θ)) in (40) is not explicitly needed because we can use an iterative algorithm to find the eigenvalues of A such as by using implicitly restarted Arnoldi methods [31].
The extension of this work to study the performance of multi-source calibration is straightforward. We need to partition θ into partitions corresponding to each direction and we can evaluate (40) in block partitioned form. Furthermore, the results can also be used to study calibration without consensus optimization by setting ρ = 0 in ((21).

V. CONCLUSIONS
We have derived closed form relations for the performance analysis of distributed radio interferometric calibration. To study the weak signals buried in the data, preservation of their statistical behavior is essential. Using this work, we are able to study the effect of calibration in possible transformations of input data and if needed, compensating for these effects. We will produce software based on this work to accompany our calibration software as future work.
APPENDIX I: Proof of (16) and (17) First note that if g 1 = Y H (AJ + BY + C) and g 2 = AJ + BY + C 2 for some arbitrary constant matrices A, B and C, then the derivatives with respect to J and Y are, grad(g 1 , J) = 1 2 A H Y, grad(g 1 , Y) = 1 2 (AJ+(B+B H )Y+ C), grad(g 2 , J) = A H (AJ + BY + C) and grad(g 2 , Y) = B H (AJ + BY + C). Using this to find the gradient of (15) we get and substitution F = I − ρBP leads to (16) and (17). II: Proof of derivative (21) Taking the differential of (17) at the solution (first using grad(L f , Y) = 0) and substituting this to the differential of (16) at the solution (using grad(L f , J) = 0) We use the chain rule to expand dgrad(g f (J f ), J f ) as dvec (grad(g f (J f ), J f )) = D J grad(g f (J f ))vec (dJ)) (46) where D J grad(g f (J f )) is found by using definition 4 of [24]. In other words, if G(J, J ⋆ ) is a matrix function of J, the derivatives satisfy dvec (G) = (D J G) dvec (J) + (D J ⋆ G) dvec (J ⋆ ) and what we need is (D J G).