Technical Summary: Gaussian Process Latent Factor Regression for Low-Data, High-Dimensional Output Problems

Problem Statement

In scientific domains, regression tasks frequently require predicting high-dimensional structured outputs ($D_y$) from a limited number of training examples ($N$). While Multi-Output Gaussian Processes (MOGPs) excel in low-data regimes by modeling output correlations, they often struggle when $D_y$ is very large (e.g., $10^3$–$10^4$). Conversely, standard "compress-then-predict" pipelines, such as Principal Component Analysis followed by Gaussian Process regression (PCA-GP), handle high dimensionality by learning a basis for output reconstruction. However, these methods suffer from a fundamental mismatch: the learned basis is optimized to minimize reconstruction error (capturing all variance) rather than to maximize predictability from the input space. When outputs contain structured "nuisance" variation that is correlated across dimensions but independent of the input, PCA-GP may allocate modeling capacity to these unpredictable directions, degrading performance.

Methodology: Gaussian Process Latent Factor Regression (GPLFR)

The authors propose Gaussian Process Latent Factor Regression (GPLFR), a probabilistic model that couples representation learning and regression into a single objective.

Core Model Structure

GPLFR assumes that each high-dimensional output $y$ is a linear-Gaussian decoding of a low-dimensional latent state $z$ drawn from a Gaussian Process (GP) prior over the input space $x$.

1. Latent Prior: For each example $i$, a latent vector $z_i \in \mathbb{R}^{D_z}$ (where $D_z \ll D_y$) is generated. The latent functions $z^{(q)}(\cdot)$ for each dimension $q$ follow independent GP priors over the inputs $X$ with specific kernels $k_q$.
2. Linear Decoder: The output is modeled as $Y = ZW^\top + E$, where $W \in \mathbb{R}^{D_y \times D_z}$ represents the decoder weights and $E$ is observation noise.
3. Analytical Marginalization: Unlike standard latent variable models that require sampling or variational inference for the weights, GPLFR analytically marginalizes out the decoder weights $W$ (assuming a matrix-normal prior). This yields a collapsed likelihood that depends only on the latent variables $Z$ and the kernel hyperparameters.

Relationship to Existing Models

The paper frames GPLFR and the Linear Model of Coregionalization (LMC) as two different marginalizations of the same underlying joint linear-Gaussian factor model:

• LMC: Obtained by marginalizing the latent factors $Z$. This results in an explicit output covariance structure but treats the input-output mapping as a black box regarding the latent geometry.
• GPLFR: Obtained by marginalizing the decoder weights $W$. This results in a likelihood where the latent variables are constrained simultaneously by the GP prior over inputs and the reconstruction of outputs.

This distinction allows GPLFR to concentrate modeling capacity on learning a similarity structure over the input space ($x$) rather than enriching the output coupling structure, which is crucial when $N$ is small and $D_y$ is large.

Regularization and Optimization

To address the statistical unreliability of estimating a rich coregionalization matrix $B$ in low-data regimes, GPLFR typically sets $B = I$ (identity), relying on the latent geometry to capture correlations. To prevent the likelihood from overstating the information content of output dimensions when residuals are correlated (model misspecification), the authors introduce likelihood tempering with an inverse-temperature parameter $\beta \in (0, 1]$. Additionally, a small latent noise term $\lambda$ is added to the GP covariances to relax the requirement that latent scores be explained exactly by smooth GP functions, stabilizing the optimization landscape.

Key Contributions

1. Joint Representation Learning: The paper introduces a model that learns a low-dimensional latent representation jointly with the regression task, biasing the representation toward structures predictable from inputs rather than just high-variance structures.
2. Analytical Marginalization: By deriving a collapsed likelihood that marginalizes decoder weights, the method avoids the need for point-estimating a specific basis (as in PCA) or complex variational inference, scaling effectively to high-dimensional outputs.
3. Theoretical Connection: The work clarifies the mathematical relationship between GPLFR and the LMC, demonstrating how different marginalizations lead to different inductive biases suitable for different data regimes.
4. Scientific Application: The authors apply GPLFR to build the first spatially resolved emulator of Global Climate Models (GCMs) for rocky exoplanets, a task characterized by high-dimensional outputs ($D_y \approx 3 \times 10^4$) and very limited training data.

Experimental Results

The authors evaluate GPLFR against PCA-GP and other baselines on three tasks:

1. Synthetic Benchmark:

   • Setup: Data generated with a predictable signal component and a structured, input-independent "nuisance" component.
   • Findings: GPLFR significantly outperforms PCA-GP in the presence of structured nuisance variation. It achieves the same error rates as PCA-GP with approximately 4$\times$ fewer training data. GPLFR's latent representation captures more signal energy (77% vs. 54% for PCA-GP at the true rank) and less nuisance energy.
   • Limitation: In the noiseless limit where all variance is relevant to prediction, PCA-GP slightly outperforms standard GPLFR, though PCA-initialized GPLFR can close this gap.

2. Biomedical Optics (PyXOpto Emulation):

   • Setup: Predicting 500-dimensional radial reflectance profiles from tissue optical properties. The noise here is largely uncorrelated, making it a "PCA-friendly" regime.
   • Findings: Even when PCA assumptions hold, GPLFR demonstrates superior sample efficiency and representation efficiency, maintaining a lower error floor with fewer latent dimensions.

3. Exoplanet Climate Emulation:

   • Setup: Predicting 3D atmospheric states (temperature, humidity, winds, etc.) for rocky exoplanets from 8 continuous planet properties and a discrete GCM label. The dataset contains only 307 high-fidelity examples (plus 1,419 auxiliary examples).
   • Findings: GPLFR achieves the lowest Root Mean Squared Error (RMSE) and Energy Score across all variables compared to PPCA-ICM, PPCA-MLP, kNN, and training-mean baselines. Improvements over the best probabilistic baseline (PPCA-ICM) range from 8% to 34%. The results suggest that the structured, input-independent variation arising from differences in simulation setups (nuisance) confounds PCA-based methods, whereas GPLFR successfully filters this out.

Significance and Claims

The paper claims that GPLFR provides a principled, end-to-end alternative to compress-then-predict pipelines for low-data, high-dimensional regression problems. Its primary significance lies in its ability to prioritize predictable structure over high-variance structure, making it particularly effective when outputs contain structured noise or nuisance variation that is not correlated with the inputs.

The authors note that the practical cost of GPLFR is a more challenging optimization landscape due to the coupling between latent variables and kernel hyperparameters. However, they demonstrate that this can be mitigated through regularization (likelihood tempering and latent noise). The successful application to exoplanet climate modeling highlights the method's utility in scientific domains where data generation is expensive and outputs are complex, spatially resolved fields.