GPz: Gaussian Processes for Photometric Redshifts

• GPz is a supervised Bayesian framework that uses sparse, non-stationary Gaussian processes to provide both point estimates and calibrated predictive variances for photometric redshifts.
• It employs heteroscedastic uncertainty estimation and cost-sensitive learning to adjust for variable noise and non-uniform training data in large cosmological surveys.
• Its efficient sparse formulation scales as O(N m²) and decomposes uncertainty into epistemic and aleatoric components, guiding targeted improvements in observational campaigns.

GPz usually denotes Gaussian Processes for photometric redshift estimation, a supervised Bayesian machine-learning framework introduced as a non-stationary sparse Gaussian process with heteroscedastic uncertainty estimation for photometric redshifts (Almosallam et al., 2016). In the original usage, the “$z$” refers to redshift, and the method was developed for survey regimes in which spectroscopy is infeasible for the full galaxy sample, but accurate point estimates and well-characterized predictive variances remain essential for cosmological analyses (Hatfield et al., 2019). Across subsequent literature, GPz became both a specific photo-$z$ algorithm and a more general sparse heteroscedastic regression framework applied to uncertainty-aware surrogate modeling, domain-shift mitigation, and hybrid probabilistic inference (Hatfield et al., 2022).

1. Observational problem and methodological rationale

GPz emerged from the requirement that forthcoming cosmology surveys such as Euclid, the Large Synoptic Survey Telescope, and the Square Kilometre Array must rely predominantly on photometric rather than spectroscopic redshifts, because spectroscopy cannot be obtained for the full imaged galaxy population (Almosallam et al., 2016). In that setting, accurate redshift point estimates are insufficient by themselves: predictive variances are needed to optimize samples for weak lensing, baryon acoustic oscillations, and supernova analyses by trading completeness against reliability (Almosallam et al., 2016).

Two empirical properties of photometric-redshift data motivated the design of GPz. First, the noise is heteroscedastic: photometric measurement errors vary across filters, depths, and source populations, so the irreducible uncertainty is input-dependent rather than constant (Almosallam et al., 2016). Second, the spectroscopic training set is typically non-uniform in colour space and redshift, so uncertainty should also respond to local training density rather than only to measurement noise (Almosallam et al., 2016). Later studies preserved this framing and repeatedly treated GPz as especially useful when predictive uncertainty must separate lack-of-coverage effects from intrinsic noise, including in radio-selected AGN samples, DES-like large catalogs, and synthetic LSST-style benchmarks (Luken et al., 2023, Abdalla et al., 13 Aug 2025, Stylianou et al., 2022).

This dual emphasis distinguishes GPz from purely deterministic regressors and from methods whose uncertainty estimates are appended post hoc. In the literature, GPz is consistently described as yielding a point estimate plus predictive variance, with the variance intended to be scientifically actionable rather than merely diagnostic (Gomes et al., 2017).

2. Statistical construction and sparse non-stationary formulation

In its Gaussian-process context, the reference regression model is

$$f \sim \mathcal{GP}(0, k_\theta(x,x')), \quad y_n = f(x_n) + \epsilon_n, \quad \epsilon_n \sim \mathcal{N}(0, \sigma^2),$$

with homoscedastic marginal likelihood

$$p(\mathbf{y}|\mathbf{X},\sigma^2) = \mathcal{N}\!\big(0, \mathbf{K} + \sigma^2 \mathbf{I}\big).$$

GPz replaces the full $O(N^3)$ GP computation with a sparse basis-function construction that is mathematically equivalent to a sparse GP with inducing basis functions (Almosallam et al., 2016, Hatfield et al., 2019). In the form used repeatedly in later work, the predictive mean and variance proxy are parameterized as

$\mu(x) = \sum_{i=1}^m \phi_i(x) w_i,$

$$\sigma^2(x) = \exp\!\left(\sum_{i=1}^m \phi_i(x) v_i + b\right),$$

with localized squared-exponential basis functions

$$\phi_i(x) = \exp\!\left( -\frac{1}{2} (x - \mu_i)^T \Lambda_i^{-1} (x - \mu_i) \right).$$

Here, $\mu_i$ are basis centers, $\Lambda_i$ are basis-specific covariance or width parameters, $z$0 govern the mean function, $z$1 govern the log-variance function, and $z$2 is a variance bias term (Hatfield et al., 2019). This construction is called non-stationary because the effective local covariance structure varies across input space, and sparse because the number of basis functions $z$3 is much smaller than the number of training samples $z$4 (Almosallam et al., 2016).

The commonly used GPz variant in later applications is GPVC, in which each basis function has its own covariance. Other documented options include GPGC, GPVD, GPGD, GPVL, and GPGL, corresponding to different covariance-sharing structures (Rivera et al., 2018). Training typically uses L-BFGS or L-BFGS-B, often with feature standardization and joint learning of a linear mean component and heteroscedastic variance function (Hatfield et al., 2019, Lima et al., 2021).

The principal computational claim repeated across the literature is that GPz scales as

$z$5

rather than $z$6 for a full GP (Hatfield et al., 2019, Hatfield et al., 2020). This scaling enabled studies with $z$7, $z$8, $z$9, and $$f \sim \mathcal{GP}(0, k_\theta(x,x')), \quad y_n = f(x_n) + \epsilon_n, \quad \epsilon_n \sim \mathcal{N}(0, \sigma^2),$$0 basis functions in different data regimes, depending on dimensionality, sample size, and desired uncertainty fidelity (Rivera et al., 2018, Lima et al., 2021, Hatfield et al., 2022).

3. Predictive uncertainty, calibration, and cost-sensitive learning

A central feature of GPz is its explicit decomposition of predictive uncertainty into epistemic and aleatoric components. In the basis-function formalism, the total predictive variance is written as

$$f \sim \mathcal{GP}(0, k_\theta(x,x')), \quad y_n = f(x_n) + \epsilon_n, \quad \epsilon_n \sim \mathcal{N}(0, \sigma^2),$$1

or equivalently

$$f \sim \mathcal{GP}(0, k_\theta(x,x')), \quad y_n = f(x_n) + \epsilon_n, \quad \epsilon_n \sim \mathcal{N}(0, \sigma^2),$$2

The epistemic term $$f \sim \mathcal{GP}(0, k_\theta(x,x')), \quad y_n = f(x_n) + \epsilon_n, \quad \epsilon_n \sim \mathcal{N}(0, \sigma^2),$$3 reflects lack of nearby training data, while the aleatoric term $$f \sim \mathcal{GP}(0, k_\theta(x,x')), \quad y_n = f(x_n) + \epsilon_n, \quad \epsilon_n \sim \mathcal{N}(0, \sigma^2),$$4 or $$f \sim \mathcal{GP}(0, k_\theta(x,x')), \quad y_n = f(x_n) + \epsilon_n, \quad \epsilon_n \sim \mathcal{N}(0, \sigma^2),$$5 reflects input-dependent noise that does not vanish with denser sampling (Hatfield et al., 2019, Gomes et al., 2017). This distinction is repeatedly used in the literature to decide whether additional spectra, simulations, or experiments would actually reduce uncertainty.

The documented expression for the epistemic term in the ICF formulation is

$$f \sim \mathcal{GP}(0, k_\theta(x,x')), \quad y_n = f(x_n) + \epsilon_n, \quad \epsilon_n \sim \mathcal{N}(0, \sigma^2),$$6

where $$f \sim \mathcal{GP}(0, k_\theta(x,x')), \quad y_n = f(x_n) + \epsilon_n, \quad \epsilon_n \sim \mathcal{N}(0, \sigma^2),$$7 is the vector of basis responses and $$f \sim \mathcal{GP}(0, k_\theta(x,x')), \quad y_n = f(x_n) + \epsilon_n, \quad \epsilon_n \sim \mathcal{N}(0, \sigma^2),$$8 is the posterior covariance of the mean weights (Hatfield et al., 2019). In photometric-redshift work, the same conceptual split underpins sample selection, PDF diagnostics, and hybrid weighting schemes (Hatfield et al., 2022).

GPz also supports cost-sensitive learning. In the photo-$$f \sim \mathcal{GP}(0, k_\theta(x,x')), \quad y_n = f(x_n) + \epsilon_n, \quad \epsilon_n \sim \mathcal{N}(0, \sigma^2),$$9 literature, one reported weighting scheme is

$$p(\mathbf{y}|\mathbf{X},\sigma^2) = \mathcal{N}\!\big(0, \mathbf{K} + \sigma^2 \mathbf{I}\big).$$0

used to emphasize lower-redshift objects in one GAMA-based study (Gomes et al., 2017). In the S-PLUS implementation, the weights were normalized as $$p(\mathbf{y}|\mathbf{X},\sigma^2) = \mathcal{N}\!\big(0, \mathbf{K} + \sigma^2 \mathbf{I}\big).$$1 (Lima et al., 2021). In the covariate-shift setting of Gaussian-mixture augmentation, the weights were defined as

$$p(\mathbf{y}|\mathbf{X},\sigma^2) = \mathcal{N}\!\big(0, \mathbf{K} + \sigma^2 \mathbf{I}\big).$$2

with $$p(\mathbf{y}|\mathbf{X},\sigma^2) = \mathcal{N}\!\big(0, \mathbf{K} + \sigma^2 \mathbf{I}\big).$$3 and $$p(\mathbf{y}|\mathbf{X},\sigma^2) = \mathcal{N}\!\big(0, \mathbf{K} + \sigma^2 \mathbf{I}\big).$$4, thereby reweighting the training objective toward the target colour–magnitude distribution (Hatfield et al., 2020). Outside photo-$$p(\mathbf{y}|\mathbf{X},\sigma^2) = \mathcal{N}\!\big(0, \mathbf{K} + \sigma^2 \mathbf{I}\big).$$5, the inertial-confinement-fusion application used

$$p(\mathbf{y}|\mathbf{X},\sigma^2) = \mathcal{N}\!\big(0, \mathbf{K} + \sigma^2 \mathbf{I}\big).$$6

to prioritize high-yield regions of parameter space (Hatfield et al., 2019).

Calibration is commonly assessed through PIT, CRPS, QQ-based adjustments, and related metrics. A well-calibrated ensemble should yield a uniform PIT histogram, whereas concavity, convexity, or slope indicate overconfidence, underconfidence, or bias (Lima et al., 2021). In GPz-specific post-processing, a constant $$p(\mathbf{y}|\mathbf{X},\sigma^2) = \mathcal{N}\!\big(0, \mathbf{K} + \sigma^2 \mathbf{I}\big).$$7 shift per photo-$$p(\mathbf{y}|\mathbf{X},\sigma^2) = \mathcal{N}\!\big(0, \mathbf{K} + \sigma^2 \mathbf{I}\big).$$8 bin was optimized from QQ plots to reduce mean bias without changing the predictive variances (Gomes et al., 2017).

4. Performance in photometric redshift estimation

The original SDSS DR12 study reported that GPz substantially outperformed other machine-learning methods for photo-$$p(\mathbf{y}|\mathbf{X},\sigma^2) = \mathcal{N}\!\big(0, \mathbf{K} + \sigma^2 \mathbf{I}\big).$$9 estimation and associated variance prediction, including TPZ and ANNz2, while providing both Matlab and Python implementations (Almosallam et al., 2016). Subsequent work refined this picture rather than overturning it: GPz consistently remained competitive, but its strengths and weaknesses depended strongly on feature sets, representativeness, and whether the task emphasized point estimates, full $O(N^3)$0, or calibration.

In the S-PLUS DR1 analysis using 30 input features from S-PLUS, unWISE, colours, and morphology, GPz achieved for the full $O(N^3)$1 sample

$O(N^3)$2

compared with $O(N^3)$3 and $O(N^3)$4 for the dense-network and BNN+MDN models, and $O(N^3)$5 for BPZ2 (Lima et al., 2021). In that study, GPz PDFs were unimodal Gaussians and the PIT was described as “near ideal,” but with heavier tails near 0 and 1 than the deep-learning alternatives, consistent with the limitations of a unimodal approximation (Lima et al., 2021).

The GAMA/SDSS/UKIDSS study examined how GPz responds to richer features and calibration procedures. Adding near-infrared $O(N^3)$6 photometry and angular size improved the normalized RMSE from $O(N^3)$7 to $O(N^3)$8 under “Normal” cost-sensitive learning and from $O(N^3)$9 to $\mu(x) = \sum_{i=1}^m \phi_i(x) w_i,$0 under “Normalized” learning, with corresponding gains in $\mu(x) = \sum_{i=1}^m \phi_i(x) w_i,$1 and mean log-likelihood (Gomes et al., 2017). The abstract’s claim of roughly 15–20 per cent improvement refers to these typical gains across bins and metrics (Gomes et al., 2017). In the same work, QQ-based post-processing improved the bias by approximately 40 per cent, and replacing SDSS/UKIDSS photometry with higher-precision HSC photometry reduced RMSE by 17.4 per cent to 19.2 per cent, depending on the weighting regime (Gomes et al., 2017).

Under non-representative training, GPz remained strong for single-point estimates but became more vulnerable in population-level $\mu(x) = \sum_{i=1}^m \phi_i(x) w_i,$2 recovery. In simulations and SDSS/GAMA-based tests, GPz was slightly better than ANNz2 for single-point photo-$\mu(x) = \sum_{i=1}^m \phi_i(x) w_i,$3 in the representative case, but ANNz2 with CDF-based estimators was mildly better for reconstructing the full redshift distribution in deeper $\mu(x) = \sum_{i=1}^m \phi_i(x) w_i,$4-band cuts (Rivera et al., 2018). Training in colour rather than magnitude space was reported to recover reliable photo-$\mu(x) = \sum_{i=1}^m \phi_i(x) w_i,$5 estimates for about 42 per cent of galaxies relative to the deeper $\mu(x) = \sum_{i=1}^m \phi_i(x) w_i,$6-band cut, whereas the magnitude-representative subset constituted about 20 per cent (Rivera et al., 2018).

In radio-selected AGN-dominated samples, GPz had the lowest overall scatter among the tested regressors but not the lowest catastrophic-outlier fraction. The GMM-segmented GPz model reported

$\mu(x) = \sum_{i=1}^m \phi_i(x) w_i,$7

whereas $\mu(x) = \sum_{i=1}^m \phi_i(x) w_i,$8NN achieved the best $\mu(x) = \sum_{i=1}^m \phi_i(x) w_i,$9 (Luken et al., 2023). Yet GPz’s uncertainty estimates were operationally valuable: applying the certainty cut

$$\sigma^2(x) = \exp\!\left(\sum_{i=1}^m \phi_i(x) v_i + b\right),$$0

reduced the GPz catastrophic-outlier rate from $$\sigma^2(x) = \exp\!\left(\sum_{i=1}^m \phi_i(x) v_i + b\right),$$1 to $$\sigma^2(x) = \exp\!\left(\sum_{i=1}^m \phi_i(x) v_i + b\right),$$2 on the “certain” subset (Luken et al., 2023).

A large-scale DES DR2 application, trained on VIPERS and filtered with a K-d-tree reliability score $$\sigma^2(x) = \exp\!\left(\sum_{i=1}^m \phi_i(x) v_i + b\right),$$3, reported $$\sigma^2(x) = \exp\!\left(\sum_{i=1}^m \phi_i(x) v_i + b\right),$$4 and a catastrophic outlier rate of approximately 3 per cent on the reliable subset (Abdalla et al., 13 Aug 2025). That work emphasized that GPz’s heteroscedastic predictive variance can be used directly in constructing stacked $$\sigma^2(x) = \exp\!\left(\sum_{i=1}^m \phi_i(x) v_i + b\right),$$5 estimates and onion-like tomographic slices for cosmological mapping (Abdalla et al., 13 Aug 2025).

5. Augmentations, hybrid methods, and reuse beyond photo-$$\sigma^2(x) = \exp\!\left(\sum_{i=1}^m \phi_i(x) v_i + b\right),$$6

A major line of development combined GPz with explicit models of domain mismatch. The Gaussian-mixture augmentation framework modeled $$\sigma^2(x) = \exp\!\left(\sum_{i=1}^m \phi_i(x) v_i + b\right),$$7 and $$\sigma^2(x) = \exp\!\left(\sum_{i=1}^m \phi_i(x) v_i + b\right),$$8 separately in colour–magnitude space, used the density ratio for weighting, altered the validation split accordingly, and partitioned the input space into mixture regions with separate GPz models (Hatfield et al., 2020). In that setting, the combined “All” strategy was reported to reduce the bias by up to a half, improve PIT calibration at $$\sigma^2(x) = \exp\!\left(\sum_{i=1}^m \phi_i(x) v_i + b\right),$$9, increase the recovered number of high-$$\phi_i(x) = \exp\!\left( -\frac{1}{2} (x - \mu_i)^T \Lambda_i^{-1} (x - \mu_i) \right).$$0 galaxies in stacked posteriors by more than threefold relative to the unweighted baseline, and reduce runtime by approximately a factor of $$\phi_i(x) = \exp\!\left( -\frac{1}{2} (x - \mu_i)^T \Lambda_i^{-1} (x - \mu_i) \right).$$1 through mixture-specific division (Hatfield et al., 2020).

Another major extension embedded GPz in hybrid Bayesian consensus models with template fitting. In the COSMOS and XMM-LSS fields, GPz was run with $$\phi_i(x) = \exp\!\left( -\frac{1}{2} (x - \mu_i)^T \Lambda_i^{-1} (x - \mu_i) \right).$$2, GPVC covariances, normalization, joint mean learning, and GMM-All regional modeling, then combined hierarchically with galaxy-template and AGN-template PDFs (Hatfield et al., 2022). GPz alone achieved

$$\phi_i(x) = \exp\!\left( -\frac{1}{2} (x - \mu_i)^T \Lambda_i^{-1} (x - \mu_i) \right).$$3

while the hierarchical Bayesian consensus improved this to

$$\phi_i(x) = \exp\!\left( -\frac{1}{2} (x - \mu_i)^T \Lambda_i^{-1} (x - \mu_i) \right).$$4

against spectroscopy (Hatfield et al., 2022). The stage-two reliability weights were driven by the template $$\phi_i(x) = \exp\!\left( -\frac{1}{2} (x - \mu_i)^T \Lambda_i^{-1} (x - \mu_i) \right).$$5 and the GPz extrapolation fraction $$\phi_i(x) = \exp\!\left( -\frac{1}{2} (x - \mu_i)^T \Lambda_i^{-1} (x - \mu_i) \right).$$6, explicitly coupling template credibility to GPz’s own uncertainty decomposition (Hatfield et al., 2022).

GPz has also been repurposed outside redshift estimation. In inertial-confinement fusion, it was used as a fast surrogate model for implosion yields, with $$\phi_i(x) = \exp\!\left( -\frac{1}{2} (x - \mu_i)^T \Lambda_i^{-1} (x - \mu_i) \right).$$7, $$\phi_i(x) = \exp\!\left( -\frac{1}{2} (x - \mu_i)^T \Lambda_i^{-1} (x - \mu_i) \right).$$8, $$\phi_i(x) = \exp\!\left( -\frac{1}{2} (x - \mu_i)^T \Lambda_i^{-1} (x - \mu_i) \right).$$9, and heteroscedastic joint learning (Hatfield et al., 2019). That study reported training times of about 5 s and prediction times of about 0.02 s on a laptop for the main 5D problem, and used the ratio $\mu_i$0 to identify regions where additional simulations would reduce epistemic uncertainty (Hatfield et al., 2019). The same paper treated GPz’s original redshift meaning as historical but demonstrated that the underlying sparse heteroscedastic regression machinery was not specific to photo-$\mu_i$1 (Hatfield et al., 2019).

6. Limitations, diagnostic metrics, and terminological ambiguity

The most persistent limitation in the GPz literature is the assumption that each object’s predictive redshift PDF is a single Gaussian mode. This can be adequate when the conditional mapping is locally unimodal, but it misses genuinely multi-modal colour–redshift degeneracies. The S-PLUS comparison attributed part of the advantage of BNN+MDN models to their ability to represent multi-modal $\mu_i$2, whereas GPz’s unimodal Gaussian PDFs produced heavier PIT tails near 0 and 1 (Lima et al., 2021).

A second recurring issue is sensitivity to training-set imperfection. In synthetic LSST-like experiments, line-confusion fractions above roughly 1 per cent caused a substantial drop in photo-$\mu_i$3 quality, and inverse redshift incompleteness below a pivot redshift of 1.5 also degraded performance (Stylianou et al., 2022). Under strong covariate shift in the StratLearn-z benchmark, GPz reached

$\mu_i$4

while the propensity-stratified conditional-density approach was markedly less affected (Moretti et al., 2024). This suggests that GPz’s heteroscedastic variance is useful but does not by itself resolve severe mismatch between $\mu_i$5 and $\mu_i$6.

Recent diagnostic studies have emphasized that some conventional metrics can be misleading for GPz under broad or miscentered Gaussian predictions. In RAIL-based stress tests with realistic survey degradations, Kullback–Leibler divergence, Wasserstein distance, and PIT were singled out as particularly informative metrics for assessing GPz under training-set imperfection, whereas iterative outlier fractions could understate failure when the model returned extremely broad Gaussians (Crafford et al., 15 Jan 2026). The same work argued that inverse redshift incompleteness alone lacks the complexity needed to represent realistic spectroscopic selection effects for GPz benchmarking (Crafford et al., 15 Jan 2026).

Model-capacity and approximation issues also remain. Too few basis functions can under-represent uncertainty at domain edges, finite $\mu_i$7 can underestimate uncertainty far from sampled regions, and poor initialization of basis locations or anisotropic covariances can affect convergence (Hatfield et al., 2019). Later hybrid and GMM-based studies can be read as partial responses to these failure modes, because they regionalize the mapping and alter training emphasis rather than relying on one global GPz fit (Hatfield et al., 2020, Hatfield et al., 2022).

Finally, the acronym itself is no longer unique. Although GPz is most strongly associated with Gaussian Processes for photometric redshifts (Almosallam et al., 2016), unrelated later papers use GPZ for Ghost Probe Zones in autonomous driving (Qu et al., 23 Apr 2025) and for GPU-Accelerated Lossy Compressor for Particle Data (Li et al., 14 Aug 2025). A plausible implication is that literature searches for “GPZ” increasingly require domain-specific disambiguation, whereas “GPz” with a lowercase final letter remains the clearest identifier for the photometric-redshift method.