# Hierarchical Bayesian inference of photometric redshifts with stellar
population synthesis models

Boris Leistedt Department of Physics, Imperial College London, Blackett
Laboratory, Prince Consort Road, London SW7 2AZ, UK  Justin Alsing Oskar Klein
Centre for Cosmoparticle Physics, Department of Physics, Stockholm University,
Stockholm SE-106 91, Sweden  Hiranya Peiris Department of Physics and
Astronomy, University College London, Gower Street, London, WC1E 6BT, UK
Oskar Klein Centre for Cosmoparticle Physics, Department of Physics, Stockholm
University, Stockholm SE-106 91, Sweden  Daniel Mortlock Department of
Physics, Imperial College London, Blackett Laboratory, Prince Consort Road,
London SW7 2AZ, UK  Department of Mathematics, Imperial College London, London
SW7 2AZ, UK  Oskar Klein Centre for Cosmoparticle Physics, Department of
Physics, Stockholm University, Stockholm SE-106 91, Sweden  Joel Leja
Department of Astronomy & Astrophysics, The Pennsylvania State University,
University Park, PA 16802, USA  Institute for Computational & Data Sciences,
The Pennsylvania State University, University Park, PA, USA  Institute for
Gravitation and the Cosmos, The Pennsylvania State University, University
Park, PA 16802, USA

###### Abstract

We present a Bayesian hierarchical framework to analyze photometric galaxy
survey data with stellar population synthesis (SPS) models. Our method couples
robust modeling of spectral energy distributions with a population model and a
noise model to characterize the statistical properties of the galaxy
populations and real observations, respectively. By self-consistently
inferring all model parameters, from high-level hyper-parameters to SPS
parameters of individual galaxies, one can separate sources of bias and
uncertainty in the data. We demonstrate the strengths and flexibility of this
approach by deriving accurate photometric redshifts for a sample of
spectroscopically-confirmed galaxies in the COSMOS field, all with 26-band
photometry and spectroscopic redshifts. We achieve a performance competitive
with publicly-released photometric redshift catalogs based on the same data.
Prior to this work, this approach was computationally intractable in practice
due to the heavy computational load of SPS model calls; we overcome this
challenge using with neural emulators. We find that the largest photometric
residuals are associated with poor calibration for emission line luminosities
and thus build a framework to mitigate these effects. This combination of
physics-based modeling accelerated with machine learning paves the path
towards meeting the stringent requirements on the accuracy of photometric
redshift estimation imposed by upcoming cosmological surveys. The approach
also has the potential to create new links between cosmology and galaxy
evolution through the analysis of photometric datasets.

photometric redshifts - galaxy surveys - cosmological parameters

††journal: APJS

##  1 Introduction

Deriving accurate redshifts and redshift distributions from photometry alone
is of central importance for the scientific exploitation of modern galaxy
surveys. In particular, cosmological analyses involving galaxy clustering and
weak gravitational lensing require exquisitely accurate estimates of the
redshift distributions of the selected galaxy populations (see e.g.,
Hildebrandt et al. 2020; Abbott et al. 2022). Photometric redshifts
(photo-$z$’s) of individual galaxies are also useful for selecting targets of
interest (e.g., at high redshift).

Photo-$z$’s are typically obtained via two types of methods: template fitting,
and machine learning111A third technique, clustering redshifts, can constrain
redshift distributions, but cannot deliver individual photo-$z$’s, see e.g.,
Schneider et al. (2006); Newman (2008); Ménard et al. (2013); McQuinn & White
(2013); Schmidt et al. (2013); Morrison et al. (2017); Gatti et al. (2021)..
Template fitting (e.g., , Benitez 2000; Ilbert et al. 2006; Brammer et al.
2008; Ilbert et al. 2006; Tanaka 2015; Acquaviva et al. 2015; Carnall et al.
2018; Battisti et al. 2019) relies on lists of grids or templates for the
spectral energy distributions (SEDs) of the galaxies considered, combined with
a prior on the parameters involved (galaxy type, magnitude, etc). Machine
learning techniques (e.g., Collister & Lahav 2004; Kind & Brunner 2013)
involve flexible models trained on spectroscopic or synthetic data. Despite
their past successes, both these approaches are unable to satisfy the accuracy
requirements imposed by upcoming Stage IV surveys, as they are not robust to
errors in the data (calibration offsets, underestimated uncertainties), in the
SED model, or to the non-representativity or sample variance of the
spectroscopic data used for training and validation (e.g., Hartley et al.
2020; Gatti et al. 2021; Myles et al. 2021; Newman & Gruen 2022). For example,
the Vera C. Rubin Observatory’s Legacy Survey of Space and Time (LSST) will
require the uncertainty on the mean of redshift distributions to be smaller
than $0.001(1+z)$ by its Year 10 data release (The LSST Dark Energy Science
Collaboration et al., 2018), which in turn imposes very stringent requirements
on the accuracy of individual photo-$z$’s.

We introduce a framework capable of addressing these needs: a hierarchical
model of galaxy photometry, consisting of three components describing a
generating process of survey data:

  * •

galaxy SEDs generated with stellar population synthesis;

  * •

a population model describing the distributions of the intrinsic parameters of
galaxies in the sample analyzed;

  * •

a noise model for generating photometric observations (i.e., a likelihood
function).

This provides a powerful combination of empirical modeling with physics-driven
components, incorporating known (astro)physics while leaving sufficient room
for data-driven corrections, addressing the aforementioned issues of
flexibility, interpretability, mis-calibration, and extrapolation outside of
training data. Thus, once applied and calibrated to large survey data-sets, it
has the potential to yield both precise and accurate inferences on intrinsic
galaxy parameters, including redshift. Building on previous work on
hierarchical modeling in photometric surveys (Leistedt et al., 2016; Jones &
Heavens, 2018; Sánchez & Bernstein, 2018; Rau et al., 2019; Alarcon et al.,
2020; Sánchez et al., 2020), our framework takes advantage of the increasing
accuracy and flexibility of SPS, and of methodological innovations in Bayesian
inference for large hierarchical models, in particular involving machine
learning-facilitated acceleration of computationally-intensive analysis steps.

In this paper, we demonstrate the efficacy of this framework by applying it to
the COSMOS2020 data set (Weaver et al., 2022). For the choice of SED model, we
consider Prospector-$\alpha$ (Leja et al., 2017, 2019; Johnson et al., 2021).
We address the high computational needs of its model calls with neural
emulators (Alsing et al., 2020) applied to both SPS total fluxes and emission
line contributions. We find that the emission line contributions are
responsible for the largest photometric residuals in our analysis and
therefore require calibration in order to derive accurate photo-$z$’s.
Together, the SED, data, and population models form a complex hierarchical
model which would normally be computationally intractable, a challenge we
overcome with neural emulator-based acceleration. We obtain redshift estimates
competitive with those publicly released with the COSMOS2020 data using this
baseline model, providing a benchmark validation test demonstrating the power
of the new methodology. The extensive 26-band wavelength coverage of
COSMOS2020 simplifies parameter inference and allows us to demonstrate the
conceptual advantages of the new methodology without the significant overhead
of carrying out a Markov chain Monte Carlo (MCMC) analysis. Application to
broadband-only photometric surveys, such as the LSST, will require a more
involved treatment of the uncertainties and selection effects, which we defer
to future work. In a companion paper (Alsing et al., 2022), we show how to
correctly include selection effects in the forward modeling of such surveys,
and in particular for the inference of galaxy redshift distributions.

The outline of this paper is as follows. In Sec. 2 we describe our framework
in the form of a generative model of photometric data. An inference formalism
is presented in Sec. 3. We describe the data in Sec. 4, and the results in
Sec. 5. We discuss the results in further detail in Sec. 6, and conclude in
Sec. 7. In what follows, a Planck 2015 (Ade et al., 2016) and a Chabrier
(2003) initial mass function are adopted for all relevant calculations.
$p(\cdot|\cdot)$ refers to the probability density of the quantity before the
“$|$”, conditioned on quantities behind it. We work with fluxes in AB units.

##  2 Generative model for galaxy photometry

Hierarchical models tackling photometric redshifts and redshift distributions
have been developed (Leistedt et al., 2016; Rau et al., 2019; Jones & Heavens,
2018; Sánchez & Bernstein, 2018; Alarcon et al., 2020; Sánchez et al., 2020),
but typically assume a fixed SED model. Leistedt et al. (2019) pioneered the
approach of inferring redshifts jointly with calibration of the hyper-
parameters of a SED model. However, like most photo-$z$ methodologies (Brammer
et al., 2008; Ilbert et al., 2006; Tanaka, 2015; Acquaviva et al., 2015;
Carnall et al., 2018; Battisti et al., 2019), the model still relied on a
curated set of SED templates, which made it difficult to construct more
structured, interpretable corrections to the SEDs or to the population model.
Our work introduces the use of full continuous stellar population synthesis
(SPS) models for redshift estimation, tailored to the multi-level calibration
and inference that will be needed for the upcoming generation of wide-deep
photometric surveys (Newman & Gruen, 2022).

Parameter | Description
---|---
| _hyper-parameters_
$\omega_{b}$ | Zero-point offset (for $b$’th band, relative to the $i$ band)
$\gamma_{b}$ | Additional flux uncertainty contribution for the $b$’th band (fraction of the total model flux)
$\alpha_{j}$ | Additional flux for the $j$’th emission line (fraction of the line strength in FSPS)
$\beta_{j}$ | Additional flux uncertainty contribution for the $j$’th emission line (fraction of the line strength in FSPS)
$\boldsymbol{\kappa}$ | hyper-parameters describing the galaxy population model
| _Latent parameters_
$\boldsymbol{\varphi}$ | Stellar population parameters describing the rest-frame spectrum (per galaxy)
$z$ | Redshift (per galaxy)
| _Derived parameters_
$\ell_{\mathrm{SPS}}(\boldsymbol{\varphi})$ | Flux energy density predicted by FSPS (per galaxy)
$\ell_{j}(\boldsymbol{\varphi})$ | Amplitude of emission line (Dirac delta model) in FSPS (per galaxy)
$L_{jb}(\boldsymbol{\varphi},z)$ | Flux of the $j$’th emission line in the $b$’th band (per galaxy), i.e., $\ell_{j}$ inserted in Eq. (2.3)
$F_{b}(\boldsymbol{\varphi},z,\boldsymbol{\alpha})$ | Total model flux in the $b$’th band, with all contributions (per galaxy), i.e., Eq. (1) inserted in Eq. (2.3)
$\Sigma_{b}(\boldsymbol{\varphi},z,\boldsymbol{\alpha},\gamma_{b},\boldsymbol{\gamma})$ | Additional flux uncertainty in the $b$’th band, with all contributions (per galaxy), see Eq. (5)
| _Data_
$\hat{F}_{b}$ | Measured flux in the $b$’th band (per galaxy)
$\sigma_{b}$ | Flux measurement uncertainty in the $b$th band (per galaxy)

Table 1: Notation for all model parameters. Vectors are indicated with bold
symbols. Our total log likelihood and graphical model presented below include
an extra index $i$ per galaxy, which we have omitted in this table and all
other equations in the paper for simplicity.

###  2.1 SED model

SPS provides a powerful way to model the SED of a galaxy. It exploits
fundamental principles to generate and sum the contributions of ‘simple’
stellar populations (as observed in, e.g., star clusters) and apply the
effects of additional complications such as dust and nebular emission in order
to create galaxy SEDs within a forward modeling framework. This step involves
a set of specific modeling choices for these components, for example an
explicit star formation history (SFH). We use the Flexible Stellar Population
Synthesis (FSPS, Conroy et al. 2009a, b; Conroy & Gunn 2010) code, accessed
through the python-FSPS binding (Foreman-Mackey et al., 2014).

For the SED model we consider a variation of Prospector-$\alpha$, which was
first introduced in Leja et al. (2017), revised in Leja et al. (2019); Johnson
et al. (2021), and has been successfully used to constrain Bayesian models for
the galaxy population Leja et al. (2020, 2021); Nagaraj et al. (2022); Whitler
et al. (2022). While these works fixed the redshift to the value measured with
spectroscopy (spec-$z$), we relax this assumption and treat it as a free
parameter. In our variation of Prospector-$\alpha$, a galaxy SED is
(deterministically) defined by 15 parameters:

  * •

8 parameters describing the SFH: 1 parameter for the total stellar mass formed
$\log M_{\star}$ [M⊙] (integral of the SFH); 1 parameter for the stellar
metallicity $\log(Z_{\star}/Z_{\odot})$, assumed to be the same for all stars
in the galaxy; and 6 parameters for the relative age SFH bins (ratios $r_{i}$
of the piece-wise SFR in adjacent temporal bins), where the seven time bins
are spaced following Leja et al. (2020).

  * •

3 parameters for the dust attenuation, following the model of Charlot & Fall
(2000) (see Leja et al. 2017 for details) with birth-cloud $\tau_{1}$, and
diffuse attenuation $\tau_{2}$ with a power law (of index $n$) from Calzetti
et al. (2000).

  * •

a gas-phase metallicity parameter $\log(Z_{\mathrm{gas}}/Z_{\odot})$ for the
nebular emission (decoupled from the stellar metallicity) modeled after the
grids of Byler et al. (2017).

  * •

two parameters, $\log{\rm f}_{\rm AGN}$ and $\tau_{\rm AGN}$, for the AGN
torus emission model of Nenkova et al. (2008).

  * •

redshift, $z$

The SED includes dust heating from stars via energy balance, via a dust SED of
fixed shape (Draine & Li, 2007).

###  2.2 Emission lines

Prospector-$\alpha$ includes a nebular emission model where the gas is ionized
by the same stars synthesized in the SED (Byler et al., 2017; Byler, 2018).
There is a large variability of the line strengths, as demonstrated by Byler
et al. (2017). In practice, previous studies have required the flexible
addition of emission lines on top of SED templates in order to obtain accurate
photo-$z$’s (Ilbert et al., 2006, 2008; Brammer et al., 2008; Alarcon et al.,
2021), so we adopt a similar approach.

We model offsets and uncertainty in the strength of the emission lines by
encoding them into parametric bias and variance parameters. For a set of SPS
parameters $\boldsymbol{\varphi}$, the total energy density as a function of
wavelength, $\ell$, consists of the base FSPS prediction
$\ell_{\mathrm{SPS}}$, with an extra additive contribution from emission
lines. In practice, in FSPS, emission lines are modeled as a delta functions,
integrated in band-passes, and added to the model photometry. We can therefore
write the SED as

$\displaystyle\ell(\lambda;\boldsymbol{\varphi},\boldsymbol{\alpha})\ =\
\ell_{\mathrm{SPS}}(\lambda;\boldsymbol{\varphi})+\sum_{j=1}^{\mathrm{N}_{\mathrm{lines}}}\alpha_{j}\ell_{j}(\boldsymbol{\varphi})\delta^{\mathrm{D}}(\lambda-\lambda_{j})\quad\
$ (1)

where $\ell_{j}$ is the amplitude of the $j$’th line, $\lambda_{j}$ the rest-
frame wavelength of the line, and $\delta^{\mathrm{D}}(\cdot)$ is the Dirac
delta function. $\alpha_{j}$ the hyper parameter for the additional
contribution from line $j$. Because line $\ell_{j}$ is already included in
$\ell_{\mathrm{SPS}}$ with a default weight, its total contribution in $\ell$
is $1+\alpha_{j}$.

###  2.3 Synthetic photometry

We can obtain total model fluxes in the $b$’th band by integrating the SED
through the band-pass filter $W_{b}(\lambda)$, such that

$\displaystyle F_{b}(\boldsymbol{\varphi},z,\boldsymbol{\alpha})=$ (2)
$\displaystyle\quad\quad\frac{(1+z)^{-1}}{4\pi
d_{L}^{2}(z)}\int_{0}^{\infty}\ell\Bigl{(}\frac{\lambda}{1+z};\boldsymbol{\varphi},\boldsymbol{\alpha}\Bigr{)}\,e^{-{\tau}(z,\lambda)}W_{b}(\lambda)d\lambda,\quad\quad$

where $d_{L}(z)$ the luminosity distance, and ${\tau}(z,\lambda)$ the
effective optical depth of the inter-galactic medium, calculated with the
model of Madau (1995).

###  2.4 Population model

So far we have a mechanism for using SPS to generate the photometry of a
single galaxy given its intrinsic properties, $\boldsymbol{\varphi}$, and
redshift, $z$. We now describe a formalism for generating a sample of galaxies
with a probability distribution
$p(\boldsymbol{\varphi},z|\boldsymbol{\kappa})$, where $\boldsymbol{\kappa}$
are the hyper-parameters which describe the population. These could in
principle be inferred as well, but this typically requires a careful treatment
of selection effects (i.e., how galaxies are selected into the sample at
hand), which must be included at the population level. We show how to treat
selection effects in a companion paper, Alsing et al. (2022). For our
demonstration with COSMOS2020 in this paper, we focus on redshift inference
for individual galaxies. Due to the constraining power of the data, the
population prior will have very little impact on the galaxy posterior
distributions, especially the redshifts. Therefore, we are able to bypass an
explicit treatment of selection effects, and also adopt a fairly uninformative
population model.

We adopt minor variations on the fiducial Prospector-$\alpha$ prior (see Leja
et al. 2019; Johnson et al. 2021), summarized in Table 2. It factorizes as

$\displaystyle
p(\boldsymbol{\varphi},z|\boldsymbol{\kappa})=p(\boldsymbol{\varphi}|z,\boldsymbol{\kappa})\
p(z)$ (3) $\displaystyle\quad\ \ =\ p(\mathrm{mass})\times
p(\mathrm{SFH})\times p(\mathrm{dust}|\boldsymbol{\kappa})$
$\displaystyle\quad\quad\times p(\mathrm{stellar\ metallicity})\times
p(\mathrm{AGN})\times p(z)$ $\displaystyle\quad\quad\times
p(\mathrm{gas\mathchar 45\relax phase\ metallicity}|\mathrm{stellar\
metallicity,\boldsymbol{\kappa}})\,.$

Most parameters have simple uniform or log-uniform priors, except in the cases
discussed below.

Table 2: Summary of parameters and priors describing the SPS model. $\mathcal{N}_{c}(l,s,m,M)$ refers to a truncated normal distribution of location $l$, scale $s$, and in the range $[m,M]$. Our fiducial values for the parameters below are $\mu_{n}=-0.095+0.111\,\tau_{2}-0.0066\,\tau_{2}^{2}$, $\sigma_{n}=0.4$, $\sigma_{Z}=3.0$, $(\mu_{1},\sigma_{1})=(1.0,0.3)$, $(\mu_{2},\sigma_{2})=(0.3,1.0)$, following Leja et al. (2019).  Parameter | Prior bounds | Prior
---|---|---
$\log M_{\star}$ [M⊙] | $[7,13]$ | Uniform($7$, $13$)
$\log(Z_{\star}/Z_{\odot})$ | $[-1.98,0.19]$ | Uniform($-1.98,0.19$)
$\\{\log r_{i}\\}_{i=1,\cdots,6}$ | $[-5,5]$ | Stu($0,0.3,2$)
$\tau_{2}$ | $[0,4]$ | $\mathcal{N}_{c}(\mu_{2},\sigma_{2},0,4)$
$\tau_{1}/\tau_{2}$ | $[0,2]$ | $\mathcal{N}_{c}(\mu_{1},\sigma_{1},0,2)$
$n$ | $[-1,0.4]$ | $\mathcal{N}_{c}(\mu_{n},\sigma_{n},-1,0.4)$
$\log{\rm f}_{\rm AGN}$ | $[10^{-5},150]$ | LogUniform($10^{-5},150$)
$\tau_{\rm AGN}$ | $[-2,0.5]$ | LogUniform($5,150$)
$\log(Z_{\mathrm{gas}}/Z_{\odot})$ | $[-2,0.5]$ | $\mathcal{N}_{c}(\log(Z_{\star}/Z_{\odot}),\sigma_{Z},-2,0.5)$
$z$ | $[0,2.5]$ | Uniform($0,2.5$)

For the gas-phase metallicity we take a normal222Since we apply normal priors
to parameters only defined in finite ranges, a renormalization is necessary,
giving rise to truncated normals, as described in Table 2. For simplicity we
omit this detail throughout this section. prior of mean at the
$\mathrm{stellar\ metallicity}$ and standard deviation $\sigma_{Z}=3.0$.
Despite the fact that they should track each other approximately, it is more
robust to treat the gas-phase metallicity as a nuisance parameter (Leja et
al., 2020). We take a normal prior on the diffuse dust components $\tau_{2}$
with mean equal to $\mu_{2}=0.3$ and standard deviation $\sigma_{2}=1.0$. For
the birth cloud component $\tau_{1}$ we take a normal prior on the ratio
$r=\tau_{1}/\tau_{2}$, with mean equal to $\mu_{1}=1$ and standard deviation
$\sigma_{1}=0.3$. The index of the dust attenuation law (for the diffuse
component) is assumed to vary as a function of the total dust attenuation,
with mean given by:
$\langle\delta\rangle=\mu_{n}=-0.095+0.111\,\tau_{2}-0.0066\,\tau_{2}^{2},$
where $\delta$ is the (negative) offset from the index of the Calzetti
attenuation curve (Calzetti et al., 2000). We take a normal prior on $\delta$,
with mean $\mu_{n}$ given above and standard deviation $\sigma_{\delta}=0.4$.
This is a simple average of the results of Leja et al. (2019), but could be
updated with a more complex prior, such as that of Nagaraj et al. (2022).

Finally, we collect the following hyper-parameters in a vector
$\boldsymbol{\kappa}=(\mu_{r},\sigma_{r},\sigma_{n},\sigma_{\delta},\sigma_{Z})$.
We will infer them jointly with the other components of the model, in order to
demonstrate how one could in principle learn about the underlying physics with
this type of forward modeling framework. However, we expect little sensitivity
to these hyper-parameters for the dataset we use. A more detailed population
prior, encoding more of the known relations from galaxy formation and
evolution, was developed by Alsing et al. (2022), and successfully reproduced
galaxy redshift distributions given broad-band photometry.

###  2.5 Noise model

The final stage of the generative model is to simulate observed photometry
$\hat{F}_{b}$ from the model flux $F_{b}$, by adding noise. The noise model
can be fully characterized using a likelihood function. We use a scaled and
translated Student’s-t distribution333The density of the scaled and translated
Student’s-t distribution is
$\displaystyle\mathrm{Stu}(x;l,s,\nu)=\frac{\Gamma(\frac{\nu+1}{2})}{\sqrt{\nu\pi}\Gamma(\frac{\nu}{2})}\left(1+\frac{t^{2}}{\nu}\right)^{-\frac{\nu+1}{2}}$
with $t=(x-l)/s$ the residuals and $\Gamma(\cdot)$ the Gamma function. $l$ and
$s$ are often referred to as location and scale parameters, and $\nu$ the
number of degrees of freedom. with two degrees of freedom, which is more
robust to outliers than a normal distribution, since it has heavier tails:

$\displaystyle p(\hat{F}_{b}|F_{b},\sigma_{b},\Sigma_{b},\omega_{b})=$ (4)
$\displaystyle\quad\quad\mathrm{Stu}\bigl{(}\hat{F}_{b}\ ;\ \omega_{b}F_{b}\
,\ \sigma_{b}^{2}+(\omega_{b}\Sigma_{b})^{2}\ ,\ 2\bigr{)}.$

Aside from the observational uncertainty $\sigma_{b}$, which is given for each
object (typically measured from the images), we introduce two additional
terms: a zero-point offset $\omega_{b}$, and the model uncertainty
$\Sigma_{b}$ added in quadrature. $\omega_{b}$ appears in front of the
uncertainty because the latter will involve the model $F_{b}$, which must be
rescaled by $\omega_{b}$ (although there is some freedom in choosing how to
define these terms self-consistently).

We write the model uncertainty as the sum (in quadrature) of a fraction of the
total model flux and a fraction of the total line fluxes,

$\displaystyle\Sigma_{b}^{2}(\boldsymbol{\varphi},z,\boldsymbol{\alpha},\gamma_{b},\boldsymbol{\gamma})=$
(5)
$\displaystyle\quad\left(\gamma_{b}F_{b}(\boldsymbol{\varphi},z)\right)^{2}+\sum_{j=1}^{\mathrm{N}_{\mathrm{lines}}}\left(\beta_{j}(1+\alpha_{j})L_{jb}(\boldsymbol{\varphi},z)\right)^{2},$

where $L_{jb}$ is the photometric flux of the $j$’th line in the $b$’th band,
i.e., the result of applying Eq. (2.3) to
$\ell_{j}\delta^{\mathrm{D}}(\lambda-\lambda_{j})$. The first term operates at
the level of the band (controlled by a hyper-parameter $\gamma_{b}$) and is
analogous to the magnitude uncertainty floor often employed in photo-$z$
techniques. The second term is similar but adds a contribution from each line.
As before, the hyper-parameter $\beta_{j}$ can be interpreted in terms of a
fraction of the total emission line flux (which is now $1+\alpha_{j}$). In
summary, this allows us to separate the extra uncertainty $\Sigma_{b}$
preferred by the data into contributions from calibration uncertainty (via
$\gamma_{b}$) and SED uncertainty (via $\beta_{j}$).

Figure 1: Graphical representation of the main parameters involved of our
model. Circles are inferred random variables, shaded circles are observed, and
dots indicate random variables that are deterministic in their inputs.

##  3 Inference formalism

We now assemble the components presented in the previous section into a
hierarchical model for generating and analysing photometric data. We describe
an inference methodology and also computational acceleration using machine
learning.

###  3.1 Hierarchical model

A summary of the parameters of our model is provided in Table 1. A graphical
representation is provided in Fig. 1. Each galaxy has 15 intrinsic SPS
parameters (including redshift), as described in Sec. 2.1, and in Secs 2.2 to
2.5 we have introduced five sets of hyper-parameters, aiming to add some
targeted sources of flexibility in the model. Furthermore, a hierarchical
modeling approach has the ability to fix subsets of parameters, and to observe
the effect on the others, offering additional robustness diagnostics.

The full posterior distribution of our model (adding a subscript $i$ for
galaxies) is

$p\bigl{(}\underbrace{\boldsymbol{\alpha},\boldsymbol{\beta},\boldsymbol{\omega},\boldsymbol{\gamma},\boldsymbol{\kappa}}_{\mathrm{hyper}},\underbrace{\\{\boldsymbol{\varphi}_{i},z_{i}\\}}_{\mathrm{latent}}\bigr{)}=\underbrace{p(\boldsymbol{\alpha},\boldsymbol{\beta},\boldsymbol{\omega},\boldsymbol{\gamma},\boldsymbol{\kappa})}_{\mathrm{global\
prior}}\prod_{i=1}^{\mathrm{N}_{\mathrm{obj}}}\underbrace{p(\boldsymbol{\varphi}_{i},z_{i}|\boldsymbol{\kappa})}_{\mathrm{population\
model}}\prod_{b=1}^{\mathrm{N}_{\mathrm{bands}}}\underbrace{p\Bigl{(}\hat{F}_{ib}\Bigl{|}\underbrace{F_{ib}(\boldsymbol{\varphi}_{i},z_{i},\boldsymbol{\alpha})}_{\mathrm{SED\
model}},\sigma_{ib},\Sigma_{ib}^{2}(\boldsymbol{\varphi}_{i},z,\boldsymbol{\alpha},\gamma_{b},\boldsymbol{\gamma}),\omega_{b}\Bigr{)}}_{\mathrm{noise\
model}}\,.$ (6)

###  3.2 Priors

In our approach it is straightforward to encode prior knowledge or constraints
on corrections such as zero-points or emission line contributions as the prior
on hyper-parameters appears explicitly as
$p(\boldsymbol{\alpha},\boldsymbol{\beta},\boldsymbol{\omega},\boldsymbol{\gamma},\boldsymbol{\kappa})$.
The priors we use are given in Table 3.

Table 3: Priors on hyper-parameters. The numbers in parentheses indicate the location and scale of the distribution. Hyper-parameter | Prior
---|---
$\\{\alpha_{j}\\}_{j=1,\cdots,\mathrm{N}_{\mathrm{lines}}}$ | Laplace(0, 10)
$\\{\beta_{j}\\}_{j=1,\cdots,\mathrm{N}_{\mathrm{lines}}}$ | Normal(0, 1)
$\\{\omega_{b}\\}_{b=1,\cdots,\mathrm{N}_{\mathrm{bands}}}$ | Normal(1, 1)
$\\{\gamma_{b}\\}_{b=1,\cdots,\mathrm{N}_{\mathrm{bands}}}$ | Normal(0, 1)

We find that the Laplace prior for each $\alpha_{j}$ is preferable to a normal
prior since it promotes sparsity of the coefficients, i.e., the coefficients
should be zero unless they significantly improve the model fits. Without this
regularization there is a risk that line contributions would mimic other
narrow features in the SED, such as absorption lines, which should instead be
associated with a separate set of stellar evolution physics.

###  3.3 Optimization strategy

Our inference strategy is to optimize the posterior distribution with respect
to all parameters. This is because the mixture of broad bands (typically
deeper) and noisier intermediate or narrow bands makes it possible to pin down
SED features accurately and unequivocally. As a result, it is possible to
derive estimates of the hyper- and latent parameters through optimization,
without the need for more involved inference strategies, such as MCMC. We
adopt the Adam optimizer with a learning rate of $10^{-3}$.

Challenges that need to be overcome at this stage involve the computational
needs of SPS calls, and the fact that they are not differentiable analytically
(although differentiable SPS models are now being developed, Alarcon et al.
2022; Hearin et al. 2021). We address these issues with neural emulators of
the model predictions, which are both $\sim 10^{4}$ faster than native SPS
calls and differentiable (Alsing et al., 2020). We return to the details of
our neural emulators in the next section.

In principle, one could determine all parameters in a single run. But for the
purpose of this demonstration study we adopt a two-stage approach (e.g.,
Laigle et al. 2016; Weaver et al. 2022) which is commonly used and leverages
the availability of spectroscopic redshifts. We first determine the hyper-
parameters with redshifts fixed at the values acquired via spectroscopy
(optimizing all SPS parameters other than the redshift). A second run is then
performed by fixing the hyper-parameters at the inferred values, and only
optimizing the latent SPS parameters (now including the redshift). The
resulting optimized redshifts are denoted $\hat{z}^{\mathrm{MAP}}$ (omitting
the galaxy index).

Finally, we employ bijectors to map the original parameters to new variables
which have zero-mean, unit-variance normal priors. These are easier to work
with than the original variables since they do not have boundaries, which
facilitates optimization (and sampling). The parameter bounds are shown
explicitly in Table 2.

###  3.4 Emulators

As previously mentioned, one drawback associated with SPS models is the
computational cost of each model call, which is the result of the many
intensive numerical steps involved in predicting flux densities and photometry
from stellar isochrones, initial mass function, and the other ingredients of
SPS. While this has been overcome with MCMC in the past and applied to medium-
sized ($\sim 10^{5}$ objects) data sets (e.g., Leja et al. 2019), photometric
galaxy surveys are approaching of order a billion objects (Abbott et al.,
2022; The LSST Dark Energy Science Collaboration et al., 2018), so brute-force
inference is not practical. Hyper-parameters or model calibration would make
the computational load even heavier. We now discuss how we employ
emulators444Note that emulators are not the only approach to accelerating
inference. For example, Hahn & Melchior (2022) and Ramachandra et al. (2021)
showed how to approximate posterior distributions directly for the estimation
of SPS parameters or redshifts from broad-band photometry. of fluxes and line
emission in order to overcome this challenge.

####  3.4.1 SPS fluxes

For each band, we use the methodology described in Alsing et al. (2020) to
train an emulator to predict the photometric flux, i.e., $\ell_{\mathrm{SPS}}$
fed into Eq. (2.3). While in principle one could train emulators to predict
the $\ell_{\mathrm{SPS}}(\lambda)$, it is advantageous to emulate photometric
fluxes directly in order to avoid evaluating Eq. (2.3). The emulator for each
band is a fully-connected 4-layer neural network (with 128 neurons in each
layer). We train with staggered learning rates ($10^{-3}$, $10^{-4}$,
$10^{-5}$, $10^{-6}$) and batch sizes (1000, 10000, 50000, all objects) to
improve convergence. We train on 4 million samples from the prior, setting 10%
of the training data aside for validation, and stopping the training when the
loss does not improve for 20 epochs.

####  3.4.2 Emission line fluxes

Our model photometry has additional contributions from emission lines, so for
any parameter vector $(\boldsymbol{\varphi},z)$ we need to calculate them in
addition to the base SPS predictions. This extra computational load can again
be alleviated with emulators. However, in this case building emulators for
individual emission lines and all bands would be prohibitive. Instead, we
adopt a different strategy. We train an emulator to model the fluxes
$\\{\ell_{j}\\}$ on the rest-frame wavelength grid $\\{\lambda_{j}\\}$,
following the procedure of Alsing et al. (2020). In our setting the size of
the wavelength grid is the number of emission lines in FSPS (all $128$ lines
in Byler et al. 2017), so the accuracy of the emulator is easier to keep under
control555For reference, a standard FSPS spectrum with the default MILES
spectral library has $\sim 6\times 10^{3}$ wavelength points.. To predict the
photometric flux of an emission line, we need to implement Eq. (2.3). Because
the band-passes $\\{W_{b}(\lambda)\\}$ are in tabular form, it would be
preferable to work with an interpolated or continuous approximation. We find
that fitting the sum of three generalized normal distributions to each $W_{b}$
offers a simple solution, since it is straightforward to fit to the tabulated
data without any tuning. This approach delivers accurate and fast results for
the final emission line fluxes.

####  3.4.3 Emulator accuracy

When comparing the initial and the emulator-reconstructed fluxes (computed
from fluxes), we find that they are accurate at a level better than 1% (99th
percentile of the magnitude differences in the range spanned by the COSMOS2020
data below). This will be comfortably covered by $\Sigma_{b}$, for example
thanks to the $\gamma_{b}$ hyper parameter being greater than 0.01. The
accuracy of emulation of emission line fluxes is typically below $0.1\%$ (99th
percentile) so the errors in this emulation can safely be ignored.

Figure 2: Photo-$z$ metrics calculated on our sample for our hierarchical
model compared with the COSMOS2020 released LePhare and Eazy photo-$z$’s. We
compute the median of
$\Delta_{z}=(z_{\mathrm{spec}}-\hat{z}^{\mathrm{MAP}})/(1+z_{\mathrm{spec}})$,
as well as $\sigma_{\mathrm{MAD}}=1.48\times\mathrm{median}(|\Delta_{z}|)$,
and the outlier fraction as the fraction of objects with $|\Delta_{z}|>0.15$.
Figure 3: Differences between the maximum a-posteriori redshift estimates and
the spectroscopic redshifts (divided by $1+z$) for the new model compared with
the COSMOS2020 released LePhare and Eazy redshifts. This highlights the areas
where they achieve similar performance in terms of scatter and outliers.
Splitting between the three sources of spectroscopic redshifts also reveals
the differences in photo-$z$ performance for samples with different selection,
as described in Sec. 4.  Figure 4: Maximum a-posteriori redshift estimates
versus spectroscopic redshifts for the three spectroscopic sources comprising
the data we analyze. Differences are best highlighted in Fig. 3, whereas these
scatter plots with a logarithmic color bar better highlight the overall
distributions and the outliers (which are rare).

##  4 Data

We now construct a dataset that demonstrates that our framework is able to
produce accurate photometric redshift inferences on a benchmark sample. This
consists of selecting broad- and narrow-band photometry in the COSMOS field,
and cross-matching with a selection of spectroscopic datasets to enable
validation of inferred redshifts.

###  4.1 Photometry

Weaver et al. (2022) presented COSMOS2020, the latest compilation of imaging
in the COSMOS field. Two sets of photometric measurements were performed on
these data: Classic, and Farmer. For our case study we choose Farmer, which is
based on the model-fitting software The
Tractor666https://github.com/dstndstn/tractor. The photo-$z$’s published with
this dataset were derived using Eazy (Brammer et al., 2008), as well as
LePhare (Ilbert et al., 2006, 2008) using slightly different subsets of bands.
While we will compare our results to both of these methods, we adopt the same
subset of selected bands as Eazy; its SED templates are also derived using
FSPS, thus providing a more direct comparison. This selection of bands exclude
the Subaru Suprime-Cam broad-bands (shallower than other filters at similar
wavelengths) and the GALEX bands (shallow and with broad PSFs) (Weaver et al.,
2022). We apply the recommended ‘combined’ mask, which retains the deepest
regions with the greatest number of available bands, and removes areas
corrupted by bright stars and artifacts in all the relevant bands.

This photometric dataset is prepared for further analysis using the code
released with the COSMOS2020 data, which applies various flux corrections
(including Galactic extinction) and unit conversions.

###  4.2 Spectroscopic redshifts

We compile the positions and spectroscopic redshifts from three large
campaigns covering the COSMOS field: zCosmos-bright (Lilly et al., 2007),
DEIMOS (Hasinger et al., 2018)], and C3R2 (Masters et al., 2017, 2019;
Stanford et al., 2021).

zCosmos-bright (Lilly et al., 2007) is a highly representative sample of 20k
bright ($i*\leq 22.5$) galaxies. DEIMOS (Hasinger et al., 2018) is a sample of
$\sim 4\times 10^{3}$ objects selected from a variety of input catalogs based
on multi-wavelength observations in the field, and thus represents a diverse
selection function. Finally, the Complete Calibration of the Color-Redshift
Relation (C3R2, Masters et al. 2017, 2019; Stanford et al. 2021) is a multi-
institution, multi-instrument survey targeting faint galaxies ($i\sim 24.5$)
populating regions of color space under-sampled by other surveys, in order to
provide validation data for upcoming weak lensing surveys. We use data from
the Data Releases 1, 2, and 3 ($\sim 4\times 10^{3}$ objects).

For all three samples, we keep objects with quality flags 3 and 4. We perform
a spatial cross-match between the positions of these objects in the
spectroscopic and photometric catalogs. Specifically, we keep the nearest
unique photometric source within a radius of 0.6 arcseconds around each
spectroscopic object, as recommended by Weaver et al. (2022). We also remove
objects classified as stars by Eazy or LePhare, as well as the few sources
with spectroscopic redshift $>2.5$.

We are left with $12,473$ objects: 6920, 3168, and 2043 for zCosmos-bright,
DEIMOS, and C3R2, respectively. This is the dataset on which we carry out our
demonstration study.

This dataset combines three samples and so has a complicated selection
function (i.e., the probability, conditional on its intrinsic properties, that
a galaxy appears in the sample). The selection function must be incorporated
into any population level analysis, as is treated explicitly in the inference
of redshift distributions presented by Alsing et al. (2022). But if the focus
is on inference of object-level parameters (i.e., individual galaxy redshifts)
then selection effects can be ignored if the data on detected objects is
sufficiently constraining: in the limit that the photometry determines the
redshift (and other parameters) perfectly, the posterior is a delta function
and hence completely decoupled from the population. While not strictly
satisfied here - the high-quality photometry and broad wavelength coverage do
yield finite parameter uncertainties - this is an excellent approximation. We
hence do not consider selection effects further.

Figure 5: Hyper-parameters obtained with the optimization procedure of Sec.
3.3. The black stars correspond to our best model, with all parameters
optimized, while the green diamonds correspond to an additional model (used
for robustness checks, see text for details) fixing the emission line offsets
$\alpha_{j}=0\ \forall j$ and error floors $\gamma_{b}=0.01\ \forall b$.

##  5 Results

We now evaluate the performance of our inference framework applied to the
cross-matched dataset described in the previous section. In what follows,
unless explicitly stated, the hyper-parameters of the population model
($\boldsymbol{\kappa}$) are not optimized, and are set to the fiducial values
defined previously.

In practice, to avoid degeneracies with other parameters of the model, we
anchor the zero-point of the HSC $i$ band to unity, and define the others
relative to this band.

###  5.1 Emission lines

We perform an initial fit to the data with our model with no contribution from
lines, i.e., setting $\alpha_{j}=\beta_{j}=0\ \forall j$. Thus, we optimize
the zero-point offsets and uncertainty hyper-parameters
$\\{\omega_{b},\gamma_{b}\\}_{b=1,\cdots,\mathrm{N}_{\mathrm{bands}}}$, as
well as the SPS parameters of each galaxy, with the procedure outlined in Sec.
3.3. This consists of a run with redshifts fixed at the spectroscopic
redshifts to estimate the hyper-parameters, and a second run now relaxing
redshifts but fixing the hyper-parameters. The resulting maximum a-posteriori
redshift estimates, $\hat{z}^{\mathrm{MAP}}$, are poor compared to the
released COSMOS2020 photo-$z$’s. We also observe sharp features in the
redshift residuals $\hat{z}^{\mathrm{MAP}}-$spec-$z$ and in the flux residuals
(data minus model at the best-fit parameters) binned in spec-$z$. Many of
these features correspond to strong emission lines, such as [O II],
demonstrating that their luminosities are not correctly calibrated for our
data set. Re-calibrating them is straightforward within our framework.

We must first determine a suitable subset of emission lines (out of the $128$
implemented in FSPS by Byler et al. 2017) to include in the model. To achieve
this we take the best-fit SEDs of our first run described above (at
$z$=spec-$z$) and calculate the contribution of each emission line to the
total photometric flux, i.e., the flux ratio $L_{jb}/F_{b}$, for all objects,
lines, and bands. We then take the average of these values over objects and
bands. This procedure yields the lines which make the strongest contributions
to the data at hand. This set includes most of the lines with intrinsic large
variance according to Byler et al. (2017). Importantly, this procedure also
filters out lines that are outside our wavelength coverage. We are left with
$44$ lines. Note that nearby emission lines are not necessarily all well-
resolved by the data, i.e., each may only unequivocally contribute to a few
objects. This could be improved with more informative priors on the lines or
on line-selection; however we do not expect this to have a significant impact
on our redshift results.

In total, following this line-selection procedure, we have
$2\times\mathrm{N}_{\mathrm{bands}}+2\times\mathrm{N}_{\mathrm{lines}}=140$
hyper-parameters (excluding $\boldsymbol{\kappa}$), and
$15\times\mathrm{N}_{\mathrm{objects}}\sim 2\times 10^{5}$ ‘latent’ SPS
parameters. We now optimize all these parameters and extract maximum
a-posteriori redshift estimates $\hat{z}^{\mathrm{MAP}}$.

###  5.2 Photo-$z$’s

We evaluate the performance of this model using some commonly-used metrics.
For every object in our catalog, a spectroscopic redshift $z_{\mathrm{spec}}$
and a redshift estimate $\hat{z}^{\mathrm{MAP}}$ are available. The main
quantity of interest is the difference between the two, divided by the
standard $1+z$ factor accounting for the expected scaling of the quality of
the estimates due to the effect of redshift on wavelength,

$\displaystyle\Delta_{z}=\frac{z_{\mathrm{spec}}-\hat{z}^{\mathrm{MAP}}}{1+z_{\mathrm{spec}}}\,.$
(7)

We compute the mean, median, and standard deviation of $\Delta_{z}$ over our
sample of galaxies. For the standard deviation, an estimator more robust to
outliers than the sample variance is $\Sigma_{\mathrm{NMAD}}$, the median
absolute deviation (the median of $|\Delta_{z}|$) multiplied by $1.48$.
Finally, the outlier fraction is defined as the percentage of objects with
$|\Delta_{z}|>0.15$.

Fig. 2 shows these metrics, calculated in redshift and magnitude bins (in the
reference HSC $i$ band, which is not affected by zero-points), for our
framework as well as the publicly-released COSMOS2020 LePhare and Eazy
redshifts.

Overall, our approach yields comparable levels of bias (as measured by the
median of $\Delta_{z}$), fewer outliers, and a slightly larger level of
scatter (as measured by $\sigma_{\mathrm{NMAD}})$. The comparison with LePhare
is not straightforward since it uses a slightly different set of bands. On the
other hand, we use the same set of bands as Eazy, also based on templates
derived from a grid of parameters using a SPS model. Eazy includes a
magnitude-redshift prior as well as a template error function (Brammer et al.,
2008). These differences may be responsible for the small improvement in
$\sigma_{\mathrm{NMAD}}$. Note that we have not attempted to tune or improve
our model in order to carry out redshift inferences on these data, beyond the
procedure described above.

While these photo-$z$ metrics are informative summaries of the results, it is
also interesting to examine the distributions of redshift estimates
themselves. The residuals $\Delta_{z}$ are shown in Fig. 3, while Fig. 4 also
shows a conventional comparison of the estimates and the spectroscopic
redshifts. These confirm the results captured by the metrics above. This also
highlights that the redshift estimates are noisier for the DEIMOS and C3R2
samples because these surveys targeted populations typically under-represented
in other spectroscopic campaigns (fainter, bluer, higher redshift). These
populations also present greater challenges in terms of accurate photo-$z$
estimation.

Figure 6: Flux residuals (black), as defined in Eq. (8), compared with the
likelihood function of Eq. (4) in blue, showing a satisfactory agreement
between the data and the model predictions.

###  5.3 Hyper-parameters

In this section, we investigate the results of our optimization procedure for
the model’s hyper-parameters, and their implications for the robustness of the
data-model. The hyper-parameters resulting from the optimization are shown in
Fig. 5. We also show the result of a run with the noise and emission line
offsets set to fiducial values $\gamma_{b}=0.01\ \forall b$ , $\alpha_{i}=0\
\forall i$. Throughout this section we will also comment on robustness checks
obtained by running variations of these models.

####  5.3.1 Zero-point offsets

The inferred zero-points $\\{\omega_{b}\\}$ are at the level of a few percent
in most bands. They are well-constrained by the data and robust to changes in
the model. The biggest changes are shifts in the intermediate bands caused by
emission line offsets $\alpha_{i}$, as shown in Fig. 5. These results are
qualitatively consistent with the Eazy zero-points released with COSMOS2020.
This is not a surprise given that we use the same set of bands and that Eazy
also involves a grid of SPS templates. Some possible explanations for the
disparities are differences in the SPS modeling, the approximation resulting
from using a grid as opposed to a continuous model, or even the fact that Eazy
calibrates the zero-points based on the whole COSMOS2020 dataset (not only
objects with spectroscopic redshifts).

####  5.3.2 Examination of photometric uncertainties

Fig. 5 shows that the inferred values of flux uncertainties $\gamma_{b}$ are
at the level of a few percent in most bands: smaller for the broad bands
(which are typically deeper), and larger for intermediate and narrow bands
(noisier and also more sensitive to the details of the SED modeling, such as
emission lines). The $u$ and ch$2$ bands are the most noisy, consistent with
previous findings (Laigle et al., 2016; Weaver et al., 2022). Just as for
zero-points, the inferred values are mostly stable when switching other hyper-
parameters on and off. However, there is no uncertainty contributed by
emission lines in this model (bottom right panel), indicating that the band-
to-band uncertainty is preferred by the data. If we set the band-to-band
uncertainty to 1% (green diamonds), then emission lines contribute to the
uncertainty budget (via non-zero $\beta$) to compensate, as expected. Thus, we
conclude that the noise model performs well, but that emission line
uncertainty is simply not required in our full model for this dataset.

####  5.3.3 Emission line offsets

We find that including extra contributions from emission lines significantly
improves the quality of the derived photo-$z$’s. This is consistent with
previous findings Ilbert et al. (2006, 2008); Brammer et al. (2008); Alarcon
et al. (2021). Note that in Fig. 5 emission lines are ordered by their
contribution relative to the total broadband flux, following the procedure
described above. Offsets are primarily given to the strongest lines, with some
of them being almost entirely cancelled ($\alpha\approx-1$). The Laplace prior
is effective at keeping most line offsets at zero, unless they bring a
significant improvement to the fits. We also note that the interpretation of
the hyper-parameters is conditioned on the sample at hand: changing the sub-
set of COSMOS2020 galaxies analysed may change the resulting values of zero-
points and emission line offsets, for example.

###  5.4 Population model

We perform one additional run also optimizing the hyper-parameters of the
population model, $\boldsymbol{\kappa}$. We find that the data do not
significantly perturb the fiducial values originally adopted; nor do these
parameters appear degenerate with other parameters of the model. This
demonstrates the robustness of the model to the prior, and additionally shows
that varying these hyper-parameters does not lead to significantly improved
fits. However, it should be noted that we have not included selection effects
in our model. If we had found that the inferred values had significantly moved
from the fiducial values, we would have been unable to determine whether this
had been caused by incomplete modeling of the selection process (Alsing et
al., 2022). However in this setting, as we have previously argued, the
selection effects should be negligible so this concern does not arise.

###  5.5 Residuals

Finally, we check the ability of the presented model to explain the
statistical properties of the data. We examine the distribution of the flux
residuals, defined as

$\displaystyle
r_{b}=\frac{\hat{F}_{b}-\omega_{b}F_{b}}{\sqrt{\sigma_{b}^{2}+(\omega_{b}\Sigma_{b})^{2}}}\,.$
(8)

These are shown in Fig. 6. The blue lines show our likelihood function, the
Student’s-t distribution described in Eq. (4). There is good agreement between
the chosen likelihood and the residuals. We have checked that this agreement
significantly deteriorates when all the hyper-parameters are not optimized
(not shown here). The deviations in the $u$ and ch2 bands are likely to be due
to residual data systematics which are not well-described by our model.

When repeating this analysis with a normal likelihood function, we find that
the residuals are worse and provide a poor fit of the likelihood function,
despite optimizing the other components of the model. This setting also
results in worse redshift estimates. This demonstrates the presence of
outliers, and shows that accommodating them with a suitable likelihood
function can significantly improve the ability of a model to fit complicated
data.

##  6 Discussion

Our results demonstrate that the hierarchical SPS inference framework
presented here (and in Alsing et al. 2022) delivers accurate redshifts on a
benchmark dataset, therefore passing a stringent validation test. We now
discuss possible extensions which relax some of its assumptions.

The model itself could be made more sophisticated in various ways in order to
better suit specific applications and increase the accuracy of the photo-$z$’s
further. Such extensions include a stellar model, which could
straightforwardly be incorporated to allow star-galaxy separation, akin to
what is done in Eazy and LePhare with stellar templates.

In the future it will be critical to include effects which are currently
neglected but affect the data (and in turn, photometric redshifts) at
significant levels for LSST, such as airmass-dependent fluxes (Graham et al.
2018). Furthermore, extending our model to include spatial information
(Sánchez & Bernstein, 2018; Jasche & Lavaux, 2019; Alarcon et al., 2020;
Sánchez et al., 2020) could provide a powerful probe of the connection between
galaxies and dark matter, and exploit more of the information offered by
photometric surveys.

Our chosen optimization strategy, paired with the posterior distribution taken
as a loss function, minimizes the residuals between the measured and the model
fluxes. Other choices may be possible, especially when including spectroscopic
redshifts. For example, one could explicitly minimize the residuals between
photometric (inferred) and spectroscpic redshifts, although this will render
the results more sensitive to selection effects in the data.

We neglected the uncertainties on hyper-parameters and SPS parameters. This is
because the extensive wavelength coverage of the 26-band COSMOS2020 data, the
small flux uncertainties, and the availability of spectroscopic redshifts
combine to yield well-constrained galaxy SEDs and in turn hyper-parameters.
Nevertheless, there are fundamental degeneracies expected between some of the
SPS parameters (Leja et al., 2017, 2019), as well as with the redshift when it
is not fixed to a spectroscopic value. Selection effects may also be
important, and can be accounted for in the posterior distributions as
described in Alsing et al. (2022). In analyzing datasets where the
uncertainties on hyper-parameters or SPS parameters (including redshift) are
not negligible (e.g., broadband-only surveys like the LSST), one can resort to
the variety of Bayesian inference techniques (i.e., MCMC, variational
inference, Laplace approximation, simulation-based inference) available for
hierarchical models. Since the impact of selection cuts and parameter
uncertainties is strongly dependent on the specific data and models being
considered, we defer further discussion and the development of effective
inference strategies to future work.

##  7 Conclusion

We have presented a hierarchical model to infer redshift and other intrinsic
galaxy properties from photometric data. This approach makes it possible to
self-consistently encode knowledge of (1) galaxy formation and evolution via a
population model, (2) astrophysics via SEDs predicted with stellar population
synthesis, (3) observational effects via a noise model and a likelihood
function. By formulating the inference within a Bayesian hierarchical model,
we are able to parameterize additional sources of bias and uncertainty in each
of these components, and to solve for them self-consistently. Thus, any manual
tuning or inversion methods can be avoided, with the added benefit that it is
possible to set informative priors on these additional terms.

Our photo-$z$’s (as measured with traditional metrics on redshift point
estimates) are competitive with the publicly released photo-$z$’s from
COSMOS2020 (Weaver et al., 2022). While this is a benchmark test on relatively
bright, low-redshift sample of spectroscopically-confirmed galaxies, these
results demonstrate the power and flexibility of our methodology, and motivate
its further development towards addressing the challenging setting presented
by stage IV cosmological surveys. We made our inference tractable using neural
emulators to accelerate model calls (Alsing et al., 2020).

In future applications of our framework, we will relax some of the modeling
assumptions made here. In particular, we demonstrate how to incorporate
selection in the companion paper (Alsing et al., 2022). Comparing inferences
from multiple SPS and population models will also provide valuable consistency
checks. In fact, refining the SED and population models directly from the
upcoming survey data is a promising avenue for directly probing the underlying
(astro)physics using large amounts of untapped statistical power. Leveraging
hybrid data sets and incorporating additional data available (e.g.,
spec-$z$’s, or even full spectra) will be instrumental for correctly
separating the different physical effects at play. This is straightforwardly
included in hierarchical models (see e.g., Alarcon et al. 2021; Nagaraj et al.
2022). Together, these attributes make the methodology presented here well-
suited for addressing the challenges of upcoming photometric surveys and the
robust scientific exploitation of their data.

Author contributions. BL: Conceptualization, methodology, software,
validation, formal analysis, writing - original draft. JA: Conceptualization,
methodology, software, validation, writing - review & editing. HVP:
Conceptualization, methodology, validation, writing - review & editing,
funding acquisition. DM: Conceptualization, methodology, validation, writing -
review & editing, funding acquisition. JL: Conceptualization, methodology,
validation, writing - review & editing.

Acknowledgements. We thank George Efstathiou for valuable input during the
course of this project, and John R. Weaver for assistance with the COSMOS2020
data. We also thank Konrad Kuijken, Hendrik Hildebrandt, Angus Wright, and
Will Hartley for useful discussions. BL is supported by the Royal Society
through a University Research Fellowship. This project has received funding
from the European Research Council (ERC) under the European Union’s Horizon
2020 research and innovation programme (grant agreement no. 101018897
CosmicExplorer). This work has also been enabled by support from the research
project grant ‘Understanding the Dynamic Universe’ funded by the Knut and
Alice Wallenberg Foundation under Dnr KAW 2018.0067. JA, HVP and DJM were
partially supported by the research project grant “Fundamental Physics from
Cosmological Surveys” funded by the Swedish Research Council (VR) under Dnr
2017-04212. The work of HVP was additionally supported by the Göran Gustafsson
Foundation for Research in Natural Sciences and Medicine. HVP and DJM
acknowledge the hospitality of the Aspen Center for Physics, which is
supported by National Science Foundation grant PHY-1607611. The participation
of HVP and DJM at the Aspen Center for Physics was supported by the Simons
Foundation.

## References

  * Abbott et al. (2022) Abbott, T., Aguena, M., Alarcon, A., et al. 2022, Physical Review D, 105, doi: 10.1103/physrevd.105.023520
  * Acquaviva et al. (2015) Acquaviva, V., Raichoor, A., & Gawiser, E. 2015, The Astrophysical Journal, 804, 8, doi: 10.1088/0004-637x/804/1/8
  * Ade et al. (2016) Ade, P. A., Aghanim, N., Arnaud, M., et al. 2016, Astronomy & Astrophysics, 594, A13
  * Alarcon et al. (2022) Alarcon, A., Hearin, A. P., Becker, M. R., & Chaves-Montero, J. 2022, Diffstar: A Fully Parametric Physical Model for Galaxy Assembly History, arXiv, doi: 10.48550/ARXIV.2205.04273
  * Alarcon et al. (2020) Alarcon, A., Sánchez, C., Bernstein, G. M., & Gaztañaga, E. 2020, Monthly Notices of the Royal Astronomical Society, 498, 2614, doi: 10.1093/mnras/staa2478
  * Alarcon et al. (2021) Alarcon, A., Gaztanaga, E., Eriksen, M., et al. 2021, Monthly Notices of the Royal Astronomical Society, 501, 6103, doi: 10.1093/mnras/staa3659
  * Alsing et al. (2022) Alsing, J., Leistedt, B., Leja, J., Peiris, H., & Efstathiou, D. M. G. 2022, The Astrophysical Journal Supplement Series
  * Alsing et al. (2020) Alsing, J., Peiris, H., Leja, J., et al. 2020, The Astrophysical Journal Supplement Series, 249, 5, doi: 10.3847/1538-4365/ab917f
  * Battisti et al. (2019) Battisti, A. J., da Cunha, E., Grasha, K., et al. 2019, The Astrophysical Journal, 882, 61, doi: 10.3847/1538-4357/ab345d
  * Benitez (2000) Benitez, N. 2000, The Astrophysical Journal, 536, 571, doi: 10.1086/308947
  * Brammer et al. (2008) Brammer, G. B., van Dokkum, P. G., & Coppi, P. 2008, The Astrophysical Journal, 686, 1503, doi: 10.1086/591786
  * Byler (2018) Byler, N. 2018, nell-byler/cloudyfsps: Initial release of cloudyFSPS, v1.0.0, Zenodo, doi: 10.5281/zenodo.1156412
  * Byler et al. (2017) Byler, N., Dalcanton, J. J., Conroy, C., & Johnson, B. D. 2017, The Astrophysical Journal, 840, 44, doi: 10.3847/1538-4357/aa6c66
  * Calzetti et al. (2000) Calzetti, D., Armus, L., Bohlin, R. C., et al. 2000, The Astrophysical Journal, 533, 682
  * Carnall et al. (2018) Carnall, A. C., McLure, R. J., Dunlop, J. S., & Davé , R. 2018, Monthly Notices of the Royal Astronomical Society, 480, 4379, doi: 10.1093/mnras/sty2169
  * Chabrier (2003) Chabrier, G. 2003, Publications of the Astronomical Society of the Pacific, 115, 763, doi: 10.1086/376392
  * Charlot & Fall (2000) Charlot, S., & Fall, S. M. 2000, The Astrophysical Journal, 539, 718, doi: 10.1086/309250
  * Collister & Lahav (2004) Collister, A. A., & Lahav, O. 2004, Publications of the Astronomical Society of the Pacific, 116, 345, doi: 10.1086/383254
  * Conroy & Gunn (2010) Conroy, C., & Gunn, J. E. 2010, The Astrophysical Journal, 712, 833, doi: 10.1088/0004-637x/712/2/833
  * Conroy et al. (2009a) Conroy, C., Gunn, J. E., & White, M. 2009a, The Astrophysical Journal, 699, 486, doi: 10.1088/0004-637x/699/1/486
  * Conroy et al. (2009b) Conroy, C., White, M., & Gunn, J. E. 2009b, The Astrophysical Journal, 708, 58, doi: 10.1088/0004-637x/708/1/58
  * Draine & Li (2007) Draine, B. T., & Li, A. 2007, The Astrophysical Journal, 657, 810, doi: 10.1086/511055
  * Foreman-Mackey et al. (2014) Foreman-Mackey, D., Sick, J., & Johnson, B. 2014, python-fsps: Python bindings to FSPS (v0.1.1), v0.1.1, Zenodo, doi: 10.5281/zenodo.12157
  * Gatti et al. (2021) Gatti, M., Giannini, G., Bernstein, G. M., et al. 2021, Monthly Notices of the Royal Astronomical Society, 510, 1223, doi: 10.1093/mnras/stab3311
  * Graham et al. (2018) Graham, M. L., Connolly, A. J., Ivezić, Ž ., et al. 2018, The Astronomical Journal, 155, 1, doi: 10.3847/1538-3881/aa99d4
  * Hahn & Melchior (2022) Hahn, C., & Melchior, P. 2022, Accelerated Bayesian SED Modeling using Amortized Neural Posterior Estimation, arXiv, doi: 10.48550/ARXIV.2203.07391
  * Hartley et al. (2020) Hartley, W. G., Chang, C., Samani, S., et al. 2020, Monthly Notices of the Royal Astronomical Society, 496, 4769, doi: 10.1093/mnras/staa1812
  * Hasinger et al. (2018) Hasinger, G., Capak, P., Salvato, M., et al. 2018, The Astrophysical Journal, 858, 77, doi: 10.3847/1538-4357/aabacf
  * Hearin et al. (2021) Hearin, A. P., Chaves-Montero, J., Alarcon, A., Becker, M. R., & Benson, A. 2021, DSPS: Differentiable Stellar Population Synthesis, arXiv, doi: 10.48550/ARXIV.2112.06830
  * Hildebrandt et al. (2020) Hildebrandt, H., Köhlinger, F., van den Busch, J. L., et al. 2020, Astronomy & Astrophysics, 633, A69, doi: 10.1051/0004-6361/201834878
  * Ilbert et al. (2006) Ilbert, O., Arnouts, S., McCracken, H. J., et al. 2006, Astronomy & Astrophysics, 457, 841, doi: 10.1051/0004-6361:20065138
  * Ilbert et al. (2008) Ilbert, O., Capak, P., Salvato, M., et al. 2008, The Astrophysical Journal, 690, 1236, doi: 10.1088/0004-637x/690/2/1236
  * Jasche & Lavaux (2019) Jasche, J., & Lavaux, G. 2019, Astronomy & Astrophysics, 625, A64, doi: 10.1051/0004-6361/201833710
  * Johnson et al. (2021) Johnson, B. D., Leja, J., Conroy, C., & Speagle, J. S. 2021, The Astrophysical Journal Supplement Series, 254, 22, doi: 10.3847/1538-4365/abef67
  * Jones & Heavens (2018) Jones, D. M., & Heavens, A. F. 2018, Monthly Notices of the Royal Astronomical Society, 483, 2487, doi: 10.1093/mnras/sty3279
  * Kind & Brunner (2013) Kind, M. C., & Brunner, R. J. 2013, Monthly Notices of the Royal Astronomical Society, 432, 1483, doi: 10.1093/mnras/stt574
  * Laigle et al. (2016) Laigle, C., McCracken, H. J., Ilbert, O., et al. 2016, The Astrophysical Journal Supplement Series, 224, 24, doi: 10.3847/0067-0049/224/2/24
  * Leistedt et al. (2019) Leistedt, B., Hogg, D. W., Wechsler, R. H., & DeRose, J. 2019, The Astrophysical Journal, 881, 80, doi: 10.3847/1538-4357/ab2d29
  * Leistedt et al. (2016) Leistedt, B., Mortlock, D. J., & Peiris, H. V. 2016, Monthly Notices of the Royal Astronomical Society, 460, 4258, doi: 10.1093/mnras/stw1304
  * Leja et al. (2017) Leja, J., Johnson, B. D., Conroy, C., van Dokkum, P. G., & Byler, N. 2017, The Astrophysical Journal, 837, 170, doi: 10.3847/1538-4357/aa5ffe
  * Leja et al. (2020) Leja, J., Speagle, J. S., Johnson, B. D., et al. 2020, The Astrophysical Journal, 893, 111, doi: 10.3847/1538-4357/ab7e27
  * Leja et al. (2019) Leja, J., Johnson, B. D., Conroy, C., et al. 2019, The Astrophysical Journal, 877, 140, doi: 10.3847/1538-4357/ab1d5a
  * Leja et al. (2021) Leja, J., Speagle, J. S., Ting, Y.-S., et al. 2021, A New Census of the $0.2<z<3.0$ Universe, Part II: The Star-Forming Sequence, arXiv, doi: 10.48550/ARXIV.2110.04314
  * Lilly et al. (2007) Lilly, S. J., Fevre, O. L., Renzini, A., et al. 2007, The Astrophysical Journal Supplement Series, 172, 70, doi: 10.1086/516589
  * Madau (1995) Madau, P. 1995, ApJ, 441, 18, doi: 10.1086/175332
  * Masters et al. (2017) Masters, D. C., Stern, D. K., Cohen, J. G., et al. 2017, The Astrophysical Journal, 841, 111, doi: 10.3847/1538-4357/aa6f08
  * Masters et al. (2019) —. 2019, The Astrophysical Journal, 877, 81, doi: 10.3847/1538-4357/ab184d
  * McQuinn & White (2013) McQuinn, M., & White, M. 2013, Monthly Notices of the Royal Astronomical Society, 433, 2857, doi: 10.1093/mnras/stt914
  * Morrison et al. (2017) Morrison, C. B., Hildebrandt, H., Schmidt, S. J., et al. 2017, Monthly Notices of the Royal Astronomical Society, 467, 3576, doi: 10.1093/mnras/stx342
  * Myles et al. (2021) Myles, J., Alarcon, A., Amon, A., et al. 2021, Monthly Notices of the Royal Astronomical Society, 505, 4249, doi: 10.1093/mnras/stab1515
  * Ménard et al. (2013) Ménard, B., Scranton, R., Schmidt, S., et al. 2013, Clustering-based redshift estimation: method and application to data, arXiv, doi: 10.48550/ARXIV.1303.4722
  * Nagaraj et al. (2022) Nagaraj, G., Forbes, J. C., Leja, J., Foreman-Mackey, D., & Hayward, C. C. 2022, The Astrophysical Journal, 932, 54, doi: 10.3847/1538-4357/ac6c80
  * Nenkova et al. (2008) Nenkova, M., Sirocky, M. M., Ivezić, Ž ., & Elitzur, M. 2008, The Astrophysical Journal, 685, 147, doi: 10.1086/590482
  * Newman (2008) Newman, J. A. 2008, The Astrophysical Journal, 684, 88, doi: 10.1086/589982
  * Newman & Gruen (2022) Newman, J. A., & Gruen, D. 2022, Annual Review of Astronomy and Astrophysics, 60, doi: 10.1146/annurev-astro-032122-014611
  * Ramachandra et al. (2021) Ramachandra, N., Chaves-Montero, J., Alarcon, A., et al. 2021, Machine learning synthetic spectra for probabilistic redshift estimation: SYTH-Z, arXiv, doi: 10.48550/ARXIV.2111.12118
  * Rau et al. (2019) Rau, M. M., Wilson, S., & Mandelbaum, R. 2019, Monthly Notices of the Royal Astronomical Society, 491, 4768, doi: 10.1093/mnras/stz3295
  * Sánchez & Bernstein (2018) Sánchez, C., & Bernstein, G. M. 2018, Monthly Notices of the Royal Astronomical Society, 483, 2801, doi: 10.1093/mnras/sty3222
  * Sánchez et al. (2020) Sánchez, C., Raveri, M., Alarcon, A., & Bernstein, G. M. 2020, Monthly Notices of the Royal Astronomical Society, 498, 2984, doi: 10.1093/mnras/staa2542
  * Schmidt et al. (2013) Schmidt, S. J., Mé nard, B., Scranton, R., Morrison, C., & McBride, C. K. 2013, Monthly Notices of the Royal Astronomical Society, 431, 3307, doi: 10.1093/mnras/stt410
  * Schneider et al. (2006) Schneider, M., Knox, L., Zhan, H., & Connolly, A. 2006, The Astrophysical Journal, 651, 14, doi: 10.1086/507675
  * Stanford et al. (2021) Stanford, S. A., Masters, D., Darvish, B., et al. 2021, The Astrophysical Journal Supplement Series, 256, 9, doi: 10.3847/1538-4365/ac0833
  * Tanaka (2015) Tanaka, M. 2015, The Astrophysical Journal, 801, 20, doi: 10.1088/0004-637x/801/1/20
  * The LSST Dark Energy Science Collaboration et al. (2018) The LSST Dark Energy Science Collaboration, Mandelbaum, R., Eifler, T., et al. 2018, The LSST Dark Energy Science Collaboration (DESC) Science Requirements Document, arXiv, doi: 10.48550/ARXIV.1809.01669
  * Weaver et al. (2022) Weaver, J. R., Kauffmann, O. B., Ilbert, O., et al. 2022, The Astrophysical Journal Supplement Series, 258, 11, doi: 10.3847/1538-4365/ac3078
  * Whitler et al. (2022) Whitler, L., Stark, D. P., Endsley, R., et al. 2022, Star formation histories of UV-luminous galaxies at $z\simeq 6.8$: implications for stellar mass assembly at early cosmic times, arXiv, doi: 10.48550/ARXIV.2206.05315