Optimal Probabilistic Catalogue Matching for Radio Sources

For any given hypothesis, we introduce the true position of the assumed latent object (i.e., the host galaxy). We consider every possible position of the host galaxy relative to the radio components, and let the available data (measured positions and their uncertainty) provide appropriate weighting via the likelihood function, e.g., a Gaussian. The computational methods for efficient evaluations are discussed in Fan et al. (2015) in detail.

The updated model is illustrated in Figure 1 in the tangent plane. For the calculations, we chose the position of the infrared detection as the origin. In this 2D coordinate system, the radio core is at $\boldsymbol{\omega}$ and the lobes are at $\boldsymbol{\omega}^{\prime}$ and $\boldsymbol{\omega}^{\prime\prime}$. The following equations describe their relations to each other using additional model parameters $\phi,k$: $k$ allows for the lobes to be at different distances from the core, and $\phi$ permits (slight) bending.

where $R_{\phi}$ is the usual 2$\times$2 matrix of the rotation by $\phi$. All model parameters $\theta=\left\{\boldsymbol{\omega},\boldsymbol{\omega}^{\prime},k,\phi\right\}$ are involved in the assessment of the associations along with their priors, which we can write as

due to the appropriate (in)dependencies. Hereafter we assume that $\rho(\boldsymbol{\omega})$ is isotropic (or uniform on the surface of the sphere) and $\rho(\boldsymbol{\omega}^{\prime}|\boldsymbol{\omega})$ depends only on the separation of the two positions. Furthermore, we will assume simple Gaussian priors for $k$ and $\phi$ with 0 mean, a combination that corresponds to a straight and symmetric configuration of the two lobes.

Following the approach of Budavári & Loredo (2015) and Shi et al. (2019) we ideally consider all possible matched catalogues with every possible combinations of associated detections. Let $(i,c)$ represent the $i$th detection (or an internal identifier) in catalogue $c$ with measured position $\boldsymbol{x}_{ic}$. The collection of all detections $D$ is the set of all $(i,c)$ measurements. We adopt the use of the Fisher (1953) distribution for the likelihood calculation,

that essentially corresponds to a Gaussian in the usual flat sky limit for small uncertainty, where the compactness parameter $\kappa=1/\sigma^{2}$, and the $\sigma$ is the positional error from the catalogue.

The final cross-matched catalogue will consist of multiple associations, which we will call objects, and index them by $o$. A particular object $o$ will have an $S_{o}$ set of $(i,c)$ detections in its association. In fact, we can think of matching radio observations as if the radio catalogue appeared multiple times according to what roles the radio detections play in a particular association: a core or one of the two possible lobes. So, for matching IR and radio observations, we order the collections as follows: $c\!=\!0$ corresponds to the IR catalogue, $c\!=\!1$ represents radio core detections, and $c\!=\!2,3$ correspond to lobes. Conceptually, the catalogues corresponding to $c\!=\!1$, $2$, and $3$ are identical copies of the radio observations. In practice, however, they really are just views of the catalogue; there is no need for duplication of the data. Naturally, a given radio detection can only appear in one association $o$, which we will enforce explicitly when looking for the optimal catalogue.

If the measurements in $S_{o}$ correspond to one object, the marginal likelihood of that $o$ is calculated by the integral

where the prior $\rho_{o}$ has an index to signify the possibility that different catalogues could have different sky coverage, hence different overlapping regions would yield different priors, see discussion in Budavári & Szalay (2008) and Budavári & Loredo (2015). The member likelihoods are simply given by the Fisher distribution at the appropriate true positions,

where $\boldsymbol{\omega}^{\prime\prime}$ is really a function of all parameters, as given by eq.(1) above. To find the optimal catalogue, we need to find the set of objects such that the product of their marginal likelihoods, $\prod\mathcal{M}_{o}$, is maximal.

In addition, we also introduce the marginal likelihood of “no association”, where each detection is considered by itself as a separate object. The integrals simplify to the familiar product

with which we can define a Bayes factor for object $o$ as

Considering that the product of the marginal likelihoods $\prod\mathcal{M}_{o}^{\textrm{NA}}$ over all $o$ in a matched catalogue always contains the same terms corresponding to every $(i,c)$ detection, the product is constant across every possible partitioning. This means that the optimal catalogue will also have maximal $\prod B_{o}$ value. Using the latter over the product of marginals has the advantage that $B_{o}\!=\!1$ for the orphan objects that consist of a single detection, and hence their contribution to the product can neglected.

To formalize the above maximization of the global cross-match catalogue likelihood (Budavári & Basu, 2016), one first has to enumerate all possible associations. Following Shi et al. (2019) we index these by $T$ and introduce a binary variable $x_{T}$ that is 1 if the association is selected for the cross-match catalogue and 0 otherwise. Since a radio detection should be used only from exactly one of the three catalogues $c=1,2,3$, we must not have any variables that correspond to a set of detections where the same detection is used from more than one catalogues $c\in\{1,2,3\}$. More precisely, we do not include a variable for any set $S_{o}$ which contains $(i,c)$ and $(i^{\prime},c^{\prime})$ such that $i\!=\!i^{\prime}$ but $c\!\neq\!c^{\prime}$.

The Bayes factors can be (numerically) evaluated for every one of the candidates, $B_{T}$. Formally, the optimization over $x\!=\!\{x_{T}\}$ can now be written as

where $w_{T}\!=\!-\log B_{T}$, but further constraints are required to ensure that every detection appears exactly once in the resulting catalogue, and each radio detection appears only in one of the three possible roles ($c\!=\!1$, $2$, $3$). Considering that orphans have no contributions to the above objective as their $w_{T}\!=\!0$, we can write inequalities that are better handled by optimization algorithms. For example,

means that every $(i,c)$ detection can appear at most in one $T$ association of a matched catalogue: the sum goes over all $T$ that include $(i,c)$. For the radio detections, we have further constraints to ensure that each and every one of them only appears (at most) in one role, i.e., $c=1$, $2$ or $3$, which we can write as

where $R$ is the set of indices of the detections in the radio catalogue, i.e., $R\!=\!\left\{\,i:(i,\!1)\!\in\!D\,\right\}$. Such an optimization problem falls in the realm of Integer Linear Programming (hereafter ILP) for which various commercial and research solvers are available off the shelf.

The above optimization considers only the marginal likelihoods of the associations, which neglects any prior information about the resulting catalogues. If we had such knowledge, the objective could be easily updated to include that. For example, one might have a distribution of the expected number of objects in a catalogue, or perhaps know the relative frequency of occurrence for the different kinds of matches, such as radio triples, doubles, etc. We expect that with more data, a hierarchical modeling approach will be eventually developed to simultaneously learn the population distributions and the assignments, but for now we will resort to the simplest possible scenario: In order to maximise the number of radio components in each association, we minimize the number of final associations. In case that, a group of three radio components is more likely to be associated with a single galaxy, than three separate galaxies, if these two scenarios each have the same calculated probability.

Our modified optimization would now have two terms,

where $\lambda$ expresses the strength of our preference for fewer objects. This objective is also linear, and hence does not add to the complexity of the problem. Our input data consists only of the directions of the detections, with inherent uncertainties, and ignores other information that could in principle be used as constraints. For example, we expect that additional shape and flux measurements will play a key role in defining scientific catalogues, but these need to be studied further with more data. With these caveats and ideas in mind, we will explore solutions using the multi-objective approach, as in eq.(11).

The cross-matching process consists of several steps. The first step is to remove highly unlikely associations, using a spatial two-way matching between the IR and radio catalogues. This helps in the efficiency of the method by reducing the number of variables $x_{T}$ in the optimization (see the discussion in Section 2.3). In particular, when two measurements are far from each other on the sky, they are unlikely to be related. We set the angular separation threshold for the preprocessing to 2′, which is a safe limit considering that the farthest components in the radio triples associated by Norris et al. (2006) are 61.5″apart.

In the second step, all possible hypotheses are enumerated as possible candidates for associations. When considering two lobes we break the degeneracy by considering the farthest lobe to be at $\boldsymbol{\omega}^{\prime}$ and the closer one at $\boldsymbol{\omega}^{\prime\prime}$, then we evaluate the marginal likelihoods for all combinations of detections. For the member likelihood functions, we choose $\sigma_{0}\!=\!0.2$″ astrometric uncertainty for SWIRE, $\sigma_{1}\!=\!1.2$″ for the radio cores and $\sigma_{2}\!=\!\sigma_{3}\!=\!2.2$″ for the lobes, whose position is harder to determine due to their more extended and somewhat asymmetric nature.

The prior probability density function in eq.(2) is assumed to be isotropic over $\boldsymbol{\omega}$ (see discussions in Budavári & Szalay (2008) and Budavári & Loredo (2015)). The conditional PDF $\rho(\boldsymbol{\omega}^{\prime}|\boldsymbol{\omega})$ is assumed to be a function of only the angular separation of $\boldsymbol{\omega}$ and $\boldsymbol{\omega}^{\prime}$, and its radial dependence has the shape of the Rayleigh distribution with $\sigma_{R}=8.5"$. The dimensionless parameter $k$ that describes the difference in the lobe separations from the core is assumed to have a Gaussian prior with a standard deviation of $\sigma_{k}\!=\!1$. The angle $\phi$ that describes how bent the geometry is also assumed to be a Gaussian with 0 mean and $\sigma_{\phi}\!=\!15^{\circ}$ degrees. The 15 degrees is training from Norris et al. (2006), of which maximal angle is 16.96 degrees.

With those in hand, our implementation numerically integrates the marginal likelihoods and computes the Bayes factors. For an efficient integration, we sample $\boldsymbol{\omega}$ from the IR core’s likelihood function and $\boldsymbol{\omega}^{\prime}$ from the farthest radio lobe’s likelihood. Given those and random samples from the Gaussian priors for $k$ and $\phi$, we can generate a random $\boldsymbol{\omega}^{\prime\prime}$ and evaluate the contribution of that sample to the integral. Figure 2 illustrates 10,000 sample points for a known radio triple, whose positions are shown by large blue crosses. Red contours around (0,0) represent the density of $\boldsymbol{\omega}$ samples, which are highly concentrated due to the small uncertainty of the IR source, and therefore appear as a red dot at the origin. The core radio component (the blue cross near the origin) is located a small distance from the origin due to its larger uncertainty – its likelihood is evaluated at $\boldsymbol{\omega}$. The black contours and the yellow points toward the bottom of the figure show the $\boldsymbol{\omega}^{\prime}$ samples from the likelihood function of the radio lobe. The contours on top for $\boldsymbol{\omega}^{\prime\prime}$ show the samples that are a combination of all model parameters, $\boldsymbol{\omega}$, $\boldsymbol{\omega}^{\prime}$, $k$ and $\phi$.

When every $w_{T}\!=\!-\log B_{T}$ is computed, we solve the minimization using the ILP implementation by Gurobi Optimization (2019) for the binary variables $x_{T}$ using the aforementioned objectives and constraints. Gurobi has builtin routines to handle the multi-objective optimization of eq.(11), for which we adopt $\lambda\!=\!2.5$ to prefer solutions with fewer objects, i.e., larger associations of detections. Considering that the average weight $\bar{w}$ is significant for selected associations, the weighting effectively is by a $\lambda^{\prime}=\lambda/\lvert\bar{w}\rvert$ which in our experiments is about one fifth.

First, some caveats: Norris et al. (2006) associate the SWIRE and ATLAS measurements by eye. Expert radio astronomers look at the IR images on top of which multiple contours show the details of radio intensity, and they annotate the catalogues based on such rich views of the data. In comparison, our current approach considers only the positions of the detections and their uncertainty but not their brightness nor their shapes. We consider the expert matched catalogue as our ground truth, hereafter called the “reference catalogue”, to which we compare our results. We hope that our approach can approximately recover many of the known associations without introducing too many spurious matches. However, we also caution that in some cases the classifications are ambiguous, and even experts will sometimes disagree over the correct classification, so the “reference catalogue” is not always reliable. The reference catalogue primarily consists of IR detections associated with radio triples, doubles and singles but also lists six complex associations, see Table 2, which do not fit the usual model. We will discuss all these scenarios in turn.

Our implementation produces an optimal cross-matched catalogue under the morphology model described in Section 2. Table 1 lists the best solution along with the number of discovered associations with radio triples, doubles and singles. In fact, we separate them into three categories based on whether they appear in the reference catalogue as common, missed, or extra. We see that our code finds all 10 of the radio triples in the reference catalogue with only two extra associations. Out of the 25 radio doubles we find 21 with an additional 10 extra. To determine how good these results really are, we now examine in turn the associations with radio triples, doubles and singles.

Since the code finds all radio triples in the reference catalogue, we study the two extra associations that our approach also identifies:

1. 1.

   The radio components C148–C151–C141 (shown in Figure 3) associated with SWIRE3_J032825.92-271701.3 also appear in the reference catalogue but not as a triple, because there is no central radio component associated with an IR source. Instead they form complex #1 with the comment: “Radio double with connecting jet”. The ILP has identified the IR source SWIRE3_J032825.92-271701.3 as the host galaxy. While this is well separated from the radio component C148, there is low-level radio emission associated with this SWIRE source, so it possible that this is the correct host galaxy, even though it does not seem to be associated with the listed central radio component.

2. 2.

   The radio components C440–C446–C437 associated with SWIRE3_J033210.74-272635.5 shown Figure 4 also appear in the reference catalogue as complex #3 with the comment “weird complex source with tails. z(g)”. It was not classified as a triple in the reference catalogue because to do so would ignore the bent tails, probably resulting from an interaction with an intra-cluster medium. However, here we focus only on the positions of the components, and in those terms this source is correctly classified as a triple.

So, we conclude that, provided only the positions of the components are considered, one of the two extra sources is correctly classified as a triple. So our procedure associates radio triples with a 92% success rate, with the only failure being a source that closely resembles a triple, although lacking a central core.

Our procedure correctly classified 21 of the 31 doubles, missed 4, and incorrectly identified a further 10 sources as doubles. However, we caution that, as discussed in Section 3.2 above, the reference catalogue can occasionally be unreliable. In some cases below, it is hard to tell, with the given data, whether the Reference catalogue or the ILP is more likely to be correct.

Our procedure misses four of the radio doubles in the reference catalogue: C089–C091, C393–C398, C500–C501, and C573–C574, as shown in Figure 5 .

C089–C091 with SWIRE3_J032731.61-275245.5 is classified in the reference catalogue as “complex region, but probably core-jet with C091 core, C089 jet”; ILP classified C089–C091 with the other SWIRE3_J032730.87-275247.6 as a “LOBE–LOBE”.

C393–C398 with SWIRE3_J033145.54-281955.0 is classified in the reference catalogue as “Probably asymmetric radio double”; ILP classified C393–C398 with the other SWIRE3_J033145.28-281956.1 as a “CORE–LOBE”.

C500–C501 with SWIRE3_J033242.82-273817.6 is classified in the reference catalogue as “Radio double”; ILP classified C500–C501 with the other SWIRE3_J033242.60-273816.0 as “LOBE–LOBE” (i.e. a double), so the only difference here is the choice of the host galaxy.

Interestingly, the first three of them actually appear as extra doubles, but in associations with different SWIRE measurements. Upon visual inspection, these mismatches appear to be a result of limited information in the positions only considered in this paper. It is possible that including brightness or shape measurements might help find the correct associations, but this needs to be studied further.

The fourth double of C573–C574 also appears in our catalogue, but it is separated into two singles due to two SWIRE detections close to the radio components, see bottom left image of Figure 5. The reference catalogue comment reads “M-test classifies it as double with C573, but C574 has good ID so C573 is probably a jet”.

Our code produces seven extra radio doubles: six of them were classified as “complex” as shown in Table 3. For example, #6 was classified as a double in the comments to the reference catalogue, and has been correctly identified as such by the ILP, but it also contains a third component, C774, as an extension to one of the lobes, and is therefore classified as complex in the reference catalogue rather than as a simple double. Thus the ILP’s classification as “double” in this case is possibly a more useful classification than “complex”.

The remaining C517–C514 is classified by ILP as a Core-Lobe configuration, but the reference catalogue lists them as two singles, as shown the bottom right image of Figure 5. Given the astrometric uncertainty, it is difficult to tell whether C517 and C514 are associated with the IR sources, so in this case it is hard to tell whether the disagreement represents an error in ILP or an error in the reference catalogue.

Singles are individual radio components associated with a SWIRE detection. Out of the 685 single associations in the reference catalogue, our approach identifies 680, with 4 extra associations.

Out of the five missed (C381, C475, C514, C517, C763), three C381, C514, C517 are listed in the reference catalogue as parts of doubles. The reference catalogue comment for C475 (see Figure 6) is “strong spitzer ID at one end of core-jet src.”, but this source is too far from that SWIRE source to pass our detection threshold (see Section 3.1), resulting in a low Bayes factor for this association.

Finally, C763 (see Figure 6) is associated with a different SWIRE object from the reference catalogue, so it is listed as an extra. The ILP associates the radio source with a weak SWIRE detection, instead of a brighter source farther away. It is interesting that the reference catalogue identification is associated with the peak of the contours, but the ILP association is with the catalogue position, which is significantly different from the peak of the contours, suggesting an error in the reference catalogue.

There are 4 extra singles, most of them are mentioned above, except for C774, which is part of the complex #6, which is discussed above in Section 5.

As discussed in Section 2.2. the ILP technique generates many slightly different hypotheses, just one of which is chosen as the optimal solution.

In addition to the best solution, our approach can also create similarly good but slightly sub-optimal associations. We studied the first 3 solutions, which show surprisingly little variation. The difference is due to the assignments of C759, C760, C761, C764 and C768, which we illustrate in Figure 7. These actually form complex #5 with comment “Centre of a linear complex (jets or grav arc?)”.

The differences are as follows: In the best solution, C759–C760 appears as a Lobe-Lobe combination and C761–C764 is a Core-Lobe with C768 being single. In the second solution, lists C761–C759–C764 as a radio triple with singles C760 and C768. In the third solution, there C761 and C759 do not appear, and C764–C768 is a Lobe-Lobe, with single C760, as summarized in Table. 4

Our algorithm will, of course, perform correctly only if there are no errors in the catalogue. In particular, we assume that the catalogue lists one position per radio component, which is true for the data being considered here, and generally true for the source finder (Selavy) used for the EMU surveys (Norris et al., 2011). However, we are aware that some source finders can generate several Gaussian components for one radio component. While we have not tried the ILP algorithm on catalogues produced by such source finders, we suspect a pre-processing stage may be need to consolidate the listed components into one per radio component.