{'fraction_non_alphanumeric': 0.035703197764669355, 'fraction_numerical': 0.009081030735796337, 'mean_word_length': 4.287238407878539, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Efficient Computation of Prolate Spheroidal Wave Functions in Radio Astronomical Source Modeling

The application of orthonormal basis functions such as Prolate Spheroidal Wave Functions (PSWF) for accurate source modeling in radio astronomy has been comprehensively studied. They are of great importance for high fidelity, high dynamic range imaging with new radio telescopes as well as conventional ones. But the construction of PSWF is computationally expensive compared to other closed form basis functions. In this paper, we suggest a solution to reduce its computational cost by more efficient construction of the matrix kernel which relates the image domain to visibility (or Fourier) domain. Radio astronomical images are mostly represented using a regular grid of rectangular pixels. This is required for efficient storage and display purposes and moreover, comes naturally as a by product of the Fast Fourier Transform (FFT) in imaging. We propose the use of Delaunay triangulation as opposed to regular gridding of an image for a finer selection of the region of interest (signal support) during the PSWF kernel construction. We show that the computational efficiency improves without loss of information. Once the PSWF basis is constructed using the irregular grid, we revert back to the regular grid by interpolation and thereafter, conventional imaging techniques can be applied.


I. INTRODUCTION
For radio interferometric imaging, we measure the spatial coherency or Fourier components of the celestial sources on different baselines in our array (also called as visibilities). Due to a finite number of receivers, we only sample the Fourier components at a finite number of positions. After calibration or correction for the corruptions due to the propagation path and instrumental effects, we have to deconvolve the source flux from the antennae response pattern to obtain its true intensity map or its radio image.
Conventionally, CLEAN [1], [2] algorithm has been used for modeling and removal of the known sources such that the faint unknown sources in the background can appear. It is a pixel based deconvolution approach. This means any point in the image is associated with a weighted delta function. CLEAN will fail to perform satisfactorily, if two point sources are too close or if the source is partially resolved, as it has analytically been proven in [3]. Moreover, we have shown that making image grid smaller (expecting that accuracy of modeling extended sources improves), is a futile effort by presenting the Cramer Rao Bound (CRB) analysis [3]. Instead, we showed that using a 2-D orthonormal basis for deconvolution achieves the theoretical noise bound or CRB, regardless of the source structure [3].
Cartesian shapelets have been used for source modeling in [3], (and references therein). It has been proven that shapelets provide a way to model extended sources by a closed form construction of an orthonormal basis. But this technique requires additional constraints such as information theoretic bounds to limit the model order. In [4], [5], we have shown that we can construct the optimal basis to model the source by minimum number of basis functions and less artifacts outside Region Of Interest (ROI) using prolate spheroidal wave functions. In addition, unlike shapelets, they do not require additional explicit constrains from information theory, as it is implicitly considered in their mathematical derivation (see [5]). However, the construction of PSWF is computationally more expensive as compared to shapelets [4].
In this paper, we present computationally efficient construction of PSWF suitable for modeling extended sources. We achieve this by downsampling both the Fourier plane visibilities as well as image pixels. Most astronomical (as well as other) images are represented by a regular 2-D array of rectangular pixels. In contrast, we propose to use Delaunay triangulation to construct an irregular grid with fewer number of pixels outside ROI and finer selection of the image inside ROI to represent the image during the construction of PSWF. We can control the size of the triangulation without losing information. Once the PSWF basis has been constructed in this manner, we revert back to the original pixel grid by interpolation.
Notation: We denote vectors in bold lowercase and matrices in bold uppercase. The matrix transpose, Hermitian, and pseudoinverse are denoted by (.) T , (.) H and (.) † respectively. Operation . gives the Frobenious norm of a matrix. In radio interferometric imaging (see Fig. 1), the relation between the intensity of the q-th pixel, f (l q , m q ) in an image and the p-th point, f (u p , v p ) in the visibility domain is described by van Cittert-Zernike theorem [6] as:

II. MATHEMATICAL PRELIMINARIES II-A. Radio interferometric imaging
where N a is the total sampling points in visibility domain. N is the total number of pixels in the image. If we rewrite them in vectorized form, we will obtain: The matrix T (size N a × N) contains e − j2π(l q u p +m q v p ) on its q-th row and p-th column.

II-B. Prolate spheroidal basis
For a complete mathematical derivation of the PSWF, the reader is referred to [5] and references therein. In the following, we will mention a few key mathematical relations of PSWF derivation for stating our problem.
We aim to obtain representative basis, p(l, m) for a source that maximizes the energy in our ROI i.e. within the boundary by which a source is recognized in an image. The larger the ratio given in (3), the less artifacts will remain outside the ROI: The vectorized form of which can be written as: is a selection matrix to select the pixels that belong to the ROI. N b is the number of pixels in the ROI. Then the PSWF basis can be obtained by eigendecomposition of the kernel, K (size N b × N b ) given as: Most of the computations are required first to find the pseudo-inverse term, (TT H ) † (size N a × N a ) and secondly, to find the eigenvalue decomposition of K (size N b × N b ). Note that the rank of (TT H ) † is at most N and usually N a ≫ N > N b . By properly downsampling both the Fourier plane sampling points (reducing N a ) as well as the image pixels (reducing N and N b ) we can significantly cut down the cost of computations. In the next section, we will explain how it has efficiently been done, followed by an example.

III. COMPUTATIONAL COST REDUCTION
We take several steps to reduce the computational cost of (5). The main criterion for this reduction is given by the Landau-Pollak theorem [7]. In its simplest form (applicable to radio interferometry) [4], Landau-Pollak theorem states that the number of degrees of freedom of any given source with compact support, N D could be no greater than the product of area of the source, A lm and the area of the source observed in the Fourier plane, A uv i.e. N D ≤ A lm × A uv .

III-A. Reducing N a
Since we already have an irregular set of N a sampling points on the Fourier plane, it is straightforward to downsample this to a lower value N ′ a , which is already presented in [5]. However, we can do even better by combining this with imaging weights. Each point on the Fourier plane, (say p) has a weight associated with it, ρ(u p , v p ) ∈ [0, 1].
Prior to taking the Fourier transform in (1), we multiply each data point f (u p , v p ) by this weight. When we downsample, we generate a uniformly distributed random number r(u p , v p ) ∈ [0, 1] and if r(u p , v p ) ≤ ρ(u p , v p ) we select the p-th point.

III-B. Reducing N but finer selection of ROI
The motivation behind reducing the number of image pixels is that not all the pixels are required to define the ROI (because of Landau-Pollak theorem) while we construct the PSWF. In Fig. 2 (a), we have shown the ROI of the source as well as the Point Spread Function (PSF). The amount of information can be loosely translated to the number of PSFs that we can pack within the ROI, as shown in Fig. 2 (b). The optimal way to pack the PSFs within the ROI boils down to a circle packing problem, which is generally understood as NP-hard. Therefore, we follow a heuristic approach: Instead of solving a packing problem, we construct an irregular grid of triangles to cover the ROI as well as the area of the image outside the ROI. Obviously, we employ Delaunay triangulation [8] for this purpose. We shall illustrate this with an example as follows: In Fig. 3, we have shown an image of an extended radio source. This image (dimensions 118 × 64) has 7552 pixels.
In Fig. 4, we have shown the ROI for this source, which has 1826 pixels to represent the amount of information present. Therefore, we have generated the Delaunay triangulation for the ROI and its exterior as shown in Fig. 5 such that the triangle size within the ROI is smaller than the exterior triangles. We take the centroids of the triangles as our reduced grid to generate the PSWF. Therefore, we end up with a total number of pixels of 1220 while the ROI has 1080 pixels.

III-C. Reducing cost of pseudoinverse
The cost of (5) mostly involves evaluation of the pseudoinverse (TT H ) † . Instead of directly evaluating this, we    Fig. 3. The triangulation inside the ROI has 1080 triangles while the triangulation outside the ROI has 140 triangles. The triangle size within the ROI is made smaller than the exterior triangles. solve a set of linear equations such as (TT H + γI)x = Ax = y (6) where γ is a small positive regularization parameter. Depending on the context, we are given y and we need to find x. The solution of the linear system for x is done using Q-less QR factorization. Given QR = A (where Q is unitary and R is upper triangular), we first solve for x and make one step correction to get the solution x.

IV. RESULTS
To demonstrate how our solution works in practice, we give results on the construction of PSWF for the extended source shown in Fig. 3. The source was observed by the LOFAR (http://www.lofar.org) radio telescope during system testing. For an 8 hours of observation, we obtain, N a = 7896636 data points on the Fourier plane. We reduce this to N ′ a = 11041 points by applying our downsampling scheme as explained in section 3.1.
We use [8] to generate conforming Delaunay triangulations for the ROI and the exterior of the image. Selection of triangle scales for inside and outside of the ROI (see Fig.5), has to be done in such a way that the least amount of information is lost. For this observation, the PSF minor diameter is about d = 4.2 × 10 −5 radians. To achieve the most appropriate triangulation, we have set different scales for the triangles inside the ROI (a × d) and outside the ROI (b × d), where a, b are the parameters that we vary.
We show the eigenvalue spectrum of the kernel K, for different values of a and b in Fig. 6. The eigenvalue spectrum gives us a measure of information a certain basis can represent. By the results shown in Fig. 6, we see that in most settings, we get about 100 basis functions that carry significant information, which is close to the limit given by Landau-Pollak theorem for this case.
Once we obtain the basis functions for the downsampled grid, we interpolate back to the original image grid to get the basis vectors corresponding to the original image. Let P be the matrix whose columns correspond to the interpolated basis vectors. Then, we study the Frobenious norm of P H P−I and the condition number of P H P in Fig. 7. This, in addition to the eigenvalue spectrum in Fig. 6 give us a benchmark to choose the most suitable triangulation (or values of a and b). If the basis vectors are orthonormal, we will obtain the lowest value for the Frobenious norm. Fig. 7 (a) shows that the optimality of the source modeling is not affected significantly from one setting to another, because the basis vectors are orthonormal with negligible errors. Thus, the CRB remains minimum [3]. Moreover, the condition number of P H P helps us to understand how its condition affects the further computation. Ultimately, these conclude that the setting with a = 0.2 and b = 2.1 is the most suitable one to fulfill the aforementioned criteria. By this setting, we downsample to N ′ = 3270 pixels in which the ROI has N ′ b = 3199 pixels. In contrast, the original pixel grid has N = 7552 and N b = 1826.

V. CONCLUSIONS
The use of PSWF for source modeling and deconvolution of extended sources was proposed earlier in [4] and [5]. In this paper, we suggested a scheme to reduce its computational cost. By using the Delaunay triangulation, we can significantly reduce the cost of construction of PSWF. We have also shown by a real example that we do not lose information by this cost reduction. Further improvement will focus on finer selection of the ROI depending on the source structure such that an area with higher intensity is represented with more triangle pixels.