2.1. Data description

As noted in Section 1, most ML-based techniques are typically focused on optically selected datasets, primarily based around the SDSS datasets, providing high source counts of stars, galaxies, and QSO with photometry (generally) in u, g, r, i, and z bands, with a spectroscopically measured redshift—generally using the SDSS Galaxy or QSO datasets, shown in Fig. 1. In this work, the data are selected specifically to better represent the data expected from the upcoming EMU Survey (Norris et al. Reference Norris2011) and Evolutionary Map of the Universe–Pilot Survey 1 (EMU-PS, Norris et al. Reference Norris2021). Towards this end, we only accept SDSS objects with a counterpart in the NRAOVLA Sky Survey (NVSS) or Faint Images of the Radio Sky at Twenty-Centimeters (FIRST) radio surveys.

Figure 2.

Cross-match completed between the NVSS Radio sample and the AllWISE Infrared sample. The blue line represents the straight nearest-neighbour cross-match between the two datasets, and the orange line represents the nearest-neighbour cross-match where 1 $^{\prime}$ has been added to the declination of every radio source. The vertical black line denotes the chosen cutoff.

This work compiles three datasets, each containing multiple features for comparison:

1. 1. Northern Sky—Radio Galaxy Zoo (RGZ) and NVSS

   1. (a) Our Northern Sky dataset contains two radio samples—the NVSS sample and the RGZ FIRST-based sample (where the RGZ sample has been cross-matched with the AllWISE sample, explained in Banfield et al. Reference Banfield2015 and Wong et al. in preparation).

   2. (b) The NVSS sample was cross-matched with AllWISE at 4 $^{\prime\prime}$ , providing 564799 radio sources—approximately 32% of the NVSS sample—with an infrared counterpart. The NVSS/AllWISE cross-match has an estimated 7%—123484 sources—misclassification rate, where the misclassification rate is quantified by shifting the declination of all sources in the NVSS catalogue by 1 $^{\prime}$ and re-cross-matching based on the new declination, following the process described in Norris et al. (Reference Norris2021). Fig. 2 shows the classification/misclassification rates as a function of angular separation and is used to determine the optimum cross-match radius.

   3. (c) The NVSS and RGZ samples were then combined, removing duplicates based on the AllWISE unique identifier, providing 613551 radio sources with AllWISE detections.

   4. (d) The northern sky radio/infrared catalogue was then cross-matched against the SDSS catalogue (providing both optical photometry and spectroscopic redshifts) based on the infrared source locations—a radio-infrared cross-match tends to be more reliable, when compared with the radio-optical cross-match (Swan Reference Swan2018)—at 4 $^{\prime\prime}$ , providing a classification/misclassification rate of 9.33%/0.06% (55716/348 sources) (Fig. 3).

   5. (e) Finally, all sources from the Stripe82 Equatorial region were removed (sources with an RA between $9^{\circ}$ and $36^{\circ}$ or $330^{\circ} $ and $350^{\circ}$ and DEC between $-1.5^{\circ}$ and $1.5^{\circ}$ ).

   6. (f) The final Northern Sky Sample contains 55452 radio-selected sources with a spectroscopically measured redshift, SDSS g, r, i, and z magnitudes measured using Model/PSF/Fibre systems, and AllWISE W1, W2, W3, and W4 infrared magnitudes, shown in Table 1.

      Table 1.

      The source count in each sample compiled in Section 2.

      
      Figure 3.

      Similar to Fig. 2, cross-matching the Northern Sky sample based on the NVSS and RGZ radio catalogues and AllWISE Infrared, with SDSS optical photometry and spectroscopic redshift. Note the different scale of the right-side y-axis.

      

2. 2. Southern Sky—ATLAS

   1. (a) Beginning with the ATLAS dataset—described in Luken et al. (Reference Luken, Norris, Park, Wang and Filipović2022)—we cross-match the SpitzerWide-Area Infrared Extragalactic Survey (SWIRE) infrared positions with AllWISE in order to gain the same infrared bands as the Northern Sky dataset. Cross-matching at 1 $^{\prime\prime}$ produces a 100%/0.34% classification/misclassification rate (1156/4 sources) (Fig. 4), and a final source count of 1156 sources, all with g, r, i, and z optical magnitudes in Auto, 2 $^{\prime\prime}$ , 3 $^{\prime\prime}$ , 4 $^{\prime\prime}$ , 5 $^{\prime\prime}$ , 6 $^{\prime\prime}$ , and 7 $^{\prime\prime}$ apertures, as well as the W1, W2, W3, and W4 infrared magnitudes.

      Figure 4.

      Similar to Fig. 2, cross-matching the southern sky ATLAS sample with the AllWISE catalogue, matching the SWIRE positions with the AllWISE positions.

      
      Figure 5.

      Similar to Fig. 2, cross-matching the DES optical photometry with the SDSS optical photometry and spectroscopic redshift catalogues.

      

3. 3. Equatorial—Stripe82

   1. (a) Along the equatorial plane, the Stripe82 field has been extensively studied by both northern- and southern-hemisphere telescopes, providing a field that contains both SDSS, and DES photometry. Cross-matching the DES catalogue with the SDSS catalogue (where both catalogues were restricted to the Stripe82 field) at 1 $^{\prime\prime}$ produces a 98.4%/3.36% (170622/5831 sources) classification/misclassification rate (Fig. 5).

   2. (b) The optical catalogues were then cross-matched against the AllWISE catalogues at 1.25 $^{\prime\prime}$ producing a 84.09%/ 0.60% (129837/932 sources) classification/ misclassification rate (Fig. 6).

   3. (c) Finally, cross-matching the AllWISE Infrared catalogue against the Hodge et al. (Reference Hodge, Becker, White, Richards and Zeimann2011) (at 4 $^{\prime\prime}$ ; see Fig. 7) and Prescott et al. (Reference Prescott2018) (at 4 $^{\prime\prime}$ ; see Fig. 8) gives us a 21.96%/ 0.42% (3946/75 sources) and 45.46%/0.54% (2180/26 sources) classification/misclassification rates, respectively.

   4. (d) After combination of the Stripe82 Radio datasets with duplicates removed (based on AllWISE ID), we have a final dataset of 3030 radio-selected sources with a spectroscopic redshift, W1, W2, W3, and W4 infrared magnitude, g, r, i, and z optical magnitudes in PSF, Fibre, and Model systems (for SDSS photometry), and Auto, 2 $^{\prime\prime}$ , 3 $^{\prime\prime}$ , 4 $^{\prime\prime}$ , 5 $^{\prime\prime}$ , 6 $^{\prime\prime}$ , and 7 $^{\prime\prime}$ apertures (for DES photometry).

To summarise, all datasets are radio-selected, and contain:

Figure 6.

Similar to Fig. 2, cross-matching the DES positions with the AllWISE Infrared catalogue. Note the different scale of the right-side y-axis.

Figure 7.

Similar to Fig. 2, cross-matching the AllWISE positions with the Hodge et al. (Reference Hodge, Becker, White, Richards and Zeimann2011) Radio catalogue. Note the different scale of the right-side y-axis.

Figure 8.

Similar to Fig. 2, cross-matching the AllWISE positions with the Prescott et al. (Reference Prescott2018) Radio catalogue. Note the different scale of the right-side y-axis.

Figure 9.

Histogram showing the density of sources at different redshifts in the combined RGZ—North, Stripe82—Equatorial, and ATLAS—South—(blue), and the Square Kilometre Array Design Survey (SKADS, Levrier et al. Reference Levrier2009) simulation trimmed to expected EMU depth (Norris et al. Reference Norris2011).

• • a spectroscopically measured redshift, taken from either the Australian Dark Energy Survey (OzDES), or SDSS;

• • g, r, i, and z optical magnitudes, taken from either the DES or SDSS;

• • W1, W2, W3, and W4 (3.4, 4.6, 12, 24 $,\unicode{x03BC}\mathrm{m}$ , respectively) infrared magnitudes, taken from AllWISE;

with a final redshift distribution shown in Fig. 9. While we are still not matching the expected distribution from the EMU survey, we are ensuring all sources have a radio counterpart (and hence, will be dominated by the difficult-to-estimate AGNs), with the final distribution containing more, higher redshift radio sources than previous works like Luken et al. (Reference Luken, Norris, Park, Wang and Filipović2022).

The primary difference between the datasets is the source of the optical photometry. Even though both the DECam on the Blanco Telescope at the Cerro Tololo Inter-American Observatory in Chile and the Sloan Foundation 2.5 m Telescope at the Apache Point Observatory in New Mexico both use g, r, i, and z filters, the filter responses are slightly different (demonstrated in Fig. 10, the DES Collaboration notes that there may be up to 10% difference between the SDSS and DES equivalent filtersFootnote c), with different processing methods producing multiple, significantly different measurements for the same sources. For ML models, a difference of up to 10% is significant and had significant effects on redshift estimations in early tests without correction (sample results with one ML algorithm shown in Appendix A).

Figure 10.

A comparison of the g, r, i, and z filter responses, used by the DES (top), and the SDSS (bottom).

For the SDSS, the three measures of magnitude used in this work (Point Spread Function (PSF), Fibre and Model) are all extensively defined by the SDSS.Footnote d Simply put, the PSF magnitude measures the flux within the PSF of the telescope for that pointing, the fibre is a static sized aperture based on a single fibre within the SDSS spectrograph (generally 3 $^{\prime\prime}$ ), and the model magnitude tries to fit the source using a variety of models.

The DES pipelines produce statically defined apertures from 2 $^{\prime\prime}$ to 12 $^{\prime\prime}$ , as well as an auto magnitude that is fit by a model.

For our purposes in finding DES photometry compatible with SDSS photometry, we only examined the DES auto, and 2–7 $^{\prime\prime}$ measurements, as the larger aperture DES measurements begin to greatly differ from any measured SDSS measurement. We find that the DES auto magnitude is most similar to the SDSS model magnitude, and hence, exclusively use this pairing.

2.2. Optical photometry homogenisation

The combined dataset discussed above (Section 2.1) contains optical photometry measured using the SDSS (in the Northern, and Equatorial fields) and the DES (Southern and Equatorial fields). As shown in Fig. 10, while the SDSS and DES g, r, i, and z filters are similar, they are not identical, and hence should not be directly compared without modification before use by typical ML algorithms. As the Stripe82 Equatorial field contains observations with both optical surveys, we can fit a third-order polynomial from the $g - z$ colour, to the difference in the SDSS and DES measured magnitude for each band for each object, and use the fitted model to homogenise the DES photometry to the SDSS photometry for the Southern hemisphere data. Fig. 11 shows four panels—one for each of the g, r, i, and z magnitudes—with the orange points showing the original difference between optical samples against the $g - z$ colour, blue points showing the corrected difference, orange line showing the third-order polynomial fitted to the original data, and the blue line showing a third-order polynomial fitted to the corrected data. While this homogenisation does not adjust for the scatter in the differences, it does shift the average difference, dropping from 0.158, 0.149, 0.061, and 0.006 to 0.004, 0.001, 0.001, and 0.007 for the g, r, i, and z mag, respectively. We explore the difference the corrections make to predicting the redshift of sources with SDSS and DES, using the kNN algorithm trained on the opposite optical survey, using corrected, and uncorrected DES photometry in Appendix A.

Figure 11.

Plot showing the effects of homogenisation on the optical photometry. Each panel shows the original difference between the DES and SDSS photometry for a given band (with the band noted in the title of the subplot), as a function of $g - z$ colour. The orange scatterplots are the original data, with the orange line showing a third-order polynomial fit for to the pre-corrected data. The blue scatterplots are the corrected data, with the blue line showing the post-correction fit, highlighting the improvement the corrections bring.

2.3. Regression and classification

The distribution of spectroscopically measured redshifts is highly non-uniform, providing additional difficulties to what is typically a regression problem (a real value—redshift—being estimated based on the attributes—features—of the astronomical object). As demonstrated in Fig. 1, it also does not follow the expected distribution of the EMU survey, partly because the optical source counts of the local universe vastly outnumber those of the high-redshift universe, and partly because high-redshift galaxies are too faint for most optical spectroscopy surveys. The non-uniform distribution means high-redshift sources will be under-represented in training samples and therefore are less likely to be modelled correctly by ML models.

In an attempt to provide a uniform redshift distribution for the ML methods to provide better high-z estimations, we quantise the data into 30 redshift bins with equal numbers of sources in each (where the bin edges, and the expected value of the bin—typically the median redshift of the bin—are shown in Table 2). While binning the data means that it is no longer suitable for regression, it allows us to use the classification modes of the ML methods and test whether treating the redshift estimation problem as a classification problem rather than attempt to estimate the redshift of sources as a continuous value aids in the estimation of sources in the high-redshift regime.

Table 2.

Example redshift bin boundaries used in the classification tests, calculated with the first random seed used. We show the bin index, the upper and lower bounds, and the predicted value for the bin.

3.1. Error metrics

As stated in Section 1, this work differs from the typical training methods that attempt to minimise the average accuracy of the model (defined in Equation (4), or Equation (5)). Instead, it is primarily focused on minimising the number of estimates that are incorrect by a catastrophic level—a metric defined as the Outlier Rate:

(1) \begin{equation} \eta_{0.15} = \frac{1}{N}{\sum_{z\in Z} [\![ |\Delta z| > 0.15(1 + z_{spec}) ]\!]}, \end{equation}

where $\eta_{0.15}$ is the catastrophic outlier rate, Z is the set of sources, $|Z| = N$ , $[\![ x ]\!]$ is the indicator function (1 if x is true, otherwise it is 0), $z_{spec}$ is the measured spectroscopic redshift, and $\Delta z$ is the residual:

(2) \begin{equation} \Delta z = z_{spec} - z_{photo}, \end{equation}

Alternative, we provide the 2- $\sigma$ outlier rate as a more statistically sound comparison:

(3) \begin{equation} \eta_{0.15} = \frac{1}{N}{\sum_{z\in Z} [\![ |\Delta z| > 2 \sigma ]\!]}, \end{equation}

where $\eta_{2\sigma}$ is the 2- $\sigma$ outlier rate, and $\sigma$ is the residual standard deviation:

(4) \begin{equation} \sigma = \sqrt{\frac{1}{N} \sum_{i=1}^N(y_i - \hat{y}_i)^2}, \end{equation}

where $\sigma$ is the residual standard deviation, $y_i$ is an individual spectroscopic redshift, and $\hat{y_i}$ is the corresponding estimate for source i. The residual standard deviation gives an indication of the average accuracy of the estimates.

The Normalised Median Absolute Deviation (NMAD) gives a similar metric to the Residual Standard Deviation, but is more robust to outliers as it relies on the median, rather than the mean of the residuals:

(5) \begin{equation} \sigma_{\mathrm{NMAD}} = 1.4826 \times (\mathrm{median}(|X_i - \mathrm{median}(X)|), \end{equation}

$\sigma_{\mathrm{NMAD}}$ is the NMAD, X is a set of residuals (where the individual values are calculated in Equation (2) as $\Delta z$ ), from which $x_i$ is an individual observation.

The Mean Square Error (MSE) is only used in regression-based tests, provides the average squared error of the estimates

(6) \begin{equation} \operatorname{MSE} =\frac {1}{N} \sum _{i=1}^{N}(y_{i}-\hat {y_{i}})^{2}, \end{equation}

where $\operatorname{MSE}$ is the Mean Square Error, $y_i$ is an individual spectroscopic redshift, and $\hat {y_{i}}$ is the corresponding estimated redshift for source i.

The accuracy is only used in classification-based tests and provides the percentage of sources predicted in the correct ‘class’, where the class is a particular redshift bin. This metric is provided for completeness only, as the accuracy is only accepting of perfect classifications, whereas the aim of this work is provide redshift estimates that are approximately correct—i.e. we are inherently accepting of classifications in nearby redshift bins, which would be considered incorrect classifications by the accuracy metric.

(7) \begin{equation} \mathrm{Accuracy}(y, \hat{y}) = \frac{1}{N} \sum_{i=0}^{N-1} [\![ \hat{y}_i = y_i ]\!], \end{equation}

where y is a vector of spectroscopic redshifts, and $\hat{y}$ is the corresponding vector of estimated redshifts.

3.2. Statistical significance

In order to measure the potential random variation within our results, all tests were conducted 100 times, with different random seeds—creating 100 different training/test sets to train and test each algorithm on. All values presented are the average of the results gained, with the associated standard error:

(8) \begin{equation} \sigma_{\bar{x}} = \frac{\sigma_{x}}{\sqrt{n}},\end{equation}

where $\sigma_{\bar{x}}$ is the standard error of $\bar{x}$ which is calculated from the the standard deviation of the 100 repetitions of the experiment using different random seeds (denoted as $\sigma_x$ ), $\bar{x}$ is the mean classification/regression error, and n is the number of repetitions—100 in this case.

We note that the classification bin distribution is calculated for each random initialisation—this means that while each of the 100 random training sets will have roughly the same redshift distribution, there will be slight differences in the bin distributions calculated for classification.

3.3. k-Nearest Neighbours

The kNN algorithm is one of the oldest (Cover & Hart Reference Cover and Hart1967), as well as one of the simplest machine learning algorithms. Using some kind of distance metric—typically Euclidean distance—a similarity matrix is computed between every source in the training set, comparing the observed photometry between sources. The photometry of sources in the test set—sources with ‘unknown’ redshift—can then be compared to the photometry in the training set and find the ‘k’ (hereafter $k_n$ ) sources with most similar photometry. The mean or mode (depending on whether regression, or classification is performed respectively) of the most similar sources redshift from the training set is taken as the redshift of the unknown source. Following Luken et al. (Reference Luken, Norris, Park, Wang and Filipović2022) who have shown that Euclidean distance is far from optimal for redshift estimation, here we use the Mahalanobis distance metric (Equation (9); Mahalanobis Reference Mahalanobis1936):

(9) \begin{equation} d(\vec{p},\vec{q}) = \sqrt{(\vec{p} - \vec{q})^\mathrm{T}S^{-1}(\vec{p} - \vec{q})},\end{equation}

where $d(\vec{p},\vec{q})$ is the Mahalanobis distance between two feature vectors $\vec{p}$ and $\vec{q}$ , and S is the covariance matrix.

The value of $k_n$ is optimised using k-fold cross-validation, a process where the training set is split into k (hereafter $k_f$ and is assigned a value of 5 for this work) subsets, allowing the parameter being optimised to be trained and tested on the entire training set.

3.4. Random Forest

The RF algorithm is an ensemble ML algorithm, meaning that it combines the results of many other algorithms (in this case Decision Trees (DTs)) to produce a final estimation. DTs to split the data in a tree-like fashion until the algorithm arrives at a single answer (when the tree is fully grown). These decisions are calculated by optimising over the impurity at the proposed split using Equation (10):

(10) \begin{equation} G(Q_m, \theta) = \frac{n_{left}}{n_m} H(Q_{left}(\theta)) + \frac{n_{right}}{n_m} H(Q_{right}(\theta)) ,\end{equation}

where $Q_m$ is the data at node m, $\theta$ is a subset of data, $n_m$ is the number of objects at node m, $n_{left}$ and $n_{right}$ are the numbers of objects on the left and right sides of the split, $Q_{left}$ and $Q_{right}$ are the objects on the left and right sides of the split, and the H function is an impurity function that differs between classification and regression. For Regression, the Mean Square Error is used (defined in Equation (6)), whereas Classification often uses the Gini Impurity (defined in Equation (11)).

(11) \begin{equation} H(X_m) = \sum_{k \in J} p_{mk} (1 - p_{mk}),\end{equation}

where $p_{mk}$ is the proportion of split m that are class k from the set of classes J, defined formally in Equation (12):

(12) \begin{equation} p_{mk} = \frac{1}{n_m} \sum_{y \in Q_m} [\![ y = k]\!] ,\end{equation}

where $[\![ x ]\!]$ is the indicator function identifying the correct classifications.

3.5. ANNz2

The ANNz2Footnote e software (Sadeh et al. Reference Sadeh, Abdalla and Lahav2016) is is another ensemble method, combining the results of many (in this work we use 100) randomly assigned machine learning models as a weighted average from the pool of NNs and boosted decision trees, using settings noted in Bilicki et al. (Reference Bilicki2018, Reference Bilicki2021). However, whereas Bilicki et al. (Reference Bilicki2021) use the ANNz functionality to weight the training set by the test set feature distributions, here we do not use this option for two reasons. First, this work is designed for larger surveys to be completed, for which we do not know the distributions, so we are unable to effectively weight the training samples towards future samples. Second, when attempted, the final outputs were not significantly different, whether the training sets were weighted or not.

3.6. GMM+GPz

The GPz algorithm is based upon Gaussian Process Regression, a ML algorithm that takes a slightly different track than traditional methods. Whereas most algorithms model an output variable from a set of input features using a single, deterministic function, GP use a series of Gaussians to model a probability density function to map the input features to output variable. The GP algorithm is extended further in the GPz algorithm to handle missing and noisy input data, through the use of sparse GPs and additional basis functions modelling the missing data (Almosallam, Jarvis, & Roberts Reference Almosallam, Jarvis and Roberts2016a; Almosallam et al. Reference Almosallam, Lindsay, Jarvis and Roberts2016b).

Following Duncan (Reference Duncan2022), we first segment the data into separate clusters using a Gaussian Mixture Model before training a GPz model (without providing the redshift to the GMM algorithm) on each cluster, the idea being that if the training data better reflects the test data, a better redshift estimate can be made. We emphasise that no redshift information has been provided to the GMM algorithm, and the clusters determined by the algorithm is solely based on the g – z, and W1-W4 optical and infrared photometry—the same photometry used for the estimation of redshift.

The GMM uses the Expectation-Maximization (EM) algorithm to optimise the centers of each cluster it defines by using an iterative approach, adjusting the parameters of the models being learned in order to maximise the likelihood of the data belonging to the clusters assigned. The EM algorithm does not optimise the number of clusters, which must be balanced between multiple competing interests:

• • The size of the data—the greater the number of clusters, the more chance the GMM will end up with insufficient source counts in a cluster to adequately train a redshift estimator

• • The number of distinct source types within the data—If the number of clusters is too small, there will be too few groupings that adequately split the data up into it is latent structure, whereas if it is too high, the GMM will begin splitting coherent clusters

This means that the number of components used by the GMM to model the data is a hyper-parameter to be fine-tuned. Ideally, the number of components chosen should be physically motivated—the number of classes of galaxy we would expect to be within the dataset would be an ideal number, so the ML model is only training on sources of the same type to remove another source of possible error. However, this is not necessarily a good option, as, due to the unsupervised nature of the GMM, we are not providing class labels to the GMM, and hence cannot be sure that the GMM is splitting the data into the clusters we expect. On the other hand, being unsupervised means the GMM is finding its own physically motivated clusters which do not require the additional—often human derived—labels. The lack of labels can be a positive, as human-decided labels may be based less on the actual source properties, and more on a particular science case (see Rudnick Reference Rudnick2021 for further discussion).

In this work, we optimise the number of components hyper-parameter, where the number of components n_comp is drawn from $ \mathrm{n_{comp}} \in \{1, 2, 3, 5, 10, 15, 20, 25, 30\}$ . We emphasise that the number of components chosen is not related to the number of redshift bins used for classification and has an entirely separate purpose. The primary metric being optimised is the Bayesian Information Criteria (BIC, Schwarz Reference Schwarz1978):

(13) \begin{equation} \operatorname{BIC} = -2\log\!(\hat{L}) + \log\!(N)d,\end{equation}

where $\log\!(\hat{L})$ is the log likelihood of seeing a single point drawn from a Gaussian Mixture Model, defined in Equation (14), and d is the number of parameters.

(14) \begin{equation}\log{\hat{L}} = {\sum_i \log{\sum_j \lambda_j \mathcal{N}(y_i|\unicode{x03BC}_j, \sigma_j)}}\end{equation}

where $y_i$ is an individual observation, $\mathcal{N}$ is the normal density with parameters $\sigma_j^2$ and $\unicode{x03BC}_j$ (the sample variance and mean of a single Gaussian component), and $\lambda_j$ is the mixture parameter, drawn from the mixture model.

The BIC operates as a weighted likelihood function, penalising higher numbers of parameters. The lower the BIC, the better.

Figure 12.

Optimisation of the Bayesian Information Criteria (top), Test Size (middle), and Outlier Rate (bottom) across a range of components.

Figure 13.

Redshift distribution of the 30 GMM components. In each case the vertical axis shows the count.

Fig. 12 shows the BIC (Equation (13); top panel) with error bars denoting the standard error of each component, the average test size of each component with the error bars denoting the minimum and maximum test set size for each component (middle), and the photometric error in the form of both the outlier rate (Equation (1)), and the accuracy (Equation (7)), with error bars denoting the standard error (bottom).

Table 3.

Regression results table comparing the different algorithms across the different error metrics (listed in the table footnotes). The best values for each error metric are highlighted in bold.

^aCatastrophic outlier rate, Equation (1).

^b2$\sigma$ outlier rate, Equation (3).

^cResidual Standard Deviation, Equation (4).

^dNormalised Median Absolute Deviation, Equation (5).

^eMean Square Error, Equation, Equation (6).

Fig. 12 shows that while the n_comp is being optimised for lowest BIC, this has the additional benefit of lowering the resulting redshift estimation error (Fig. 12; middle and bottom panels)—showing that the clusters being identified by the GMM algorithm are meaningful in the following redshift estimation. A value of 30 is chosen for the n_comp, despite the BIC continuing to decline beyond this point. However, the number of sources in the smaller clusters defined by the GMM becomes too small to adequately train a GPz model.

Once the data are segmented into 30 components (an example from one random seed is shown in Fig. 13), a GPzFootnote f model is trained for each component. The GPz algorithm is based around sparse GPs, which attempt to model the feature space provided using Gaussian components.

3.7. Training method

ML algorithms are typically set up and trained following one of two procedures:

1. 1. Training/Validation/Test Splits

   • • The data is split into training, testing, and validation sets (for this work, the data are split into 50%/20%/30% subsets). The ML algorithm is trained on the training set, with model hyperparameters optimised for the validation set. Once optimised, the test set is used to estimate the model’s generalisatibility.

   • • This method is utilised by the ANNz and GMM algorithms

2. 2. k-Fold Cross-Validation

   • • The dataset is split into two sets (for this work, the data are split into 70%/30% subsets), used as training and test sets. Differing from the first method, this method trains and optimises the ML algorithms on the training set alone, before testing the optimised models on the test set.

   • • The training set is split into $k_f$ subsets. $k_f$ models are trained on $k_f - 1$ subsets, and hyperparameters optimised and validated against the remaining subset.

   • • In this work, we use a value of 5 for $k_f$ , with the k-fold Cross-Validation algorithms used to optimise the hyperparameters of the kNN and RF algorithms.

The externally developed software (ANNz and GPz) both operate using training/validation/test split datasets. This is preferable for large, mostly uniform distributions, as it greatly reduces training time. However, for highly non-uniform distributions, the under-represented values are less likely to be involved in all stages of training, validation and testing. Hence, for the kNN and RF algorithms, where we control the training process, we choose the k-fold cross-validation method of training and optimising hyperparameters, to best allow the under-represented high-redshift sources to be present at all stages of training.

3.7.1. Photometry used in training

All algorithms use the same primary photometry—g, r, i, z optical magnitudes and W1, W2, W3, and W4 infrared magnitudes. However, the different algorithms vary in how they treat the uncertainties associated with the photometry. For the simple ML algorithms (kNN and RF), the uncertainties are ignored. ANNz computes their own uncertainties using a method based on the kNN algorithm, outlined in Oyaizu et al. (Reference Oyaizu2008), and GPz uses them directly in the fitting of the Gaussian Process.

Table 4.

Classification results table comparing the different algorithms across the different error metrics (listed in the table footnotes). The best values for each error metric are highlighted in bold.

^aCatastrophic outlier rate, Equation (1).

^b2$\sigma$ outlier rate, Equation (3).

^cResidual Standard Deviation, Equation (4).

^dNormalised Median Absolute Deviation, Equation (5).

^eAccuracy, Equation (7).