The PAU Survey: narrowband photometric redshifts using Gaussian processes

PAUS is a narrowband photometric galaxy survey aimed at mapping the large-scale structure of the Universe up to $i\sim 23.0$. Using $40$ narrow bands spaced by $100$ Å in the range between $4500$ to $8500$ Å (filter responses visualised in Eriksen et al., 2019, and Fig. 4), PAUS aims to achieve redshifts with a precision of $\sigma_{\textrm{RMS}}<0.0035(1+z)$ for galaxies with $i_{\textrm{auto}}<22.5$. PAUS uses the PAUCam instrument (Padilla et al., 2019) on the $4$ m William Herschel Telescope (WHT) at Observatorio del Roque de los Muchachos (ORM) in La Palma. It has observed more than $50$ deg² of sky since the beginning of $2016$, and observations to full depth in all narrow bands for $100$ deg² are planned.

The PAUS forced-aperture coadded photometry has its aperture defined by using the $50$ per cent light radius ($r_{50}$), the point spread function (PSF), ellipticity and Sérsic index of COSMOS morphology, such that the fluxes measure a fixed fraction of light. The reader is referred to Eriksen et al. (2019) for detailed information on how the PAUS fluxes are measured. In this work we used the early data release from PAUS (objects are observed at least five times, using an elliptical aperture with $62.5$ per cent light radius), and select objects with $i_{\textrm{auto}}\leq 22.5$, entries with no missing measurement, and the COSMOS flag TYPE=0 (extended objects).

The Cosmic Evolution Survey (COSMOS, Scoville et al., 2007) covers a sky area of $2$ deg² ($149.47^{\circ}\leq\alpha\leq 150.7^{\circ}$, $1.62^{\circ}\leq\delta\leq 2.83^{\circ}$) and is known for its high sensitivity, depth and an exceptionally low and uniform Galactic extinction ($E_{B\text{-}V}\sim 0.02$).

In this work we used photometry from the COSMOS$2015$ Catalogue (Laigle et al., 2016); it is a highly complete mass-selected sample to very high redshifts, highly optimised for the study of galaxy evolution and environments in the early Universe. The COSMOS$2015$ Catalogue provides $30$ band photometry ranging from near UV to near infrared wavelengths, all these have been observed through multiple facilities, two of which are the Canada-Hawaii-France Telescope (CFHT) and Subaru Telescope (Miyazaki et al., 2002). From this catalogue we only use the CFHT $u^{*}$-band (Boulade et al., 2003) and Subaru $B$, $V$, $r$, $i^{+}$ and $z^{++}$ bands (Miyazaki et al., 2002), in conjunction with the narrowband photometry of PAUS. For simplicity, these bands will be referred to collectively as the $uBVriz$ bands; the superscripts are dropped for easier reading.

The zCOSMOS Survey (Lilly et al., 2007) targets galaxies in the COSMOS field using the Visible Multi-Object Spectrograph (VIMOS, Le Fèvre et al., 2003). zCOSMOS-Bright observed $20\,689$ galaxies in a sky area of $1.7$ deg², these galaxies have magnitudes $15<i_{\textrm{auto}}<22.5$ and redshifts in the range of $0.1<z<1.2$, its spectral range is in the red (rest-frame wavelength $5550$ Å to $9650$ Å) to follow strong spectral features around the $4000$ Å break to as high redshifts as possible.

In this work we use data from zCOSMOS-Bright DR$3$¹¹1http://www.eso.org/qi/catalog/show/65. Galaxies with redshift confidence class $3$ and $4$ (spectroscopic verification rate of $99$ and $99.8$ per cent, respectively) are selected and cross-matched with PAUS objects.

Using the aforementioned selection cuts, we cross-matched within $1^{\prime\prime}$ the $40$-narrowband photometry from PAUS, six broadband photometry ($uBVriz$) from COSMOS, and highly reliable redshifts from zCOSMOS to obtain a data sample of $8406$ galaxies, which is divided randomly into half for training and testing respectively. This sample uses a total of $46$ bands, and flux calibration between the broad and narrow bands is required as they are obtained from different surveys with different flux measurements. The calibration between these fluxes will be discussed in Section 4.

The colour-magnitude diagram of this sample is shown in Fig. 1, in comparison with the COSMOS$2015$ sample (all objects with TYPE=0 and detected in $r$ and $i$). The slight incompleteness in $i$ magnitude is due to the selection effects in brightness of the spectroscopic redshifts available.

The sample size may seem small, but is sufficient for the GP to work, since the GP essentially creates $4000+$ flux-redshift ‘templates’ to produce photo-$z$’s for objects in the testing set. However, we note that such a small training size has a major effect on the results of annz$2$ as this training size is close to the lower limit threshold suggested by Bonfield et al. (2010). We also note that the sample we have chosen is very similar to that of Eriksen et al. (2019), the only difference being that they have a more relaxed cut in the number of bands (N_BANDS), being 35<N_BANDS<40 (workable for a template code like bcnz$2$), while we used N_BANDS=40²²2The relaxed cut resulted in Eriksen et al. (2019) having a larger sample size of $10\,801$ objects.. When comparing results between delight and bcnz$2$, we will only compare photo-$z$’s of the exact same objects. Note that we have used the same broad bands as used by Eriksen et al. (2019).

delight ³³3https://github.com/ixkael/Delight (Leistedt & Hogg, 2017) is a hybrid template-based and machine learning photo-$z$ algorithm, which was constructed to combine the advantages, and minimise the disadvantages, of both types of algorithms. delight constructs a large collection of latent SED templates (or physical flux-redshift models) from training data, with a template SED library as a guide to the learning of the model. This conceptually novel approach uses Gaussian processes (GPs) operating in flux-redshift space. delight was featured in the results of the LSST Photo-$z$ Data Challenge $1$ (Schmidt et al., 2020), where it was found to have a low photo-$z$ bias but slightly broader PDFs.

A GP is a supervised learning method, which finds a distribution over the possible functions $f(x)$ that are consistent with the observed data $x$. Consider Fig. 2: suppose we have a set of observed variables $y=f(x)$, we can fit it using a GP, denoted as $f\sim\mathcal{GP}\left(\mu,k\right)$ , which assumes that the probability of all $f(x)$ is jointly Gaussian and representable by a mean function $\mu(x)$ and a covariance matrix $\Sigma(x)=k(x_{i},x_{j})$. $k(x_{i},x_{j})$ is the kernel function, which relates one variable $x_{i}$ to another $x_{j}$. An example case would be $\mu\equiv 0$ and a kernel function that takes the form of a squared exponential,

where $\sigma_{f}^{2}$ is the maximum allowable covariance between data (set by the errors on the observation), and $l$ is the tunable correlation length that determines the smoothness of the GP. In this simplistic case, the GP will try to find a marginalisation of all possible functions, but $\mu$ and $k$ can be modified if an underlying model of the data we want to fit is known. The covariance function is defined such that a smooth function is to be predicted.

Assuming that we have a set of training data $\{x_{i},f(x_{i})\}$ and would like to find the prediction $\{x_{*},f_{*}(x_{*})\}$, the GP models $f$ and $f_{*}$ as jointly Gaussian, $\mathcal{N}(\mu,\Sigma)$, and therefore

where $\Sigma=k(x_{i},x_{j})$ is the covariance between the training data, $\Sigma_{*}=k(x_{*},x_{i})$ the covariance between training and the predicted data (superscript $T$ denotes the transpose of the matrix), while $\Sigma_{**}=k(x_{*},x_{*})$ is the variance of the predicted data.

It follows from the above that the posterior $p(f_{*}|x_{*},x_{i},f_{i})$ is also Gaussian, therefore a predicted point $f_{*}(x_{*})$ is plotted (green dot in Fig. 2) is modelled by a Gaussian function (smooth blue line) which runs across all points, with its $95\%$ confidence interval ($\pm 1.96\sigma_{f_{*}}$) represented by the navy shaded area.

In the context of delight, GPs are used to calculate the predicted fluxes $\hat{F}$ at a certain redshift $z$ for a training object $i$ with fluxes $F_{i}$ and redshift $z_{i}$. This could be better understood by first defining the posterior photo-$z$ distribution $p(z|\hat{F})$ of an object in the testing set. For machine learning methods, it has the form

where $p(\hat{F}|z,z_{i},F_{i})$ is the prediction for fluxes of the training galaxy at a different redshift $z$, while $p(z|z_{i},F_{i})$ and $p(z_{i},F_{i})$ are the priors that provide the redshift distributions and abundances, generated from the training data, which are multiplied to give the combined probability $p(z,z_{i},F_{i})$ for a given redshift $z$ and training object with redshift $z_{i}$ and fluxes $F_{i}$. This is analogous to the one derived from template-based methods,

where $t_{i}$ is the template, $p(z|t_{i})p(t_{i})=p(z,t_{i})$ is the prior and $p(\hat{F}|z,t_{i})$ is the probability of the predicted flux $\hat{F}$ at redshift $z$ and for template $t_{i}$. Both equations are easily differentiated by the fact that for template-based methods, $p(z|\hat{F})$ is derived using a list of templates $t_{i}$, while for machine learning methods it is derived using the individual training set objects with fluxes $F_{i}$ and spectroscopic redshift $z_{i}$.

delight differs a little from the usual machine learning method in the sense that instead of finding a direct empirical relationship between the fluxes and redshifts of the training objects, it uses a GP to model the predicted fluxes of a training galaxy at different redshifts with the help of SED templates. This creates a latent flux-redshift template for each training object, where for a given set of fluxes in the testing set, it could be compared to several training templates to find the best predicted redshift.

The algorithm first fits a best-fit SED template to a particular training object $i$ with redshift $z_{i}$ and fluxes $F_{i}$ (multiple bands); the best-fit SED template is then used to formulate the mean function and kernel of a GP to build a flux-redshift template which could predict the expected fluxes of certain band filters when this object is redshifted to a different $z$. With each training object now becoming a flux-redshift template, the final photo-$z$ posterior distribution of a testing set object is determined by making a pairwise comparison of every training-testing pair, and a weighted solution is obtained based on the best fits of each pair.

In other words, we are computing the probability that the target galaxy has the same SED as the training galaxy but at a different redshift. delight is thus a hybrid template-machine-learning photo-$z$ algorithm in the sense that SED templates are used to ‘guide’ the creation of flux-redshift templates based on the training objects, or, if seen from another perspective, the GP ‘corrects’ the SED templates by using training data. We refer the reader to Leistedt & Hogg (2017) for more on Gaussian processes, and also for the full expressions of the $\mu$ and $k$ in relation to the filter responses, flux normalisations, linear mixtures of physical SED templates, and the manually configurable SED residual function of emission lines.

delight is advantageous over many other photo-$z$ algorithms as its output is less dependent on representative training data, and it does not strictly require the training set to use the same photometric bands. However, it still requires accurate spectroscopic redshifts, high quality training fluxes and representative templates to produce high quality photo-$z$ probability density functions (PDFs), or $p(z)$. As such, given a few photometric bands, delight is able to predict missing bands or fluxes in an entirely different set of photometric bands, and this function is utilised in Section 4.1 to predict and calibrate the flux values between two surveys.

The optimisation settings of delight used in this work are as follows. For the GP setup, the number of Gaussians to fit the filter curves (numGpCoeff) was set to $7$ instead of the default $20$, appropriately selected to accommodate the smaller full width half maximum (FWHM) of the narrowband filters. Other than that, we have mainly used the default hyperparameter settings for delight with the exception of the widths of the luminosity and redshift priors $\sigma_{\ell}$ and $\sigma_{z}$ (ellPriorSigma and zPriorSigma, see Leistedt & Hogg, 2017), which have been lowered to $0.2$ and $0.1$ respectively as they produced better results.

As mentioned earlier, the mean function and the kernel of the GP are modelled after the choice of emission lines and SED template sets. We replaced the $3$ default emission lines in delight with the list provided by Eriksen et al. (2019), although we note that the change in result for this is insignificant. As for the templates, we used the Brown et al. (2014) high-quality templates, which consist of $129$ SEDs derived from real nearby galaxies. These templates have wavelengths covering the ultraviolet to mid-infrared, and encompass a broad range of galaxy types including ellipticals, spirals, merging galaxies, blue compact dwarfs and luminous infrared galaxies. In this work we have also tested the performance of various other template sets (Coleman et al., 1980; Kinney et al., 1996; Bruzual & Charlot, 2003; Ilbert et al., 2006; Polletta et al., 2007); however they do not perform as well as those of Brown et al. (2014): the root-mean-square photo-$z$ errors could range between $21$ to $112$ per cent higher when these templates are used. Therefore, the results from these tests are not shown in this work.

We note that delight requires all magnitudes $m_{i}$ and magnitude errors to be converted into fluxes $F_{i}$ and flux variances, with a zero-point adjustment of $26.4$ in magnitude (i.e. $F_{i}=10^{-0.4(m_{i}-26.4)}$). We have also added a $3$ and $6$ per cent flux error in quadrature to the flux variances for the narrow and broad bands, respectively, to account for other flux errors from both the data and the model (values estimated via trial and error). It is also worth mentioning that while delight is capable of processing negative fluxes (non-detections), the reference band (referenceBand) used for flux normalisation only handles fluxes with positive values. In this work we have selected the narrow band $nb625$ as the reference band, or the COSMOS $r$-band in cases where narrow bands were not used.

Throughout this work, we use $z_{\textrm{map}}$ (the maximum a posteriori of the PDF) to represent the best point estimate photo-$z$ produced by delight. The output photo-$z$ PDF bins were set to be linear instead of logarithmic, with a stepsize of $0.001$, and a range of $0.02<z<1.65$, keeping close to the limits of the spectroscopic redshifts.

We are also interested in how delight compares to other common template-based or machine-learning-based methods besides bcnz$2$ and deepz. Therefore two other photo-$z$ algorithms, annz$2$ and bpz are also used in this work, using the same training and template sets, to be compared with the performance of delight. In the following paragraphs we briefly introduce the two algorithms and their optimisation settings.

annz$2$ ⁴⁴4https://github.com/IftachSadeh/ANNZ (Sadeh et al., 2016) is a machine-learning-based photo-$z$ algorithm which has been widely used in recent works (Bonnett et al., 2016; Jouvel et al., 2017; Bilicki et al., 2018; Soo et al., 2018; Schmidt et al., 2020) due to its high customisability and its ability to produce PDFs. It uses the Toolkit for Multivariate Data Analysis (TMVA, Hoecker et al., 2007) with root (Brun & Rademakers, 1997), which allows it to run multiple different machine learning algorithms for training, and outputs photo-$z$’s based on a weighted average of their performance. In this work we ran annz$2$ with a mixture of $3$ machine learning methods, namely artificial neural networks (ANNs), boosted decision trees (BDTs) and $k$-nearest neighbours (KNNs), see Hoecker et al. (2007) for detailed descriptions of these machine learning algorithms. An architecture of $N$:$\frac{2N+1}{3}$:$\frac{N+2}{3}$:$1$ was used for the ANN; the bagging method was used to boost the decision trees; a polynomial kernel was used for the KNN; while the other hyperparameters for each method were individually optimised for best performance. annz$2$ version $2.3.1$ was used in this work, and the mean value of the PDF, $z_{\textrm{pdf}}$ was chosen to represent the photo-$z$ point estimate.

bpz ⁵⁵5http://www.stsci.edu/~dcoe/BPZ/ (Benítez, 2000), on the other hand, is one of the longest-standing template-based photo-$z$ algorithms, and still widely used today (Martí et al., 2014; Bundy et al., 2015; Cavuoti et al., 2017; Tanaka et al., 2018; Joudaki et al., 2020; Raihan et al., 2020). Other than sharing the usual attributes of a template-based code, bpz uses Bayesian inference, prior information of redshift distributions and template interpolation to improve photo-$z$ results. bpz version $1.99.3$ was used in this work, and similar to delight the Brown templates were used, with the interpolation parameter set to $2$. We assumed the same functional form for the Bayesian priors as those used by COSMOS (Laigle et al., 2016). The peak of the PDF, $z_{\textrm{b}}$, was used as the best photo-$z$ point estimate.

Other than annz$2$ and bpz, the results of delight are also compared to the results of bcnz$2$, which was developed specifically for the PAUS data (Eriksen et al., 2019). bcnz$2$ is able to compute a linear combination of SED templates and is designed to deal with emission lines, extinction, and adjust zero-points between narrow and broad bands, all of which are crucial in the context of PAUS. The introduction of the code bcnz$2$ and its early demonstration of PAUS photo-$z$ can be found in Eriksen et al. (2019).

In the first step we define a parameter $R_{\textrm{g}}$, a correction factor estimated for each galaxy to be multiplied with all of its six $uBVriz$ broadband fluxes. Ideally, this factor is estimated by first finding the best-fit Brown template for each galaxy using only $40$ narrowband fluxes from PAUS and its true redshift. The best-fit template is then used to generate the predicted $uBVriz$ fluxes, and a weighted mean of the ratios between the predicted flux and the original COSMOS flux $R_{\textrm{g},b}$ is calculated for each band $b$, given by

where the sum is over the six COSMOS broad bands, and $\sigma_{R_{\textrm{g}},b}^{2}$ is the variance of $R_{\textrm{g},b}$. Here we have assumed that the Brown templates are sufficiently representative, and therefore the predicted flux derived from it is the true flux of the broad bands. We have also assumed that $R_{\textrm{g},b}$ should be almost the same across each band for each galaxy. This calibration is motivated by the fact that each galaxy requires a calibration between fixed-size and adaptive aperture photometry dependent on its apparent size.

We now explore three different methods to determine $R_{\textrm{g}}$ from the photometric data only.

After calibrating the COSMOS broadband fluxes for each galaxy, we proceed to calibrate the broadband magnitude offsets within each band. We perform a weighted least-squares fit between the predicted broadband fluxes (produced by delight using $40$ PAUS narrowband fluxes, the respective best fit Brown templates and zCOSMOS spec-$z$’s) and the original COSMOS $uBVriz$ fluxes in the training set, by using a simple linear equation,

where $F_{\textrm{p},b}$ is the predicted flux for band $b$, $F_{\textrm{g},b}$ the COSMOS broadband flux after undergoing the per-galaxy calibration, and $a_{b}$ and $c_{b}$ are constants to be optimised. The values of $a_{b}$ and $c_{b}$ estimated for each band using the training set are now used to calibrate the fluxes in the testing set, and these values are tabulated in Table 1. A weighted fit was implemented, with the inverse variances of the fluxes used as the weights, since we expect that objects which are brighter to have relatively lower variances, and by accounting for the variances of objects the fainter objects would be upweighted.

As expected from the table, the values of $a_{b}$ and $c_{b}$ are very close to $1$ and $0$ respectively, since the calibrated flux for aperture correction $F_{\textrm{g},b}$ is already very close to the predicted flux $F_{\textrm{p},b}$. Essentially, this process ‘straightens’ the correlation line, providing minor yet essential improvements to the overall calibration.

Figure 5 shows the correlation between the broadband fluxes predicted by delight (using spectroscopic redshifts, PAUS $40$ narrow bands and Brown templates) and the COSMOS broadband fluxes for our training set, both before and after calibration (red and blue, respectively). The figure only shows the result of the flux calibration method, as the other two methods look very similar graphically (which translates to a small difference in photo-$z$ results shown later in Section 5).

The RMS values displayed in Figure 5 show that for all bands, the scatter between the original fluxes with respect to the predicted fluxes has reduced by $63$ to $88$ per cent after the two-step calibration was done. The scatter at low fluxes for the $u$ and $B$-bands remains evident, which originated from the high uncertainty in flux measurements. Despite the large decrease in scatter, we note that the RMS value here is not a metric of improvement for calibration as we do not have the true values of the broadband fluxes in the matched apertures. However, the calibration of the broadband fluxes did translate into an improvement in photo-$z$ scatter and $68$th percentile error by about $70$ to $80$ per cent, as shown in Section 5.

In this work we use three metrics to quantify the performance of the photo-$z$ point estimates: the root-mean-square error ($\sigma_{\textrm{RMS}}$), the $68$th percentile error ($\sigma_{68}$) and the outlier fraction rate ($\eta_{\textrm{out}}$). With $\Delta z\equiv\frac{z_{\textrm{phot}}-z_{\textrm{spec}}}{1+z_{\textrm{spec}}}$, the above metrics are defined as follows:

Here $N$ is the total number of galaxies, while $Q$ is a percentile of the distribution. Since $\sigma_{\textrm{RMS}}$ is calculated without the outliers removed, it measures the overall scatter of the sample, whereas $\sigma_{68}$ measures the scatter with reduced sensitivity to outliers.

With similar motivations as Martí et al. (2014) and Eriksen et al. (2019), we hope to achieve an overall photo-$z$ error of $\sigma_{68}\leq 0.0035(1+z_{\textrm{spec}})$ for at least $50$ per cent of the testing sample after applying an appropriately chosen quality cut. We use the Bayesian ODDS ($\Theta$) parameter (Benítez, 2000) in delight, similar to its implementation in annz$2$ by Soo et al. (2018). $\Theta$ can be estimated from the photo-$z$ PDF, $p(z)$ using the equation

where $z_{\textrm{p}}$ is the peak of $p(z)$ and $k=0.01$. $\Theta$ ranges between $0$ and $1$, the higher the value the lower the $p(z)$ width, which implies a more precisely predicted photo-$z$ (though not necessarily accurate). The value of $k$ is arbitrary, appropriately selected such that not too many objects end up having $\Theta=1$. Therefore, an $x$ per cent quality cut on the sample keeps the top $x$ percent of objects with the highest values of $\Theta$.

To assess the quality of the $p(z)$, we use probability integral transform (PIT) plots and the continuous ranked probability score (CRPS). The PIT is the cumulative distribution function (CDF) at $z_{\textrm{spec}}$ while asserting the $p(z)$ to have an area of unity. Since the photo-$z$ CDF is $C(z)=\int^{z}_{0}p(z^{\prime})\ dz^{\prime}$, PIT is defined to be

A PIT distribution tells us on average if the $p(z)$ produced are ‘adequately shaped’: the shape of the PIT distribution can tell us if the $p(z)$ produced are generally too wide/narrow, or if the $p(z)$ are over/under-predicting the true redshift.

The CRPS on the other hand tells us how well the $p(z)$ encapsulates or predicts the true redshift ($z_{\textrm{spec}}$). The CRPS of a $p(z)$ can be expressed as

where $\mathcal{H}(z-z_{\textrm{spec}})$ is the Heaviside step function with

In this work, we use the symbol $\rho_{\textrm{CRPS}}$ to represent the average CRPS value of all galaxies in the testing sample, in which the smaller the value, the better the $p(z)$ are at predicting their true redshifts. We refer the reader to Polsterer et al. (2016) for a detailed description of both PIT and CRPS.

Finally, we also assess the quality of the redshift distribution $n(z)$. We can find how similar the spec-$z$ distribution $n_{\textrm{spec}}(z)$ is compared to the photo-$z$ distribution $n_{\textrm{phot}}(z)$ by estimating $n_{\textrm{RMS}}$, the root-mean-square difference between the distributions:

$n_{\textrm{RMS}}$ provides us a quantitative measure to compare the performances of photo-$z$ with distributions produced by different codes.

Rows $1$ and $2$ from Table 2 shows the photo-$z$’s produced when only trained using the broad and narrow bands individually, and we find that by combining both broad and narrow bands (rows $4$ to $6$), we have achieved at least $34$ and $20$ per cent improvement in the photo-$z$ scatter and $\sigma_{68}$, respectively (visualised in Fig. 6).

Rows $3$ to $7$ proceed to show the metrics for each calibration method, and on average, the performance of each method is quite similar, all within $4$ to $16$ per cent difference in $\sigma_{\textrm{RMS}}$ and $\sigma_{68}$, respectively. Statistically, the flux calibration method seems to perform slightly better compared to the remaining ones, with the exception of the photo-$z$ calibration method having better values of $\eta_{\textrm{out}}$ and $n_{\textrm{RMS}}$. This suggests that while the photo-$z$’s produced by training with only $40$ narrow bands are not as competitive as when trained with all $46$ bands and calibrated broad bands (see Table 2 and Fig. 6), it is however sufficient to guide the calibration process. Note that we have also included the results of delight run as a pure template code when calibrated using the flux calibration method for comparison, and we see that without the help of the GP, the photo-$z$ results are similar for most metrics except a degradation in scatter of up to $-33.5$ per cent. Therefore the good results of delight shown here are mainly due to the use of the Brown templates, the flux calibration, the combination of broad and narrow bands, and also the work of the GP.

As the three calibration methods presented in Section 4 all result in very similar photo-$z$ performance, we will only show results for the flux calibration method in the following. It is notable however that in all cases, the photo-$z$ requirement of $\sigma_{68}<0.0035(1+z)$ is achievable for all objects at $i_{\textrm{auto}}<20.0$, or objects with a $40$ per cent $\Theta$ cut at $i_{\textrm{auto}}<22.5$. All three methods also shows that despite such high percentage $\Theta$ cuts being implemented, a significant number of high photo-$z$ objects still remain in the sample.

Since the delight results for each of the three calibration methods are very similar to each other, we decided to select only the flux calibration method to be compared to the results obtained by the two other algorithms used in this work, annz$2$ and bpz. We also include the point estimates from Eriksen et al. (2019). The values of $\sigma_{\textrm{RMS}}$, $\sigma_{68}$ and other relevant metrics obtained from these algorithms are shown in rows $8$ to $11$ of Table 2, and visualised in Fig. 7.

From the figure, it is found that annz$2$, being a purely machine-learning based algorithm, is underperforming compared to the other algorithms. This machine-learning method is unable to make full use of the extra information provided by the $40$ narrow bands, and is shown to perform better without them. This is partially due to the problem of the curse of dimensionality (Bellman, 1957), sharply diluting the pattern recognition power of the algorithm as the number of inputs increases. Besides, the very small training sample size may have heavily affected the potential of annz$2$. Here we note however that the deep learning code deepz is shown to work well on a similar sample (Eriksen et al., 2020), therefore we hope to do follow-up evaluations of annz$2$ on PAUS data in the future when a larger training set is available.

In terms of the quality of the point estimate photo-$z$’s, delight is shown to fare well against bcnz$2$ and bpz (Fig. 7), both of which are purely template-based methods. As both delight and bpz used the same template sets in this case (i.e. the Brown templates), we find that the Gaussian process contributed to $25$ and $9$ per cent improvement in the scatter and $\sigma_{68}$, respectively, as compared to the pure template fit of bpz.

Despite the similarities in the point estimates for the entire sample (Table 2), when we cut the sample in percentages of $\Theta$ (Figs. 8 and 9), we see two major differences. Firstly, the cut in $\Theta$ for bpz does not systematically remove objects with high uncertainties (especially for objects brighter than $i_{\textrm{auto}}=21$); and secondly, the cut in $\Theta$ for bpz selectively removes objects with lower photo-$z$. In both cases, delight is shown to not only perform better in this regard as compared to bpz, but also better than all other algorithms shown.

A selection of sample $p(z)$ produced by each algorithm is shown in Fig. 10, while the overall quality of the $p(z)$ produced are visualised in the PIT plots as shown in Fig. 11. Once again we see delight on average producing superior $p(z)$ compared to annz$2$ and bpz: it is obvious from the PIT plots that the $p(z)$ produced by annz$2$ are too narrow (a U-shaped distribution), while those by bpz are too wide (a significant central peak). In terms of $\rho_{\textrm{CRPS}}$ (see Table 2), delight once again performs better than both bpz and annz$2$, where the adequate shapes and accurately positioned peaks of the $p(z)$ provide good predictions of the true redshift.

We note that the $p(z)$ produced by annz$2$ are ragged compared to bpz and delight, this is due to the limited training sample size and the low number of network committees used. We intend to look into several methodologies to smoothen machine-learning based $p(z)$ which are limited by such conditions; this is left for future work. The limited testing size has also produced an $n_{\textrm{spec}}(z)$ distribution which is not smooth, thus despite annz$2$ producing an $n(z)$ closest to the spectroscopic distribution (lowest $n_{\textrm{RMS}}$), it may have experienced overfitting. Having said that, for the different delight runs shown in Table 2, the values of $n_{\textrm{RMS}}$ are consistent with the other metrics. Therefore, we leave the analysis of $n(z)$ to future work when a large enough testing sample is available.

As we compared the photo-$z$ results, we discovered that there are some galaxies that have similar delight and bcnz$2$ photo-$z$ values, however these redshift values are far from their respective zCOSMOS spectroscopic redshifts or broadband photo-$z$’s. Since both bcnz$2$ and delight utilise the PAUS narrow bands, we expect that the photo-$z$’s they produce are more sensitive to emission lines as compared to photo-$z$’s produced using only broad bands. Therefore, we suspect that objects that have similar photo-$z$ values for delight and bcnz$2$ but have disagreeing spec-$z$ values to be an indication of either having (1) a catastrophic zCOSMOS spectroscopic redshift⁶⁶6While we have already selected to use only secure spectroscopic redshifts in this work, we still deem this as a possibility, since a $1$ per cent outlier rate in $4000+$ spec-$z$ measurements may still yield $40$ objects, which is within the same order of number of objects being investigated in this section. Our results later in this section however have verified that most of the outliers are not caused by catastrophic spectroscopic redshifts., (2) outlier broadband or narrowband fluxes, or (3) misidentification of close neighbours.

For the purpose of this inquiry, we have selected $30$ objects from the sample which are photo-$z$ outliers in $z_{\textrm{Delight}}$ vs. $z_{\textrm{spec}}$ or $z_{\textrm{BCNz}}$ vs. $z_{\textrm{spec}}$, yet are not outliers in $z_{\textrm{Delight}}$ vs. $z_{\textrm{BCNz}}$. Mathematically, they satisfy the following conditions:

1. 1.

   $\frac{\left|z_{\textrm{Delight}}-z_{\textrm{spec}}\right|}{1+z_{\textrm{spec}}}\geq 0.15$ or $\frac{\left|z_{\textrm{BCNz}}-z_{\textrm{spec}}\right|}{1+z_{\textrm{spec}}}\geq 0.15$, and

2. 2.

   $\frac{\left|z_{\textrm{Delight}}-z_{\textrm{BCNz}}\right|}{1+\frac{z_{\textrm{Delight}}+z_{\textrm{BCNz}}}{2}}<0.15$.

Note that the $z_{\textrm{Delight}}$ used here refers to the photo-$z$ produced using the flux calibration method, trained using $46$ bands guided by the Brown templates.

These $30$ objects are visualised in the redshift-redshift plots in Fig. 12. Note that in the following paragraphs, we will define a photo-$z$ to be catastrophic if it is found to be an outlier with respect to its spec-$z$, as defined mathematically above. These objects are found to have faint magnitudes ($i_{\textrm{auto}}>19.75$) and small angular sizes ($r_{50}<60$ ACS pixels, or $1.8^{\prime\prime}$), which describe most galaxies of interest for PAUS. We study several different attributes of these objects, namely their respective photo-$z$’s by delight, bcnz$2$ and lephare, photo-$z$ PDFs, best-fit templates (Brown and GP), spectra and images. We summarise important observations according to their respective attributes below.

Photo-$z$’s. While these $30$ objects have been identified as outliers when trained using $46$ bands, we find that two-thirds of these objects have non-catastrophic photo-$z$’s when trained with either only the broad or narrow bands, respectively. In other words, only one-third of these objects have catastrophic photo-$z$’s regardless of which bands were used in the training or fitting process. This suggests that most of the time, outlier fluxes in the broad or narrow bands may have caused a degradation in photo-$z$ quality when trained together (more on this in the templates paragraph below). We have also made a comparison between delight photo-$z$’s with those produced by lephare for the COSMOS2015 catalogue (Laigle et al., 2016), and found that in fact half of the $30$ objects have non-catastrophic lephare photo-$z$’s. This suggests that the infrared $yJHK$ bands could have played a role in improving the PAUS photo-$z$’s, and could be incorporated in future trainings in case the PAUS photometry is problematic⁷⁷7We note that these additional bands will not be available over most of PAUS, which targets Canada-France-Hawaii Telescope Legacy Survey (CFHTLS) wide fields W1 to W4. There is however some infrared data on these fields provided by the Wide-field InfraRed Camera (WIRCam) and the VISTA Kilo-degree Infrared Galaxy Survey (VIKING)..

Photo-$z$ PDFs. We inspected the secondary/tertiary peaks of the PDFs for all delight runs (trained with $6$ broad bands, $40$ narrow bands, or both), and find that less than $20$ per cent of these secondary/tertiary peaks coincide with their respective spec-$z$’s. We deduce that despite the importance of secondary PDF peaks in redshift distributions, they do not significantly influence the photo-$z$ quality of these $30$ objects.

Templates. delight utilises the $129$ Brown et al. (2014) templates and the $4203$ training objects to guide the GP to produce the same number of new flux-redshift templates, which are used to produce photo-$z$’s for the objects. In the training process, delight would always choose one best-fit Brown template for each training galaxy to be trained by the GP. Here we inspected two different kinds of best-fit Brown templates to these $30$ outliers: one fixed at the spec-$z$, and the other with the redshift as a free parameter. In both cases, we examined

1. 1.

   if the objects fit to the same templates when trained with only broad bands, only narrow bands, or both, respectively;

2. 2.

   if there are any trends in galaxy morphological types, based on the galaxy type classification indicated by the template;

3. 3.

   if there is any correlation between the $\chi^{2}$ value of the best-fit templates and the quality of photo-$z$’s; and

4. 4.

   if any outlier narrowband fluxes can be identified as the cause of the degradation of photo-$z$.

As expected, we find that $70$ per cent of the outlier objects have different best-fit Brown templates between the fits at fixed photo-$z$ and spec-$z$, which contrasts with the case for non-outliers at only $35$ per cent. We also find that only slightly more than a third of both the outlier and non-outlier objects were fitted to the same templates when trained using broad bands as compared to trained with all $46$ bands. The high percentage of objects with different template fits at different reference redshifts (photo-$z$ or spec-$z$) and flux combinations (broad bands, narrow bands, or both) also resulted in no trend in galaxy morphological types among the outliers.

However, it was found that up to $60$ per cent of the objects have their best-fit template $\chi^{2}$ value correlating with the quality in photo-$z$, which further affirms the usage of this as a metric to remove unreliable photo-$z$’s (see Section 6.2), as also attempted by Eriksen et al. (2019) and Eriksen et al. (2020).

Perhaps a more significant finding from the study of the best-fit templates is the ability to identify outlier narrowband fluxes. Fig. 13 shows an example which highlights the importance of identifying outlier narrowband fluxes, which is shown to significantly affect the photo-$z$ results. It was found that a third of the $30$ objects contained outlier narrowband fluxes, which results in entirely different template fits and photo-$z$’s when trained with narrow bands, as compared to when trained with broad bands only. Among these $10$ objects, $8$ of them are shown to have worse photo-$z$ as compared to training without the narrow bands. We find indications for a significant fraction of narrowband flux outliers also for galaxies without catastrophic redshift failures. Forthcoming PAUS data reductions will therefore implement methods to identify and correct flux outliers.

Images. We inspect the individual object images compiled by zCOSMOS DR$3$, these are $5^{\prime\prime}\times 5^{\prime\prime}$ images observed by the Hubble Space Telescope/Advanced Camera for Surveys (HST/ACS) in the F$814$W filter (Koekemoer et al., 2007). Among the $30$ outlier objects, we find $63.3$ and $26.7$ per cent of them having bright neighbours within $5^{\prime\prime}$ and $3^{\prime\prime}$ of the primary source, respectively. Having said that, we have not found any correlation between the presence of bright neighbours to the other attributes that we have studied thus far. In fact the opposite is true: we find that $60$ per cent of the objects with outlier narrowband fluxes actually have primary sources without any bright neighbours in vicinity.

Spectra. So far we have assumed that the zCOSMOS spectra obtained are reliable, as only entries with high-confidence quality flags have been selected for training (see Section 2.3). In order to probe further, we examined the one-dimensional spectra obtained by the VIMOS spectrograph, which is processed by the VIMOS Interactive Pipeline and Graphical Interface (VIPGI, Scodeggio et al., 2005) to produce the zCOSMOS spec-$z$’s used in this work. The spectra have a range between $5500$ Å and $9450$ Å, measured with a resolution of $R\sim 600$ at $2.5$ Å per pixel (Lilly et al., 2009).

We used the redshift measurement tool ez (Garilli et al., 2010) to inspect the spectra of the $30$ outlier objects, and compared our best fits to the spectroscopic redshift produced by zCOSMOS, and also the photo-$z$’s produced by delight, bcnz$2$, deepz, lephare (COSMOS2015) and those of Alarcon et al. (2021).

Upon inspection, we find that up to $10$ of these objects ($33$ per cent) have disputable zCOSMOS spec-$z$ (e.g. two possible redshift values, different best-fit redshift values, line confusion, and low signal-to-noise). However, most of these potential spec-$z$ failures could be forced-fitted to the zCOSMOS spec-$z$ and still look satisfactory, which leaves only $2$ ($6.7$ per cent) of these objects having truly catastrophic spec-$z$’s. Both these objects are found to have better ez fits at redshift values within $10$ per cent uncertainty from the photo-$z$’s produced by delight and other algorithms. The spectrum of one of these objects is shown in Fig. 14. We have also found one isolated case where the spectra belonged to a bright neighbour and has been mismatched to the PAUS photometry.

Generally, the higher-redshift objects are identified by clear O II ($3727.1$ Å) emission lines, while the lower-redshift objects are identified by clear H$\alpha$ ($6564.6$ Å) emission lines. We therefore conclude that although catastrophic spec-$z$’s played a role in this situation, our results did not provide enough evidence to say that it is a major cause for catastrophic photo-$z$’s produced by bcnz$2$ and delight. This is not surprising since we have only selected secure spectroscopic redshifts from COSMOS to be used in this work. However this highlights the usefulness of multiple PAUS photo-$z$’s being used to determine failure rates in insecure spectroscopic redshifts.

To summarise this part, we believe that the potentially important source for catastrophic photo-$z$’s in the context of PAUS are the outlier narrowband fluxes, with weak evidence for the existence of a small number of spec-$z$ failures. We leave the tackling of outlier narrowband fluxes to future work, but in the following section, we attempt to improve our process to identify and remove these outlier photo-$z$’s.

In Figs. 8 and 9 we have used the Bayesian odds ($\Theta$) to cut the sample, and the aim of this was to keep as many objects as possible while achieving the goal of $\sigma_{68}\leq 0.0035(1+z)$. Here, we extend our previous results further towards that goal by introducing several new metrics to better separate the photo-$z$ outliers from the sample. These metrics are motivated by the inspection of the $30$ outliers in Section 6.1, and they are defined as follows:

1. 1.

   The Delight-BCNz2 metric ($\Delta_{\textrm{DB}}$),

   $$\Delta_{\textrm{DB}}\equiv\frac{\left|z_{\textrm{Delight}}-z_{\textrm{BCNz}}\right|}{1+\frac{z_{\textrm{Delight}}+z_{\textrm{BCNz}}}{2}},$$

   (16)

   a metric used to identify the similarity between delight and bcnz$2$ photo-$z$’s. It is plausible that, in general, the closer the photo-$z$’s between the two algorithms, the more reliable they are;

2. 2.

   The Delight photo-$z$ standard deviation ($\sigma_{\textrm{D}}$), which is the standard deviation between all delight photo-$z$ runs regardless of calibration method and number of bands. Smaller deviations could indicate more reliable photo-$z$’s;

3. 3.

   The chi-squared value of the best-fit Brown template ($\chi^{2}_{\textrm{t}}$), where we identified a trend that the better the fit, the more reliable the photo-$z$; and

4. 4.

   The broadband-narrowband complementary metric ($\rho^{2}$),

   $$\rho^{2}\equiv\int p_{\textrm{BB}}(z)p_{\textrm{NB}}(z)\,dz,$$

   (17)

   where $p_{\textrm{BB}}(z)$ and $p_{\textrm{NB}}(z)$ are the $p(z)$ produced by delight when trained with only broad bands and only narrow bands, respectively. By multiplying these two $p(z)$ and summing over the distribution at each step $i$, we can identify the consistency between the broadband and narrowband $p(z)$. A higher value of $\rho^{2}$ means a larger overlap, which indicates more reliable photo-$z$’s.

Together with $\Theta$ and the delight photo-$z$ error ($\delta z$), we yield a total of $6$ metrics to experiment with. Using the results from the flux calibration method, we generate and test the individual performance for each of these metrics. For each metric, we measure the $\sigma_{\textrm{RMS}}$ and $\sigma_{68}$ after systematically removing objects with the worst metric values, $10$ per cent of the total sample size each time, until we reach a sample size of only $40$ per cent.

We also repeat the exercise by using combined cuts on several metrics, testing all $57$ combinations of the $6$ metrics. We note that we do not combine the metrics by averaging or multiplying them, as it would have diluted the impact of the individual metrics. Instead, we rank the values for each metric individually (from best to worst), and remove objects rank by rank, starting with metric values lying in the worst rank. E.g. for the combination of metrics $\Theta+\Delta_{\textrm{DB}}$, we first remove all objects which share the worst values of $\Theta$ and $\Delta_{\textrm{DB}}$, then remove all objects sharing the second worst values of them, and so on, until we reach a required sample size percentile ($90$, $80$, etc), where we output the values of $\sigma_{\textrm{RMS}}$ and $\sigma_{68}$. We visualise the performance of these metric cuts at several percentiles for $\sigma_{68}$ with respect to $i_{\textrm{auto}}$ (cumulative) in Fig. 15.

We find that each performance metric cuts the sample differently: while metric cuts of $\sigma_{\textrm{D}}$ and $\rho^{2}$ reduce the scatter ($\sigma_{\textrm{RMS}}$) significantly, metric cuts of $\Theta$ and $\Delta_{\textrm{DB}}$ reduce the $\sigma_{68}$ instead. The metric $\chi^{2}_{\textrm{t}}$, however does not seem to bring any significant improvement to the results. We have also plotted a cut in $\Delta z=\frac{\left|z_{\textrm{phot}}-z_{\textrm{spec}}\right|}{1+z_{\textrm{spec}}}$ (bottom-left panel in Fig. 15), which is the theoretical ‘best metric’, providing an upper limit to be compared with the performance of each of the metrics. Here we noticed that even with the theoretical best metric, a cut of slightly lesser than $70$ per cent (blue line) on the sample is still necessary to fulfil the PAUS target of $\sigma_{68}<0.0035(1+z)$ (dotted line) for delight.

Therefore, we select the $60$ per cent cut (navy line, retaining $60$ per cent of galaxies) as a benchmark to assess the performance of these metrics, we do so by locating where this line cuts the dotted line (i.e., finding the maximum value of $i_{\textrm{auto}}$ where the photo-$z$’s achieves the PAUS target at $60$ per cent cut). From Fig. 15, it is clear that cutting in all $6$ metrics does not necessarily outperform the performance when cutting with only $\Theta$, so we searched for the best combination of metrics for $\sigma_{\textrm{RMS}}$ and $\sigma_{68}$ separately.

For $\sigma_{\textrm{RMS}}$, the best combination of metrics is $\Delta_{\textrm{DB}}+\sigma_{\textrm{D}}+\rho^{2}$, and this combination achieves $\sigma_{\textrm{RMS}}<0.0035(1+z)$ at $i_{\textrm{auto}}<19.27$ at $60$ per cent cut, a significant improvement to the case when only $\Theta$ was used, where it did not cut the line at all. For $\sigma_{68}$, the best combination of metrics is $\Theta+\Delta_{\textrm{DB}}$ where it reached $\sigma_{\textrm{68}}<0.0035(1+z)$ at $i_{\textrm{auto}}<21.25$ at $60$ per cent cut, which is also a significant improvement as compared to $\Theta$ at $i_{\textrm{auto}}<20.88$. Here we note that in fact using $\Delta_{\textrm{DB}}$ alone, the target can be reached at a higher limit of $i_{\textrm{auto}}<21.50$, which highlights the significance of a synergy between delight and bcnz$2$ in selecting a high quality photo-$z$ sample.

Finally, we also show the performance of the metrics in terms of the completeness with respect to the photo-$z$ (using delight’s flux calibration method), visualised in Fig. 16. We find that metrics like $\sigma_{\textrm{D}}$ and $\rho^{2}$ tend to selectively remove high photo-$z$ objects, while $\Theta$, $\chi^{2}_{\textrm{t}}$ and $\Delta_{\textrm{DB}}$ tend to remove mid-ranged photo-$z$ objects. In general, a cut using all $6$ performance metrics at $60$ per cent cut shows a balanced result in the completeness, keeping a sufficient number of high redshift objects in the sample.

To summarise the performance of the individual metrics,

• •

  $\chi^{2}_{\textrm{t}}$ is the least-performing metric here; it does not bring significant positive impact to the results;

• •

  Cuts in $\sigma_{\textrm{D}}$ and $\rho^{2}$ help to improve the scatter, however, they tend to selectively remove higher photo-$z$ objects from the sample;

• •

  $\Theta$ and $\delta z$ show very similar results, however $\Theta$ tends to keep more high photo-$z$ objects in the sample; and

• •

  $\Delta_{\textrm{DB}}$ is the best-performing metric here, and we recommend the use of such a metric to remove outlier photo-$z$’s from a sample.