Estimating Galaxy Parameters with Self-organizing Maps and the Effect of Missing Data

Our analysis is based on two publicly available mock catalogs (described in Laigle et al. 2019) with observed-frame photometry extracted from the Horizon-AGN lightcone ⁵ (Dubois et al. 2014). Horizon-AGN is a cosmological hydrodynamic simulation that contains 1024³ dark matter particles in a box of L_box = 100 h⁻¹ Mpc, corresponding to a dark matter mass resolution of 8 × 10⁷ M_⊙. The simulation is run with the adaptive mesh refinement code ramses (Teyssier 2002), following the evolution of the gas, including the effects of gravity, hydrodynamics, gas cooling and heating processes, star formation and stellar feedback, and feedback from black holes (see Kaviraj et al. 2017, for details). Gas is heated from an uniform UV background after the reionization epoch, according to Haardt & Madau (1996), and it cools down to 10⁴ K via H, He, and metals, per Sutherland & Dopita (1993). Star formation is activated only in regions with gas number density n > 0.1 H cm⁻³, and according to the Schmidt law (Schmidt 1959), ${\dot{\rho }}_{\mathrm{SFR}}=\epsilon {\rho }_{\mathrm{gas}}/{\tau }_{{ff}}$, where ${\dot{\rho }}_{\mathrm{SFR}}$ is the SFR mass density, is the star formation efficiency, ρ_gas is the gas mass density, and τ_ff is the freefall time of the collapsing gas cloud. Stellar feedback is provided by stellar winds and type Ia and II supernovae, and AGN feedback is a combination of the radio or quasar mode based on the black hole accretion rate.

In Laigle et al. (2019), the 1 deg² lightcone is built by running the AdaptaHOP halo finder (Aubert et al. 2004) on the Horizon-AGN lightcone over 0 < z < 4 and selecting structures with a density threshold 178 times the average matter density at that redshift, and then considering the total stellar mass within each halo. The extracted lightcone contains 789,354 galaxies with stellar masses >10⁹ M_⊙ in a redshift range z = 0–4. For each galaxy, Laigle et al. (2019) compute its SED by adopting single stellar population models from Bruzual & Charlot (2003) and a Chabrier IMF (Chabrier 2003). They adopt an empirical conversion relation from the gas-phase metals to dust and estimate the spatial distribution of that dust from the gas density in each cell. Integrated along the line of sight, this provides the column density of dust. Dust attenuation is then applied to the spectra by assuming the dust follows the R_V = 3.1 Milky Way dust model of Weingartner & Draine (2001). They do not take into account scattering in or out of the line of sight, but only dust absorption. This results, for example, in higher (by ≈0.8dex) estimated magnitudes in the UV (Laigle et al. 2019). They also account for absorption by the IGM, but not for foreground extinction by the Milky Way, which is expected to be corrected for in any observed survey that would be compared with this simulated data set.

Apparent magnitudes are obtained by convolving the thus-computed galaxy spectra with the desired filter profiles. Of the two mock catalogs ⁶ we use, one is the COSMOS-like mock where the photometry matches the Laigle et al. (2016) COSMOS2015 catalog. Specifically, this catalog includes the u* band from MegaCam on the Canada–Hawaii–France Telescope CHFT (Boulade et al. 2003), 6 broad optical bands (B, V , r, i⁺, z⁺⁺), 12 median, and 2 narrow bands from the Suprime-Cam on the Subaru telescope (Taniguchi et al. 2007, 2015), as well as the Y band from the Hyper Suprime-Cam/Subaru (Miyazaki et al. 2012), which is slightly bluer than the Y filter from VIRCAM. The COSMOS2015 also includes near-IR YJHK_s  band from VIRCAM on the VISTA telescope (Beard et al. 2004; McCrackene et al. 2012), and the H and K_s bands from WIRCam/CFHT (Puget et al. 2004; McCracken et al. 2010). Finally, the catalog includes the mid-IR Spitzer IRAC channels centered at 3.6 and 4.5 μm (ch1 and ch2, respectively; Fazio et al. 2004). For our analysis, we need a simulated mock of the HSC-Deep survey (Aihara et al. 2017) and its ancillary surveys, as described in more detail in Section 2.2. To provide a similar photometric coverage, we use the u, r, i, z, Y, J, H, K_s , ch1, and ch2 photometric data points from the Horizon-AGN COSMOS-like mock. However, the COSMOS-like mock is missing the g band. Therefore, we also use the public Euclid+LSST-like mock that includes the LSST g filter, which is very close to the HSC g filter. From these simulated photometric catalogs, we use only the photometry without any photometric errors applied, since we add photometric errors ourselves for the specific surveys we are considering (see Section 4.3). In cross-matching the two mocks, ≈100 objects are lost, leaving us with 789,292 galaxies. In the simulation, the u band exhibits a tail that extends to highly unrealistic values (corresponding to this band dropping shortward of the Lyman break). To mitigate this issue, we consider only galaxies with m_u < 28, leaving us with 691,486 objects. This cut limits the redshifts to z ∼ 3.4, as shown in Figure 1 (top), but has no other effects on the data set.

Zoom In Zoom Out Reset image size

2.1.1. Color Calibration

We expect the magnitude and color distributions of this lightcone to not be quite right, because of the simplicity of the dust attenuation treatment, as well as the fact that this idealized photometry does not take into account all the systematic effects that one would have from extracting photometry from an image (Laigle et al. 2019). We correct for this systematic offsets in the photometry following the method of Pearl et al. (2022). This method adjusts simulated magnitudes to match the correlation between M/L_ν,obs ratios and sSFR_UV,obs, ensuring consistency between simulations and observations.

We performed the color calibration using the UltraVISTA catalog (UVISTA; Muzzin et al. 2013), which has similar photometric coverage as well as all necessary derived quantities including redshift. To ensure maximum reliability for our color calibration, we selected galaxies from the UltraVISTA catalog based on the following criteria: 0.2 < z _peak < 4.0, 5 $\lt \mathrm{log}{M}_{* }\lt 12$ [M_⊙], star=0, contamination=0, nan _contam < 3, use=1, K_s < 23.5 AB. This results in a K_s -limited catalog with only ∼100 objects at z > 3.

We derived the SFR_UV for both the simulated data set and UVISTA following the Kennicutt & Evans (2012) conversion relation from L_NUV, which was available for both data sets. For the sake of the calibration, we excluded simulated galaxies with log(sSFR_UV) < −13 and log(sSFR) < −13, as such values are not typically observed.

We perform the color calibration in 70 redshift bins. The first 69 bins cover the range 0 < z ≤ 2.5, each with ≈8 K simulated galaxies, making k correction negligible within these narrow redshift bins. The remaining 132,281 simulated galaxies at z > 2.5 are all grouped in the final bin. The choice of this binning is motivated by the use of the K_s -selected UVISTA catalog for the calibration, indicating that there are few observed objects at z > 2.5. Consequently, dividing this redshift range into smaller bins would yield too few objects for a reliable fit in the log(M/L_ν )–log(sSFR_UV) plane. Figure 2 shows as an example the M/L_ν,u ratio for the u band versus the sSFR_UV for the redshift bin 0.8 < z < 0.9. We chose the u band as an example because the necessary calibration is by far the strongest in this band. As shown in Figure 2, prior to color calibration, the slopes of the fitted lines in the log(M_*/L_u )–log(sSFR_UV) plane differ between the Horizon-AGN mock and UVISTA data. To align the simulated fitted line with the observations, we perform linear fits to both distributions and calculate Δ_m (slope difference) and Δ_b (y-intercept difference). The new luminosity in each band and redshift bin is derived as

$$\begin{eqnarray}&&{L}_{\nu ,\mathrm{new}}=\displaystyle \frac{{L}_{\nu ,\mathrm{old}}}{{10}^{{{\rm{\Delta }}}_{b}}\times {\mathrm{sSFR}}_{\mathrm{UV}}^{{{\rm{\Delta }}}_{m}}},\end{eqnarray} \tag{ 1 }$$

where L_ν,old is the old luminosity in a specific band ν. We then derive new apparent magnitudes from

$$\begin{eqnarray}&&{m}_{\nu ,\mathrm{new}}={m}_{\nu ,\mathrm{old}}-2.5\times \mathrm{log}({L}_{\nu ,\mathrm{new}}/{L}_{\nu ,\mathrm{old}}).\end{eqnarray} \tag{ 2 }$$

Zoom In Zoom Out Reset image size

In this way, we obtain the color-calibrated version of the mock data set. For the rest of this paper, we refer to this mass-limited, color-calibrated mock as the mass-limited simulated data set, and it is what we use to train the SOM (see Section 3) as well as to generate the observations-like sample as described in the next subsection.

From the mass-limited simulated data set described above, we construct a simple K_s -limited sample that is meant to mimic the Hyper Supreme-Cam (HSC) Deep survey and ancillary data.

The HSC Deep survey is one of the three layers of the HSC survey (Aihara et al. 2017). The HSC is a 1.8 deg² field-of-view imaging camera (Miyazaki et al. 2017) mounted on the 8.2 m Subaru telescope on the summit of Maunakea in Hawaii, operated by the National Astronomical Observatory of Japan (NAOJ). The product of its wide field of view and the large collecting area of the telescope makes HSC second only to the anticipated LSST Camera (Aihara et al. 2017). The Deep layer, which is our focus, covers 27 deg² spread across four well-studied extragalactic fields (XMM-LSS, COSMOS, ELAIS_N1, and DEEP2-3). The HSC photometry includes five broad bands (g, r, i, z, and Y) and three narrowband filters (NB387, NB816, and NB921).

The HSC Deep joint data set augments the above through the addition of several complementary surveys. The u-band coverage is provided by the CLAUDS survey (Sawicki et al. 2019), which uses the Megaprime camera on the CFHT telescope. The near-IR J, H, and K_s bands are supplied from a variety of surveys, including UKIDSS/DXS and UDS (Lawrence et al. 2007), VIDEO (Jarvis et al. 2012), UltraVISTA (Muzzin et al. 2013), and DUNES ⁷ (in preparation). Spitzer IRAC imaging in channels 1 and 2 (3.6 and 4.5 μm) is supplied from the DeepDrill survey (Lacy et al. 2020), the SHIRAZ survey (Annunziatella et al. 2023), as well as a number of surveys within COSMOS (e.g., Sanders et al. 2007).

At the time of writing this paper, the precise properties of the multi-band photometric catalog for this HSC-Deep joint data set had not been determined. To test the degree to which our HSC-like mock catalog depends on the precise catalog construction method, we explored two extremes in such a catalog construction. The first extreme, which results in the smallest number of galaxies included, required a 5σ detection in each band from the u band through to IRAC ch2. Our second extreme simulates a scenario where one uses a reference image for source detection and performs forced photometry on all other images, thus producing a photometry value regardless of whether the brightness of any given source exceed 5σ in these other bands or not (e.g., McCrackene et al. 2012). This method would result in simpler selection function (a single magnitude cut) and a larger number of galaxies than the first method. For this method, we chose to make a simple K_s cut of K_s < 23.5 without any additional cuts. We find that 75% of the K_s -limited sample are also within the first HSC-like mock where we required a 5σ cut in all bands. Thus, for the bulk of the galaxies, our results are independent of the specifics of the catalog construction. While we ran our analysis on both HSC-like mocks, we only show in this paper the results of the runs with the K_s -limited mock, because of its more extended redshift and mass coverage. There was no significant difference in our findings on parameter estimation accuracy or effects of missing data between the two HSC-like mocks. Our adopted HSC-like K_s -limited data set contains 232,118 galaxies, and its stellar mass versus redshift distribution is shown in the bottom panel of Figure 1.

The bottom panel of Figure 1 shows that the K_s -magnitude cut introduces an incompleteness in stellar mass as a function of redshift, with the stellar mass limit increasing with increasing redshift. As a result, low-mass galaxies are excluded from the sample at z ≳ 1, in particular those with log(M_*/M_⊙) < 10 are entirely absent by z ∼ 3. We caution that the simulated lightcone we started with represents 1 deg², whereas the HSC-Deep survey covers 27 deg². Rarer populations such as log(M_*/M_⊙)>11 galaxies at z > 3 can be present in the real data even though they are largely missing here.

A Self-Organizing Map (SOM; Kohonen 1982) is an unsupervised artificial neural network that performs a dimensional reduction of a multidimensional parameter space to a lower-dimensional space, while preserving the topology of the data. It thus consists of a lower-dimensional grid where each pixel is characterized by a weight vector that represents the mapping from the higher-dimensional to the lower-dimensional space. The training involves determining these "weights" so that objects that are similar in the high-dimensional parameter space are grouped together, maintaining the intrinsic structure of input data. For this reason, the SOM is a powerful visualization tool, especially for a final 2D grid, as considered in this work. We use the Python library SomPY ⁸ (Moosavi et al. 2014) for constructing and training our SOM. Below, we describe in more detail how the algorithm works.

The first step is selecting a representative training data, which is then normalized by SomPY to the unit variance with a mean of zero. The weight vector of each SOM pixel, or neural network, in SomPY is initialized using principal component analysis (PCA; Chatfield & Collins 1980) to bring the initial weights closer to the input data (SomPY also allows for random initialization, but it extends the training duration). At this point, the Euclidean distance between each neurons and input data is calculated. The neuron with the smallest distance to an input data becomes the Best Match Unit (BMU). The latter is thus dragged closer to the data point, and neurons around the BMU, specifically those within a "neighborhood radius," are also dragged in the same direction. This process is repeated in an iterative way until all neurons are as close as possible to the input data. This process is unsupervised, as it does not require any a priori labeling of input data, unlike, for example, traditional neural networks where the weights are optimized to match some particular outputs vector. After this optimization, similar objects in the high-dimensional parameter space will be grouped together in the final lower-dimensional (usually 2D) space, maintaining the intrinsic topology of input data.

In this work, the input parameter space has 10 dimensions representing observed-frame colors for galaxies. These colors, based on the expected photometry available for the HSC-Deep joint data set (see Section 2.2) are: u − g, g − Y, g − J, r − i, r − z, r − Y, i − z, J − H, K − ch2, and ch1 − ch2. We select these colors for their sensitivity to galaxy parameters such as redshift, normalized M_* and SFR, and sSFR. To do so, we computed the distance correlation for all possible two-color combinations (55 in total) between the 11 available bands. This analysis identified the 10 colors with the strongest nonlinear relationships with our target parameters. While we initially explored using adjacent colors for SOM training, our final selection minimized parameter spread within SOM pixels (Section 3.4), leading to more accurate parameter estimations, in particular for M_* and SFR.

For our training set, we selected a random subset from the mass-limited data set. The training set consists of 228,524 galaxies or about 1/3 of the total. The performance of SOM depends on the user set map's size and geometry, which need to be verified by the user depending on their specific case (see Davidzon et al. 2019). The choices of size and geometry of the map should represent a compromise between good sampling of the data and a high resolution (less quantization error). We achieve this, following the Davidzon et al. (2019) procedure, by training different SOMs, starting with a square 20 × 20 map and gradually increasing the size by adding 10 cells in both dimensions. We kept track of sampling versus error, finding that an 80 × 80 SOM represents a good compromise, similar to the findings of Davidzon et al. (2019). We also tested the effect of different aspect ratios, but found that a square map gave the best performance. Finally, in our final 80 × 80 SOM, we verified that the 10-color distributions in the input data set and across SOM pixels are consistent. This ensures that an SOM trained via this method is representative of the input data. ⁹

Figure 3 shows the distribution of each of the 10 normalized colors across the map. The distribution is typically smooth except for a sharp discontinuity, in the middle of the SOM, that is more evident for some of the colors such as u − g or ch1 − ch2. This is because, in this case, the algorithm organized the input data in a ring-like structure, which causes the colors to change in a counterclockwise direction. The point where the loop closes marks the discontinuity. However, this discontinuity is not random, but rather it has a physical meaning, denoting the transition between galaxies with low and high redshift and other stellar population properties. For instance, for ch1 − ch2 we expect redder colors at z > 1.7 when both IRAC channels sample blueward of the 1.6 μm bump. This means that lower values of ch1 − ch2 on the map would correspond to z < 1.7 and vice versa. The u − g colors also switch from low to high across the same discontinuity, likely due to the effects of redshift and SFR. Lower redshifts correspond to lower-mass, passive galaxies, while higher redshifts correspond to redder colors, indicating more massive, star-forming galaxies.

Zoom In Zoom Out Reset image size

At z > 1.7, the u − g colors redden further as a result of the Lyman break, which enters the u filter at z > 2.2. In general, both redshift and galaxy stellar population properties affect these distributions. Figure 3 clearly demonstrates that galaxies with similar 10-color sets are grouped together in the trained SOM.

Considering our ultimate goal is to replicate the methodology presented in this work with observed data, we conduct an additional sanity check, shown in the Appendix. This check involves verifying the similarity of the component maps of two different SOMs, one trained with an observed data set (UVISTA) and the other with simulated data where we match as closely as possible the selection function of the observed data set. We find encouragingly consistent component maps.

The color distributions across the SOM (Figure 3) already allowed us to identify trends with redshift and stellar population parameters. To quantify these trends and use SOM for parameter estimation, we need to label each SOM pixel with the galaxy properties based on its corresponding 10-dimensional color combination. There are multiple ways to approach this task. Masters et al. (2015) use the subset of their galaxies with spectroscopic redshifts to label each pixel by the expected redshift. Hemmati et al. (2019) use simple empirical template libraries covering a range of parameters such as age of the stellar population and dustiness. In our case, following the approach of Davidzon et al. (2019), we label SOM pixels with the median of each parameter. The latter are the Horizon-AGN parameters from the training set galaxies falling within a given pixel. This works for redshifts and sSFRs, which only depend on the observed-frame colors. Stellar mass M_* and SFR, however, require information about the overall amplitude of a galaxy's spectrum, which is not modeled in our color-trained SOM. To account for this amplitude dependence, we normalize M_* and SFR values from our simulated data set to a reference magnitude J = 23 mag:

$$\begin{eqnarray}&&{nX}=X\times {10}^{-0.4\times (23.0-J)},\end{eqnarray} \tag{ 3 }$$

where X is the galaxy parameter that needs to be normalized and J is the J magnitude for each particular galaxy.

The left panels of Figure 4 show how redshift (z), normalized stellar mass (n M_*), normalized star formation rate (nSFR), and sSFR vary across the SOM. Consistent with our initial visual inspection, the considered parameters increase in a counterclockwise trend. Specifically, more distant galaxies at z > 2.5 are also more massive and active. Note that the SOM's morphology is not universal; it varies depending on the specific colors and data set used for training. However, its ability to reveal trends remains consistent. For instance, comparing these plots with Figure 3, it is evident that redder galaxies are the farther, more massive ones.

Zoom In Zoom Out Reset image size

The right panels of Figure 4 show the distribution width δ for each parameter. The width is quantified by the difference between the 84th and 16th percentiles among galaxies within each pixel. The δ/2 represents the systematic uncertainties σ_sys of the derived redshifts and stellar population pararameters using SOMs. Large parameter spread within a pixel can be due to degeneracies, although the use of 10 different colors tends to mitigate this to a large degree. We also notice that large spreads in our SOM tend to concentrate the edges. This could be due to (1) boundary effects, as Davidzon et al. (2019) explains, arising from galaxies with extreme colors being placed on the borders of the map by the algorithm; or (2) poorly defined areas, where there are not enough objects to properly train a cell.

One of the main strengths of SOMs is its ability to provide interpretable visualization of high-dimensional data. In particular, once an SOM is trained using the mass-limited data set, it is possible to project on the map different samples derived by applying different cuts. Figure 5 shows the projection of the K_s -limited sample, or HSC-Deep joint-like sample, reaffirming the conclusion of Figure 1. Indeed, a lack of objects at z ≳ 2.5 with log(M_*/M_⊙)<10 is present on the SOM. This is due to the stellar mass incompleteness introduced by the K_s -magnitude cut. Additionally, we notice that the "missing" galaxies in the HSC-Deep joint-like sample are predominantly higher sSFR galaxies with log(sSFR)≳−9.2 yr⁻¹. In observed terms, these are galaxies with blue optical (e.g., g − Y and r − i), but red near/mid-IR (e.g., J − H and ch1 − ch2) colors. This projection of the K_s -limited sample highlights the effectiveness of SOMs in visualizing the impact of both magnitude and color cuts.

Zoom In Zoom Out Reset image size

We test the use of SOMs in deriving redshift and stellar population parameters by projecting the K_s -limited data set onto the SOM. For each galaxy, we assign the median of z, log(n M_*), log(nSFR), and log(sSFR) of the SOM pixel in which it resides. To denormalize M_* and SFR, we use Equation (3) and the observed J-band magnitude for each galaxy. Note that, while we employ the K_s -limited sample for predictions, the SOM was trained and its pixels labeled with a different data set, specifically the training data set, which was randomly drawn from the mass-limited data set. Although the distributions of these two data sets, in terms of z, M_*, and SFR may differ, we expect the relation between the mapping color space and physical parameter space to be consistent.

Figure 6 shows the comparison between the true (in values from the HorizonAGN simulation) and the SOM-derived values for z, M_*, SFR, and sSFR. For each parameter, we also show four statistics: normalized median absolute deviation (σ_NMAD), outlier fraction η, root mean squared error (RMSE), and bias. The outlier fraction η is the percentage of objects that satisfy $| {z}_{\mathrm{SOM}}-{z}_{\mathrm{sim}}| /(1+{z}_{\mathrm{sim}})\gt 0.15$ in the case of z_phot, and $| \mathrm{log}{X}_{\mathrm{SOM}}-\mathrm{log}{X}_{\mathrm{sim}}| \gt 2\sigma$ in the case of the other properties X. Here, σ is the standard deviation of the difference log(predicted) – log(true).

Zoom In Zoom Out Reset image size

We find that redshift outperforms other statistics, with σ_NMAD and RMSE of 0.01 and 0.03, respectively. Stellar mass estimation is also accurate, with σ_NMAD and RMSE of 0.06 and 0.07, respectively. As is typical of all parameter estimation methods, the estimation of SFR and sSFR has higher uncertainties. The spread between the true and SOM-derived parameters in Figure 6 is due to the intrinsic spread δ in these parameters within an SOM pixel as shown in Figure 4. This spread does not include the effect of photometric errors. Therefore, the performance statistics presented here represents the ideal scenario.

To consider photometric uncertainties, we adopted relative errors per magnitude derived from the HSC DR2 preliminary u2K catalog ¹⁰ for the u through K_s bands (Aihara et al. 2019). We focused on the XMM-LSS field, since its IRAC ch1 and ch2 fluxes and uncertainties were available from Krefting et al. (2020) and Nyland et al. (2023). We computed the relative photometric errors σ_ν in each band by fitting a line between log(σ_ν /f_ν ) and magnitude in that band. This approach accounts for larger relative errors for the more abundant faint objects and vice versa.

For each galaxy in our simulated HSC-like data set, we drew a new random flux value from a Gaussian distribution centered at the true flux and a width given by the photometric error. We projected these new perturbed data onto the SOM trained with noiseless colors. Afterward, we repeated the procedure to estimate parameters as described in Section 4.2. Note that, to denormalized M_* and SFR, we used the J-band magnitudes with added noise.

Figure 7 shows the comparison of true and SOM-derived parameters, accounting for photometric uncertainties. As expected, parameter uncertainties increase with photometric errors, but the estimations remain reasonably accurate. For instance, the new redshift σ_NMAD is 0.026 with only a 1.2% outlier fraction. The true versus SOM-derived z_phot distribution in Figure 7, however, shows horizontal stripes where the z_phot derivation is less reliable, with increased numbers of outliers, at specific redshift ranges. The most evident stripe is at z_SOM ∼ 1–2. We investigated the primary cause of these effects by examining the position on SOM of these outliers. The latter are all located near the area of the abrupt change in K − ch2 and ch1 − ch2 (Figure 3), as a result of the move across the 1.6 μm bump. The larger photometric errors in the IRAC channels can cause galaxies to fall into neighboring SOM pixels associated with different redshift values. We validated this hypothesis by repeating the same analysis of this section but removing the IRAC photometric errors. As a result, the horizontal stripes visible in the upper right panel in Figure 7 disappeared, confirming that the colors associated with the IRAC-channels were indeed the primary contributors. This illustrates how photometric noise in specific bands impacts derived parameter uncertainties.

Zoom In Zoom Out Reset image size

Due to the degeneracies between photometric redshift and shape of the SED, the outliers of the horizontal stripes are also responsible for the outliers in mass and SFR derivation visible in Figure 7. For example, if objects at somewhat higher z are incorrectly placed at lower z (right side of the stripe), their redder colors are associated with more quiescent galaxies, leading to an underestimation of SFR, and vice versa. Similarly, if an object is placed at a redshift lower than its true value, the observed J band suggests a lower stellar mass than in reality. This results in stripes below the 1:1 relation in M_* (see Figure 7). Since sSFR is the ratio of SFR and M_*, these effects are canceled out. This explains why our estimates for sSFR changed the least between the no-error and error cases (Figures 6 versus 7). We tested this interpretation by running the analysis without any IRAC channel errors, and as a result, the horizontal stripes in z_phot, in M_* ($\mathrm{log}{M}_{* ,\mathrm{SOM}}/{M}_{\odot }\sim 9.5\mbox{--}10$), and in SFR (log(SFR)_SOM ∼ 0) disappeared.

Figure 7 gives an overall sense of parameter uncertainties when photometric errors are considered. These parameter uncertainties can be statistical, resulting from photometric errors, and systematic, resulting from parameter value distribution within a single SOM pixel. Below, we explain how we compute both statistical and systematic parameter uncertainties for individual galaxies.

For any galaxy, the likelihood ${ \mathcal L }$ of it occupying an SOM pixel with coordinates [x, y] is computed as

$$\begin{eqnarray}&&{{ \mathcal L }}_{[x,y]}={e}^{-\tfrac{{\chi }_{[x,y]}^{2}}{2}},\end{eqnarray} \tag{ 4 }$$

with

$$\begin{eqnarray}&&{\chi }_{[x,y]}^{2}=\displaystyle \sum _{i=1}^{10}{\left(\displaystyle \frac{{C}_{i}-\overline{{C}_{i}}}{{\sigma }_{{C}_{i}}}\right)}_{[x,y]}^{2},\end{eqnarray} \tag{ 5 }$$

where, in the pixel [x, y], C_i are the 10 colors of the test galaxy, $\overline{{C}_{i}}$ are the 10 median colors in that pixel, and ${\sigma }_{{C}_{i}}$ represent the uncertainties of C_i . The likelihood surface is normalized by dividing ${{ \mathcal L }}_{[x,y]}$ by its sum. Figure 8 (left panels) shows the likelihood surfaces for two test galaxies. These test galaxies are of similar redshift (z ∼ 1.2) and were chosen to represent a high-S/N (signal-to-noise ratio) galaxy (pixel = [61, 21], S/N = 82.3 in g band, m_g = 23.4) and a low-S/N galaxy (pixel = [71, 29], S/N = 4.0 in g band, m_g = 27.7). The likelihood surface of the high-S/N test galaxy is very concentrated around its true position on the SOM, while for the low-S/N test galaxy it is more spread out.

Zoom In Zoom Out Reset image size

From these likelihood surfaces, combined with the pixel labels for a given parameter, we can derive probability distribution functions for any of our parameters. For example, the right-hand panels of Figure 8 show the redshift probability distribution functions P(z) derived from these likelihood surfaces. The P(z) of the low-S/N test galaxy exhibits a larger spread compared to the high-S/N test galaxy. This confirms the idea that the accuracy of parameter estimation with SOMs depends on a galaxy's S/Ns, as indeed is the case for all methods used in deriving galaxy parameters.

We define the statistical uncertainty for the given parameter, σ_stat,param, for any individual galaxy as the standard deviation of the probability distribution function, such as the P(z) functions discussed above. To estimate the parameter systematic uncertainty σ_sys,param, we take half of the 〈δ〉 = 84%–16%, the dispersion of galaxy parameter within the best-fit pixel for that galaxy (see right panels of Figure 4). The total uncertainty is

$$\begin{eqnarray}&&{\sigma }_{\mathrm{tot},\mathrm{param}}=\sqrt{{\sigma }_{\mathrm{stat},\mathrm{param}}^{2}+{\sigma }_{\mathrm{sys},\mathrm{param}}^{2}}.\end{eqnarray} \tag{ 6 }$$

For example, to derive the uncertainties in the redshift estimation for the high-S/N test galaxy, we look at its pixel [61, 21], which has a redshift spread 〈δ〉 = 0.057, hence σ_sys = 0.029, leading to a total redshift uncertainty of ${\sigma }_{\mathrm{tot},z}=\sqrt{{0.033}^{2}+{0.029}^{2}}=0.044$. To compare this with the nominal σ_NMAD in Figure 7, we need σ_tot,z /(1 + 1.2) = 0.02, which is consistent with the σ_NMAD of 0.026.

The previous section's analysis assumes flux measurements are available across all bands. In reality, the surveys making up the HSC-Deep joint data set do not perfectly overlap. Several square degrees of the HSC-Deep joint survey in fact will have no data outside the core bands of grizY. Additionally, even in overlapping areas, there can be instances of missing data in individual bands, due to issues such as bad data arising from image artifacts. The simplest approach would be to train an SOM excluding any missing bands. We tested this scenario for multiple individual missing bands or sets of missing bands. Figure 9 shows the z_phot σ_NMAD in a case where the u band, J band, H band, K_s  band, JHK_s , K_s +IRAC, and uJHK+IRAC are missing. The σ_NMAD gets worse compared to the baseline of complete data (dashed line) with σ_NMAD = 0.026. Note that this analysis takes into account photometric uncertainties in the same manner as done in Section 4.3, where each flux is perturbed by Gaussian noise based on the given flux uncertainty.

Zoom In Zoom Out Reset image size

The degradation in parameter estimation seen in Figure 9 is due to the loss of information resulting from the missing bands. To address this, we assume that galaxies with missing data fall within the same color distribution as the training sample. This is justifiable when the missing data is simply due to the inhomogeneous coverage of the various multiband surveys. To recover the colors involved with missing bands, we perform 500 random draws from the color distribution of these same colors, as defined in the training set. Our color distribution again includes photometric uncertainties as discussed above. The random draws involve either 1D, 2D, or 3D color distributions, depending on whether the missing band is on the edge (e.g., u band) or in the middle (e.g., H band), or if a set of adjacent bands is missing (e.g., JHK). With our 500 random draws for the "missing band" colors combined with the existing colors for each galaxy, we are able to essentially build an SOM likelihood surface for that galaxy in a manner similar to that discussed in Section 4.4. The adopted parameter value is then simply based on the maximum-likelihood SOM pixel. The resulting redshift σ_NMAD values for each of our missing data scenarios are shown in Figure 9. We note that, in most cases, this band recovery procedure allows us to achieve significantly better σ_NMAD values than by simply training the SOM without the missing data. We emphasize that this improvement is based on the assumption that the color distributions determined from the data with complete coverage are representative of the data with incomplete coverage as well.

The only cases where our recovered data perform worse than the SOM trained by omitting the missing data are those involving IRAC channels. We tested that this worsening does not depend on the adopted photometric error. This effect is likely due to a more uncertain extrapolation when more than one color needs to be recovered right on the edge of the wavelength coverage, as is the case when the IRAC channels are missing.

The missing-data analysis also provides us with information on the relative importance of having particular bands present for the recovery of our parameters. Table 1 goes beyond z_phot, shown in Figure 9, to explore how the σ_NMAD for z_phot, M_*, and SFR all vary depending on which bands are missing and filled in with our recovery method (Section 4.5). We reach reasonably good accuracy for when only a single band is missing. As more bands are missing and filled in, the values of σ_NMAD increase. For example, the worst parameter estimation occurs when all ancillary bands outside the HSC-Deep core bands (grizY), namely uJHK+IRAC, are missing and filled in with our method. If at least the u band, among the ancillary bands, is not missing (i.e., in the case of recovered missing JHK+IRAC) the parameter σ_NMAD improves, especially for z_phot estimation, with a decrease in σ_NMAD of ∼32%. These findings show the critical importance of the ancillary surveys providing the u-band, the near-IR, and the Spitzer IRAC ch1 and ch2 coverage.

The variation of parameter σ_NMAD with specific bands depends on their different sensitivities to galaxy parameters at different ranges. Figure 10 illustrates this for the cases when u-band (top panel), H-band (middle panel), and all ancillary to HSC-Deep data are missing (bottom panel). We chose these specific cases because the u band is crucial for accurate z_phot estimation at z ≳ 2.5 (as this filter enters the Lyman break) and traces the rest-frame ultraviolet, hence it is relevant for the SFR estimation. The top panel in Figure 10 shows that indeed, when the u band is missing, our redshift recovery gets worse beyond z ∼ 2.5, and the SFR recovery is poor. As seen in Table 1, the SFR recovery here is the worst compared to any other single (or couple of) missing-band cases.

Zoom In Zoom Out Reset image size

The H band, in our redshift range, traces the rest-frame optical to near-IR, therefore it is sensitive to the stellar mass. It also affects the z_phot determination, since that is sensitive to the 1.6 μm bump at lower z, as well as to the 4000 Å break at z > 3. The middle panel in Figure 10 shows that indeed, when the H band is missing, our redshift recovery is much worse at z > 2. The underestimated redshifts in this regime also lead to a pronounced bulge of underestimated stellar masses—and corresponding underestimated SFRs.

The last case we explored in more detail was the worst-case scenario (see Table 1), which was the missing uJHK+IRAC scenario, or when we have no data outside the HSC-Deep core bands (grizY). In this case, we use the scenario where the missing data are simply omitted rather than recovered, since the latter shows worse parameter recovery (see Figure 9). We see, in the bottom panel of Figure 10, that the redshift estimation is still reasonable at z < 1–1.5, where the optical bands trace the 4000 Å break, but there is essentially no handle on the redshift beyond that (up to the z ∼ 3 regime explored here). The longer-wavelength (J+) bands are critical for the z ∼ 1.5–3 estimation, due to their sampling both the Balmer break and the 1.6 μm stellar bump.

Besides cases of missing data, we can also have upper limits if the source is too faint relative to the noise in a given image. The naive first approach to handle such upper limits is to draw random fluxes from uniform distributions representing the upper limit. However, this method can lead to unrealistic colors that place the galaxy outside the SOM surface, which only represents the 10D color space present in the training set. Dealing with upper limits is, in general, a significant challenge for machine-learning methods, as it requires extrapolating outside the training set.

While a more comprehensive treatment of upper limits in SOMs is beyond the scope of the present work, we can adapt our missing data approach for one particular scenario involving upper limits. This is the case where upper limits are a result of some part of the survey having significantly shallower coverage in some band relative to another deeper portion of the survey. In this case, we can train the model using the full color set from the deep portion of the survey. The assumption is that the deeper-data training set contains realistic colors for sources with upper limits in the shallower portion. In this case, we can implement a modification of the missing data approach described in the previous section. Instead of drawing from the full color distribution, we draw from a truncated color distribution that accounts for the detection limit in the shallower data. Figure 11 illustrates this approach with an example of a test galaxy having only an upper limit in the H band. The test galaxy is drawn from a fiducial data set with H-band coverage two orders of magnitude shallower coverage than in the training set. In this case, the recovered colors are realistic and fall within the panchromatic space used to train our SOM, as shown by the gray region in Figure 11. Note that the tightness of the SOM likelihood survey of this upper limit galaxy is due to the fact that each SOM pixel involves 10 colors, which helps compensate to a large degree for any missing/upper limit data.

Zoom In Zoom Out Reset image size

This paper was originally inspired by Davidzon et al. (2019), who also explored the derivation of z_phot and SFR using SomPY on the Horizon-AGN simulated lightcone. They added photometric errors corresponding to the COSMOS photometry. Their training color set excluded g-band colors and colors from nonadjacent bands (e.g., g − Y, r − z). In this paper, we expand on this work, first by applying it to a different photometric set (the HSC-Deep ones). Next, we explore the effect of different color choices for training the SOM, and consider the effects of missing data and how to handle them. Fundamentally, we achieve a photo-z σ_NMAD of 0.026 and an outlier fraction of 1.2%, while Davidzon et al. (2019) reported a σ_NMAD of 0.044 and an outlier fraction of 6.1%. While the approach is very similar, this is a significant improvement in performance. Our tests suggests that the key aspect behind this improvement is the specific choice of colors. We use more widely separated colors, such as r − z, which is sensitive to the redshift, and we also add the g band (which is omitted in Davidzon et al. (2019)), which also improves performance, especially in redshift regimes where it helps pinpoint the Balmer break. In addition, our improved performance is aided by our uncertainties being ≈0.2 mag lower in the u band, comparable in the visible bands, and ≈0.1 mag lower in the NIR, although our IRAC channels errors are worse by ≈0.2 mag. Davidzon et al. (2019) do not quote detailed statistics for SFR, but they quote σ_NMAD < 0.2, which is consistent with our results.

In a more recent paper, Davidzon et al. (2022) used their SOM method to derive M_* from the observed COSMOS data and compare their SOM-derived values with traditional SED-fitting. Their results indicate that this method works equally well as traditional SED fitting, but it is computationally much faster, making it suitable for large data sets. Their distribution of SOM-derived versus SED-fitting-derived stellar masses appears visually wider than our SOM-derived versus true stellar mass. However, this is not a fair comparison, since there are uncertainties associated with the SED-fitting-derived values as well.

In this paper, we explore for the first time in the astronomical context the handling of missing data in SOMs. This issue has been explored in SED template fitting codes, e.g., by modifying the χ² method, usually at the expense of increased computational cost (Sawicki 2012; Boquien et al. 2019). A recent paper by Rejeb et al. (2023) does discuss how to train SOMs with incomplete data; however, this paper did not appear until the final stages of our current work, therefore we have not had a chance to implement it. In future works, we plan to explore the possible use of their algorithm through their publicly available code. ¹¹ One needs to be cautious, however, about including incomplete data in the training, especially in cases of large variance, as discussed in Cottrell & Letrémy (2007). They argue that, when a sufficient amount of data are available to build a robust model, it is preferable to train on the complete data and only afterward consider the effects of missing data. Indeed, this is the approach we adopted in our study. In our paper, we outline a method to train the SOM on the complete data set and "recover" the missing data where the data are incomplete by assuming the SOM training set is otherwise representative of the galaxies with missing bands. In line with our finding that it is better to fill in the missing data rather than simply omit it, Chartab et al. (2023) achieved better results in SED fitting after first using a random forest model to fill in any missing data.