Personalized Federated Hierarchical Gaussian Processes for Privacy-Preserving Modeling of Heterogeneous Distributed Systems

# Personalized Federated Hierarchical Gaussian Processes for Privacy-Preserving Modeling of Heterogeneous Distributed Systems
Source: [https://arxiv.org/html/2609.19337](https://arxiv.org/html/2609.19337)
Hao Yan††thanks:Arizona State University\. Email:haoyan@asu\.edu\.

###### Abstract

We present Personalized Federated Hierarchical Gaussian Processes \(pFedHGP\) for probabilistic regression and classification when data are distributed across heterogeneous clients\. Each client’s latent function decomposes into \(i\) a shared global component, \(ii\) a client\-specific deviation that shares the global kernel structure, and \(iii\) a flexible local residual\. Sparse inducing\-variable approximations and federated variational inference keep raw data local while the server synchronizes only low\-dimensional statistics for the shared component\. Full predictive distributions support uncertainty\-aware decisions\. In application studies, pFedHGP attains perfect fault classification in press tonnage monitoring using 13\.77% of labeled cycles and recovers geographic zones in federated air\-quality modeling without centralizing station\-level time series\. An Instantaneous Linear Mixing Model viewpoint links the hierarchy to multi\-output Gaussian processes for correlated sensors\.

Keywords:Federated Learning, Gaussian Process, Hierarchical Model, Uncertainty Quantification, Distributed Sensing, Personalized Modeling

Funding:This work is funded by DOE NEUP DE\-NE0009383\.

## 1Introduction

Privacy\-preserving learning from spatially and institutionally distributed measurements has become central to modern industrial and urban analytics\([Kontar et al\., 2021](https://arxiv.org/html/2609.19337#bib.bib16)\)\. With the rapid proliferation of Internet of Things \(IoT\) devices—distributed environmental sensors, manufacturing equipment monitors, and smart city infrastructure—many deployments naturally adopt federated architectures where data remain on edge devices\. In this setting, models must capture population\-level regularities while respecting data locality, bandwidth limits, and heterogeneous operating conditions\. Representative application domains include smart healthcare\([Xu et al\., 2021](https://arxiv.org/html/2609.19337#bib.bib32)\), predictive maintenance in industrial manufacturing\([Pruckovskaja et al\., 2023](https://arxiv.org/html/2609.19337#bib.bib25)\), and connected vehicle networks for smart cities\([Chellapandi et al\., 2023](https://arxiv.org/html/2609.19337#bib.bib4)\)\.

A fundamental challenge in such distributed systems is handling heterogeneity: physical systems often exhibit site\-specific behaviors due to differing operational conditions, environmental factors, equipment variations, or localized disturbances\. When learning shared models from non\-i\.i\.d\. federated clients, this heterogeneity can significantly degrade prediction quality and yield poorly calibrated uncertainty\([Ye et al\., 2023](https://arxiv.org/html/2609.19337#bib.bib36)\)\. The core tension is to retain a population\-level representation of common dynamics while allowing client\-specific systematic shifts and flexible local residuals\.

To address heterogeneity, a stream of research in personalized learning, including personalized federated learning\([Deng et al\., 2020](https://arxiv.org/html/2609.19337#bib.bib8)\), personalized PCA\([Shi and Al Kontar, 2024](https://arxiv.org/html/2609.19337#bib.bib27)\), and personalized Tucker decomposition\([Hu et al\., 2025](https://arxiv.org/html/2609.19337#bib.bib12)\), has focused on decomposing client data into shared global components and local personalized components\. These approaches typically extract a joint low\-rank subspace to represent the global structure while allowing client\-specific variations in complementary subspaces\. Essentially, they perform latent subspace decomposition by identifying shared bases along with client\-specific coefficients, thereby achieving a balance between collective representation and individualized adaptation\.

Inspired by this decomposition paradigm, our objective is to develop a federated probabilistic framework that separates shared structure from client\-specific variation in settings with correlated multi\-sensor outputs\. Such scenarios arise in manufacturing systems with multiple strain channels, urban environmental monitoring networks, and other distributed sensing stacks\. We leverage Gaussian processes \(GPs\), which naturally encode dependence through kernels, provide principled uncertainty quantification \(UQ\), and remain flexible in low\-data regimes\([Seeger, 2004](https://arxiv.org/html/2609.19337#bib.bib26)\)\. Rigorous uncertainty estimates are particularly important for monitoring and planning under safety constraints\. GPs have strong records in surrogate modeling, spatiotemporal forecasting, robotics, and sensor networks\([Deisenroth et al\., 2013](https://arxiv.org/html/2609.19337#bib.bib7);[Wilson and Adams, 2013](https://arxiv.org/html/2609.19337#bib.bib31)\)\.

Despite existing efforts, three major challenges remain in developing GP\-based models for federated heterogeneous environments:

- •Lack of explicit hierarchical calibration structure:Many federated GP methods rely on a single shared kernel with local fine\-tuning, without a transparent decomposition into shared global structure, structured client deviations, and flexible local residuals\.
- •Inadequate support for correlated multi\-sensor clients:Existing federated GP methods often treat outputs as independent or rely on global pooling, limiting explicit modeling of cross\-sensor dependence within each client\.
- •Insufficient uncertainty\-aware adaptation:Many federated methods do not deliver calibrated uncertainty that balances global consistency with client\-specific adaptation, which matters for risk\-aware decisions under heterogeneity\.

To address these gaps, we propose personalized Federated Hierarchical Gaussian Process \(pFedHGP\), trained with privacy\-preserving federated variational inference\. pFedHGP embeds a three\-level hierarchical decomposition: \(1\) a global GP capturing population\-level structure, \(2\) client\-specific deviation GPs that share the global kernel geometry, and \(3\) client\-specific local residual GPs with independent kernels and inducing geometries\. An Instantaneous Linear Mixing Model \(ILMM\) viewpoint\([Bruinsma et al\., 2020](https://arxiv.org/html/2609.19337#bib.bib2)\)clarifies how correlated multi\-sensor outputs arise from latent mixing; see Section[3\.9](https://arxiv.org/html/2609.19337#S3.SS9)\.

We summarize our contributions below:

- •Hierarchical decomposition for heterogeneous federated clients:We decompose each client’s latent function into shared global structure, structured deviations on a shared basis, and flexible local residuals, enabling interpretable monitoring under non\-i\.i\.d\. clients\.
- •Privacy\-preserving federated variational inference:We develop a two\-stage variational procedure in which the server synchronizes only compact global statistics while clients update local variational factors on private data\.
- •Uncertainty\-aware multi\-sensor extension:Full Gaussian predictive laws support marginal and joint uncertainty across outputs; the ILMM connection highlights how vector\-valued sensing maps to the hierarchy\.

The federated formulation is needed for three reasons\. First, raw sensor trajectories are naturally owned by local clients and may contain proprietary or governance\-restricted information\. Second, high\-resolution time series create communication and storage costs when naively centralized\. Third, the clients are statistically heterogeneous, so a purely pooled model can obscure local systematic shifts\. pFedHGP addresses these issues by keeping raw observations local, transmitting only model\-level updates for the shared representation, and preserving client\-specific behavior through the deviation and local components\. We use the term privacy\-preserving in this data\-locality sense; the present work does not provide formal differential privacy or secure aggregation guarantees\.

The remainder of this paper is organized as follows\. Section[2](https://arxiv.org/html/2609.19337#S2)reviews related works\. Section[3](https://arxiv.org/html/2609.19337#S3)introduces pFedHGP\. Section[4](https://arxiv.org/html/2609.19337#S4)presents a multi\-output synthetic study with ablations, robustness checks under assumption mismatch, and uncertainty metrics\. Section[5](https://arxiv.org/html/2609.19337#S5)presents realistic case studies\. Section[6](https://arxiv.org/html/2609.19337#S6)discusses limitations\. Section[7](https://arxiv.org/html/2609.19337#S7)concludes\.

## 2Background and Literature Review

The evolution of Gaussian Processes \(GPs\) from centralized models to distributed and federated architectures has been driven by the dual imperatives of computational scalability and data privacy\. We categorize the existing literature into five major domains: GP\-based calibration and model\-discrepancy methods \(Section[2\.1](https://arxiv.org/html/2609.19337#S2.SS1)\), which assume centralized data access; distributed GP methods \(Section[2\.2](https://arxiv.org/html/2609.19337#S2.SS2)\), which prioritize scalability over privacy; federated GP approaches \(Section[2\.3](https://arxiv.org/html/2609.19337#S2.SS3)\), which preserve privacy but struggle with client heterogeneity; personalized federated learning \(Section[2\.4](https://arxiv.org/html/2609.19337#S2.SS4)\), which combines global and local adaptations; and advanced federated architectures \(Section[2\.5](https://arxiv.org/html/2609.19337#S2.SS5)\), which handle multi\-output and multi\-fidelity data at the cost of increased complexity\.

### 2\.1Gaussian Processes for Calibration and Model Discrepancy

Gaussian Processes provide a probabilistic framework for model calibration, uncertainty quantification, and correction of discrepancies between simulation models and physical observations\.[Chakraborty et al\. \(2020\)](https://arxiv.org/html/2609.19337#bib.bib3)highlight the role of GP surrogates in bridging the gap between high\-fidelity simulations and noisy sensor data\. In structural dynamics,[Kessels et al\. \(2022\)](https://arxiv.org/html/2609.19337#bib.bib15)use GPs as inverse mapping models to update physical parameters and quantify calibration uncertainty\. These centralized approaches leverage the flexibility of GPs to model the systematic bias, or model\-discrepancy, function\([Kennedy and O’Hagan, 2001](https://arxiv.org/html/2609.19337#bib.bib14)\)that captures the difference between an idealized physical model and observed data\. However, traditional calibration methods typically assume centralized data access, which becomes problematic in distributed industrial settings where data privacy is paramount\.

### 2\.2Distributed Gaussian Processes

Early distributed GP methods were primarily designed to address the cubic computational complexity of standard GPs\([Kontoudis and Stilwell, 2022](https://arxiv.org/html/2609.19337#bib.bib17);[Liu et al\., 2018](https://arxiv.org/html/2609.19337#bib.bib19)\)\. These approaches typically partition data across multiple nodes to parallelize computation\.[Tresp \(2000\)](https://arxiv.org/html/2609.19337#bib.bib30)introduced the Bayesian Committee Machine \(BCM\), which aggregates predictions from independent GP experts by weighing them with their inverse covariance\.[Deisenroth and Ng \(2015\)](https://arxiv.org/html/2609.19337#bib.bib6)extended this to a distributed variational framework \(DgGP\), allowing for the robust aggregation of expert predictions\. Similarly,[Ng and Deisenroth \(2014\)](https://arxiv.org/html/2609.19337#bib.bib22)and[Nguyen and Bonilla \(2014\)](https://arxiv.org/html/2609.19337#bib.bib23)proposed mixture\-of\-experts \(MoE\) models that assign data to local experts based on hierarchical or information\-theoretic criteria\.

In signal processing,[Xu et al\. \(2019\)](https://arxiv.org/html/2609.19337#bib.bib33)and[Yin and Gunnarsson \(2017\)](https://arxiv.org/html/2609.19337#bib.bib37)developed distributed recursive algorithms for wireless traffic and sensor networks, focusing on efficient hyperparameter synchronization\.[Peng et al\. \(2017\)](https://arxiv.org/html/2609.19337#bib.bib24)and[Gal et al\. \(2014\)](https://arxiv.org/html/2609.19337#bib.bib9)further improved scalability by leveraging asynchronous updates and stochastic variational inference, enabling GPs to handle massive datasets\.

The primary limitation of distributed GPs is their assumption that data is either i\.i\.d\. or artificially partitioned for computational convenience\. They often rely on shared global kernels or inducing points, lacking mechanisms to model the intrinsic heterogeneity \(non\-i\.i\.d\. distributions\) found in natural federated settings\. Furthermore, these methods prioritize scalability over privacy, often requiring the exchange of gradients or statistics that may leak sensitive information\([Yin et al\., 2020](https://arxiv.org/html/2609.19337#bib.bib38)\)\.

### 2\.3Federated Gaussian Processes

With the rise of privacy concerns, the focus shifted to Federated Learning \(FL\), where data remains local\. Federated GPs aim to learn a global GP model without sharing raw data\.[Yu et al\. \(2022\)](https://arxiv.org/html/2609.19337#bib.bib39)proposed Federated Bayesian Neural Regression \(FedBNR\), utilizing random features to approximate a scalable global kernel\.[Yurochkin et al\. \(2019\)](https://arxiv.org/html/2609.19337#bib.bib41)introduced Bayesian Nonparametric Federated Learning \(BNFed\) to match neural network weights across clients probabilistically\.

Recent works have explored decentralized approaches\.[Kontoudis and Stilwell \(2022\)](https://arxiv.org/html/2609.19337#bib.bib17)developed fully decentralized GPs for multi\-agent systems, while[Llorente et al\. \(2025\)](https://arxiv.org/html/2609.19337#bib.bib20)introduced robust, online decentralized GPs adaptable to dynamic environments\.[Thorgeirsson and Gauterin \(2020\)](https://arxiv.org/html/2609.19337#bib.bib28)provided a probabilistic extension to the standard FedAvg algorithm, aggregating moments of client distributions to quantify global uncertainty\.

While these methods improve privacy, they generally enforce a single global model across all clients\. This ”one\-size\-fits\-all” approach struggles when clients exhibit significant statistical heterogeneity\. A single global kernel cannot capture the diverse, client\-specific patterns inherent in many real\-world applications, leading to suboptimal performance for individual participants\.

### 2\.4Personalized Federated Learning

To address heterogeneity, Personalized Federated Learning \(PFL\) integrates global knowledge with client\-specific adaptations\. While Gaussian processes offer one avenue for personalization, broader PFL frameworks have explored various decomposition strategies\.[Shi and Al Kontar \(2024\)](https://arxiv.org/html/2609.19337#bib.bib27)introduced Personalized PCA \(PerPCA\), which explicitly decouples client data into shared global features and unique local features, enabling robust dimensionality reduction in heterogeneous settings\. Similarly,[Hu et al\. \(2025\)](https://arxiv.org/html/2609.19337#bib.bib12)proposed Personalized Tucker Decomposition \(perTucker\) for tensor data, modeling commonalities via a shared core tensor while capturing peculiarities through client\-specific factor matrices\.

In the context of GPs, the seminal work on pFedGP by[Achituve et al\. \(2021\)](https://arxiv.org/html/2609.19337#bib.bib1)combines Deep Kernel Learning \(DKL\) with local GPs\. In this framework, a neural network feature extractor \(the deep kernel\) is learned globally, while each client maintains a local GP on the extracted features\. This allows for shared representation learning while preserving local flexibility\. Subsequent research has refined this paradigm; for instance,[Yang et al\. \(2024\)](https://arxiv.org/html/2609.19337#bib.bib34)incorporated reinforcement learning to optimize personalization strategies in heterogeneous environments\.

Non\-GP personalized methods like PerPCA and perTucker are powerful for feature extraction but lack the built\-in uncertainty quantification needed for safety\-critical monitoring and decision\-making tasks\. On the other hand, personalized federated GPs like pFedGP often suffer from high computational burdens on client devices \(O⁡\(N3\)O\(N^\{3\}\)complexity\) and communication overheads associated with deep kernel parameters\.

### 2\.5Advanced Federated Architectures

Recent advancements have extended federated GPs to handle more complex data structures, such as multi\-output and multi\-fidelity data\.[Chung and Al Kontar \(2024\)](https://arxiv.org/html/2609.19337#bib.bib5)introduced Federated Multi\-output Gaussian Processes \(FedMGP\), which learn shared latent functions to model correlations across multiple outputs \(e\.g\., spatial or temporal tasks\)\.[Gao and Chung \(2024\)](https://arxiv.org/html/2609.19337#bib.bib10)further enhanced this by developing automated mechanisms for selecting latent variables in federated settings\.

In reliability engineering,[Jeong et al\. \(2025\)](https://arxiv.org/html/2609.19337#bib.bib13)proposed Fed\-Joint, a framework coupling nonlinear degradation signals with failure events using federated multi\-output GPs\. Similarly,[Yue et al\. \(2024\)](https://arxiv.org/html/2609.19337#bib.bib40)explored federated data analytics for linear models, providing a foundation for scalable multi\-fidelity modeling\.

These advanced architectures introduce significant structural complexity\. They often require sophisticated variational inference schemes\([Bruinsma et al\., 2020](https://arxiv.org/html/2609.19337#bib.bib2)\)to remain tractable, and the need to capture correlations across multiple outputs or fidelities increases the communication overhead and the difficulty of optimizing the global variational objective\.

## 3Methodology

We present the Personalized Federated Hierarchical Gaussian Process \(pFedHGP\) framework for privacy\-preserving federated regression and calibration under client heterogeneity\. Section[3\.1](https://arxiv.org/html/2609.19337#S3.SS1)reviews sparse Gaussian processes and the Variational Free Energy framework that enables scalable inference\. Section[3\.2](https://arxiv.org/html/2609.19337#S3.SS2)introduces the hierarchical generative model decomposing observations into global, client\-specific deviation, and local residual components\. Section[3\.3](https://arxiv.org/html/2609.19337#S3.SS3)establishes the low\-rank inducing variable representation that provides computational efficiency\. Section[3\.4](https://arxiv.org/html/2609.19337#S3.SS4)presents the variational sparse GP framework with mean\-field assumptions and additive covariance decomposition\. Section[3\.5](https://arxiv.org/html/2609.19337#S3.SS5)derives the federated variational inference algorithm with explicit ELBO decomposition\. Section[3\.6](https://arxiv.org/html/2609.19337#S3.SS6)details the pFedHGP training algorithm alternating between local optimization and global aggregation\. Section[3\.7](https://arxiv.org/html/2609.19337#S3.SS7)covers implementation details including Cholesky factorization, whitened parameterization, and hyperparameter optimization\. Section[3\.8](https://arxiv.org/html/2609.19337#S3.SS8)describes practical usage for classification and clustering tasks\. Finally, Section[3\.9](https://arxiv.org/html/2609.19337#S3.SS9)connects pFedHGP to the Instantaneous Linear Mixing Model framework\.

### 3\.1Preliminaries: Sparse Gaussian Processes

Standard Gaussian processes \(GPs\) suffer from cubic computational complexityO⁡\(N3\)O\(N^\{3\}\), making them unsuitable for large\-scale or federated learning\. Sparse Gaussian processes \(SGPs\) address this by introducing a set ofM≪NM\\ll Ninducing variables𝐮\\mathbf\{u\}at locations𝐙=\[𝐳1,…,𝐳M\]\\mathbf\{Z\}=\[\\mathbf\{z\}\_\{1\},\\dots,\\mathbf\{z\}\_\{M\}\]\. For a single GP with priorf∼𝒢​𝒫​\(0,k\)f\\sim\\mathcal\{GP\}\(0,k\), the joint prior over function values𝐟=f⁡\(𝐗\)\\mathbf\{f\}=f\(\\mathbf\{X\}\)and inducing values𝐮=f⁡\(𝐙\)\\mathbf\{u\}=f\(\\mathbf\{Z\}\)is Gaussian, and the conditional distribution

p⁡\(𝐟∣𝐮\)=𝒩⁡\(𝐊n​m​𝐊m​m−1​𝐮,𝐊n​n−𝐊n​m​𝐊m​m−1​𝐊m​n\)p\(\\mathbf\{f\}\\mid\\mathbf\{u\}\)=\\mathcal\{N\}\\Big\(\\mathbf\{K\}\_\{nm\}\\mathbf\{K\}\_\{mm\}^\{\-1\}\\mathbf\{u\},\\;\\mathbf\{K\}\_\{nn\}\-\\mathbf\{K\}\_\{nm\}\\mathbf\{K\}\_\{mm\}^\{\-1\}\\mathbf\{K\}\_\{mn\}\\Big\)\(1\)is*exact*\. Here𝐊n​n=k⁡\(𝐗,𝐗\)\\mathbf\{K\}\_\{nn\}=k\(\\mathbf\{X\},\\mathbf\{X\}\),𝐊m​m=k⁡\(𝐙,𝐙\)\\mathbf\{K\}\_\{mm\}=k\(\\mathbf\{Z\},\\mathbf\{Z\}\), and𝐊n​m=k⁡\(𝐗,𝐙\)\\mathbf\{K\}\_\{nm\}=k\(\\mathbf\{X\},\\mathbf\{Z\}\)\.

Titsias’s variational formulation views the inducing variables as variational parameters\. Introducing a variational distribution of the form

q\(𝐟,𝐮\)=p\(𝐟∣𝐮\)q\(𝐮\),q\(𝐮\)=𝒩\(𝐦,𝐒\),\\displaystyle q\(\\mathbf\{f\},\\mathbf\{u\}\)=p\(\\mathbf\{f\}\\mid\\mathbf\{u\}\)\\,q\(\\mathbf\{u\}\),\\quad q\(\\mathbf\{u\}\)=\\mathcal\{N\}\(\\mathbf\{m\},\\mathbf\{S\}\),yields the evidence lower bound \(ELBO\) under the Variational Free Energy \(VFE\) framework\.

###### Proposition 1\(Sparse GP ELBO\)

Under the VFE framework with inducing variables𝐮\\mathbf\{u\}at locations𝐙\\mathbf\{Z\}, the ELBO takes the form:

ℒVFE=𝔼q⁡\(𝐮\)\[logp\(𝐲∣𝐮\)\]−KL\(q\(𝐮\)∥p\(𝐮\)\)−12​σ2Tr\(𝐊n​n−𝐐n​n\),\\mathcal\{L\}\_\{\\mathrm\{VFE\}\}=\\mathbb\{E\}\_\{q\(\\mathbf\{u\}\)\}\\big\[\\log p\(\\mathbf\{y\}\\mid\\mathbf\{u\}\)\\big\]\-\\mathrm\{KL\}\\big\(q\(\\mathbf\{u\}\)\\\|p\(\\mathbf\{u\}\)\\big\)\-\\frac\{1\}\{2\\sigma^\{2\}\}\\operatorname\{Tr\}\\big\(\\mathbf\{K\}\_\{nn\}\-\\mathbf\{Q\}\_\{nn\}\\big\),\(2\)wherep⁡\(𝐲∣𝐮\)=𝒩⁡\(𝐊n​m​𝐊m​m−1​𝐮,σ2​𝐈\)p\(\\mathbf\{y\}\\mid\\mathbf\{u\}\)=\\mathcal\{N\}\(\\mathbf\{K\}\_\{nm\}\\mathbf\{K\}\_\{mm\}^\{\-1\}\\mathbf\{u\},\\sigma^\{2\}\\mathbf\{I\}\),𝐐n​n=𝐊n​m​𝐊m​m−1​𝐊m​n\\mathbf\{Q\}\_\{nn\}=\\mathbf\{K\}\_\{nm\}\\mathbf\{K\}\_\{mm\}^\{\-1\}\\mathbf\{K\}\_\{mn\}, andσ2\\sigma^\{2\}is the observation noise variance\. The first two terms correspond to a low\-rank approximating model with covariance𝐐n​n\+σ2​𝐈\\mathbf\{Q\}\_\{nn\}\+\\sigma^\{2\}\\mathbf\{I\}, and the trace term penalizes the discrepancy between the full prior covariance𝐊n​n\\mathbf\{K\}\_\{nn\}and its low\-rank approximation\.

The VFE method is originally proposed in\([Titsias, 2009](https://arxiv.org/html/2609.19337#bib.bib29)\), and a proof of Proposition[1](https://arxiv.org/html/2609.19337#Thmproposition1)is provided in Appendix[B](https://arxiv.org/html/2609.19337#A2)\.

The key viewpoint we adopt is that SGPs are not defined by modifying the GP prior, but by introducing inducing variables and optimizing a variational bound of the form \([2](https://arxiv.org/html/2609.19337#S3.E2)\)\. Our pFedHGP method extends this variational sparse GP view to a federated setting by decomposing the latent function space into shared and local components, each equipped with its own inducing variables\.

### 3\.2Motivation and Formulation

For clarity, Table[1](https://arxiv.org/html/2609.19337#S3.T1)summarizes the main notation used throughout the methodology\.

Table 1:Summary of main notation used in pFedHGP\.We consider a federated learning scenario involvingTTclients, each of which holds private data inaccessible to a central server\. Let clientiipossess a local dataset\{\(𝐱i,j,yi,j\)\}j=1ni\\\{\(\\mathbf\{x\}\_\{i,j\},y\_\{i,j\}\)\\\}\_\{j=1\}^\{n\_\{i\}\}, where𝐱i,j∈ℝd\\mathbf\{x\}\_\{i,j\}\\in\\mathbb\{R\}^\{d\}is the input andyi,j∈ℝy\_\{i,j\}\\in\\mathbb\{R\}is the corresponding scalar response\. We assume that each observation is composed of three latent components—a shared global function, a client\-specific deviation, and a local residual—plus additive Gaussian noise\.

###### Assumption 1\(pFedHGP Generative Model\)

The observationyi,jy\_\{i,j\}for clientiiis decomposed as

yi,j=fg​\(𝐱i,j\)\+fδ,i​\(𝐱i,j\)\+fi​\(𝐱i,j\)\+εi,j,y\_\{i,j\}=f\_\{g\}\(\\mathbf\{x\}\_\{i,j\}\)\+f\_\{\\delta,i\}\(\\mathbf\{x\}\_\{i,j\}\)\+f\_\{i\}\(\\mathbf\{x\}\_\{i,j\}\)\+\\varepsilon\_\{i,j\},\(3\)wherefgf\_\{g\}is the shared global function,fδ,if\_\{\\delta,i\}is the client\-specific deviation from the global trend, andfif\_\{i\}is the local site\-specific variation\. The noise variables are i\.i\.d\.εi,j∼𝒩⁡\(0,σ2\)\\varepsilon\_\{i,j\}\\sim\\mathcal\{N\}\(0,\\sigma^\{2\}\)\. We endow these components with independent GP priors:

fg∼𝒢𝒫\(0,kg\),fδ,i∼𝒢𝒫\(0,kδ\),fi∼𝒢𝒫\(0,ki\),i=1,…,T,f\_\{g\}\\sim\\mathcal\{GP\}\\big\(0,k\_\{g\}\\big\),\\quad f\_\{\\delta,i\}\\sim\\mathcal\{GP\}\\big\(0,k\_\{\\delta\}\\big\),\\quad f\_\{i\}\\sim\\mathcal\{GP\}\\big\(0,k\_\{i\}\\big\),\\quad i=1,\\dots,T,\(4\)withfgf\_\{g\},\{fδ,i\}\\\{f\_\{\\delta,i\}\\\}, and\{fi\}\\\{f\_\{i\}\\\}mutually independent\. In our implementation we takekδ=ϕ​kgk\_\{\\delta\}=\\phi\\,k\_\{g\}so that the deviation processes share the same structure as the global process but are scaled by a factorϕ\>0\\phi\>0\.

The three latent components have distinct statistical roles\. The global processfgf\_\{g\}represents the population\-level structure shared across clients\. The deviation processfδ,if\_\{\\delta,i\}represents a client\-specific systematic departure from this shared structure, but it is constrained to use the same kernel family and inducing basis asfgf\_\{g\}throughkδ=ϕ​kgk\_\{\\delta\}=\\phi k\_\{g\}\. The parameterϕ\\phiis a globally shared hyperparameter\. It serves as a shared prior scale for the expected magnitude of systematic deviations across clients\. Althoughϕ\\phiis shared, the realized deviation functionfδ,if\_\{\\delta,i\}is client\-specific because each client has its own variational posterior\. It therefore captures smooth calibration shifts or persistent site\-specific offsets that remain aligned with the shared functional structure\. By contrast, the local processfif\_\{i\}uses a client\-specific kernelkik\_\{i\}and inducing locations𝐙i\\mathbf\{Z\}\_\{i\}, allowing it to absorb idiosyncratic residual behavior, local disturbances, or shorter\-scale variation that should not be forced into the shared global basis\. Thus,fδ,if\_\{\\delta,i\}handles structured client\-level calibration relative to the global component, whereasfif\_\{i\}provides flexible local residual adaptation\.

##### Rationale for hierarchical priors\.

The scaled deviation kernelkδ=ϕ​kgk\_\{\\delta\}=\\phi k\_\{g\}follows the model\-discrepancy paradigm of[Kennedy and O’Hagan \(2001\)](https://arxiv.org/html/2609.19337#bib.bib14): systematic client shifts are expressed on the same functional basis as the shared component, with population\-level intensityϕ\\philearned globally \(Section[3\.6](https://arxiv.org/html/2609.19337#S3.SS6)\)\. Sharing inducing locations𝐙g\\mathbf\{Z\}\_\{g\}betweenfgf\_\{g\}andfδ,if\_\{\\delta,i\}keeps both layers in a common subspace so the asymmetric KL penalties in the ELBO apply to the same basis coefficients\. Independent local kernelskik\_\{i\}and inducing sets𝐙i\\mathbf\{Z\}\_\{i\}absorb idiosyncratic residuals that need not follow the global geometry\. Section[4\.3](https://arxiv.org/html/2609.19337#S4.SS3)\(Scenario B\) stress\-tests deviation\-kernel mismatch when the fitted model retains an RBFkgk\_\{g\}but the simulator uses a periodic deviation kernel\.

The shared kernel formkδ=ϕ​kgk\_\{\\delta\}=\\phi k\_\{g\}is a regularizing design choice rather than an assertion that all clients deviate in exactly the same way\. It constrains systematic deviations to the same functional basis as the global component, which makes the deviation interpretable as a calibration offset around the shared structure\. Using separate inducing locations or a completely independent deviation kernel would increase flexibility, but it would also weaken the competition between the global and deviation layers and increase communication and optimization complexity\. We therefore use the shared basis for the structured deviation layer and reserve the independent kernelkik\_\{i\}and inducing set𝐙i\\mathbf\{Z\}\_\{i\}for idiosyncratic local residuals\. Section[4\.3](https://arxiv.org/html/2609.19337#S4.SS3)evaluates this choice under mismatch scenarios in which the true deviation kernel, additivity, or component independence is violated\.

##### Identifiability and regularized attribution\.

Becausekδ=ϕ​kgk\_\{\\delta\}=\\phi k\_\{g\}, the marginal prior for a single client’s sumfg\+fδ,if\_\{g\}\+f\_\{\\delta,i\}has covariance\(1\+ϕ\)​kg\(1\+\\phi\)k\_\{g\}\. Hence, if a client were analyzed in isolation and only this summed process were observed, the split betweenfgf\_\{g\}andfδ,if\_\{\\delta,i\}would not be uniquely identifiable\. This gap is about attribution, not about whether the client’s observable function can be predicted\. Squared\-exponential kernels are universal on compact input domains, so an RBF hierarchy can still approximatefg\+fδ,if\_\{g\}\+f\_\{\\delta,i\}as a single continuous function ofxx\. What a single client cannot do is uniquely name the two layers inside that sum\. The separation in pFedHGP instead relies on the multi\-client federated hierarchy:fgf\_\{g\}is one shared process used by all clients, whereasfδ,if\_\{\\delta,i\}is independent across clients\. Therefore, fori≠ji\\neq j,

Cov⁡\{yi​\(x\),yj​\(x′\)\}=kg​\(x,x′\),\\operatorname\{Cov\}\\\{y\_\{i\}\(x\),y\_\{j\}\(x^\{\\prime\}\)\\\}=k\_\{g\}\(x,x^\{\\prime\}\),while the deviation process contributes only to within\-client covariance\. Under this generative specification, off\-diagonal cross\-client covariance therefore depends only on the shared global processfgf\_\{g\}\(not on a claim thatfgf\_\{g\}is uniquely identified from data\)\. The ELBO further regularizes this attribution through an asymmetric complexity cost: a common pattern represented byfgf\_\{g\}incurs one global KL penalty, whereas representing the same pattern redundantly through the client\-specific deviations incurs∑iKL\(q\(𝜹i\)∥p\(𝜹i\)\)\\sum\_\{i\}\\mathrm\{KL\}\\bigl\(q\(\\bm\{\\delta\}\_\{i\}\)\\\|p\(\\bm\{\\delta\}\_\{i\}\)\\bigr\)\. This encourages functional patterns common across clients to be represented byfgf\_\{g\}and nonredundant client\-specific departures to be represented byfδ,if\_\{\\delta,i\}\. We therefore interpret the decomposition as a regularized hierarchical attribution supported by multi\-client data, rather than as an unconstrained unique decomposition from a single\-client marginal process\.

For clientiiwith inputs𝐗i=\[𝐱i,1,…,𝐱i,ni\]⊤\\mathbf\{X\}\_\{i\}=\[\\mathbf\{x\}\_\{i,1\},\\dots,\\mathbf\{x\}\_\{i,n\_\{i\}\}\]^\{\\top\}, we define the stacked latent vectors𝐟g,i=fg​\(𝐗i\)\\mathbf\{f\}\_\{g,i\}=f\_\{g\}\(\\mathbf\{X\}\_\{i\}\),𝐟δ,i=fδ,i​\(𝐗i\)\\mathbf\{f\}\_\{\\delta,i\}=f\_\{\\delta,i\}\(\\mathbf\{X\}\_\{i\}\), and𝐟i=fi​\(𝐗i\)\\mathbf\{f\}\_\{i\}=f\_\{i\}\(\\mathbf\{X\}\_\{i\}\)\. The total latent effect is denoted by𝐟itot=𝐟g,i\+𝐟δ,i\+𝐟i\\mathbf\{f\}\_\{i\}^\{\\text\{tot\}\}=\\mathbf\{f\}\_\{g,i\}\+\\mathbf\{f\}\_\{\\delta,i\}\+\\mathbf\{f\}\_\{i\}, with observation vector𝐲i\|𝐟itot∼𝒩⁡\(𝐟itot,σ2​𝐈ni\)\\mathbf\{y\}\_\{i\}\\mid\\mathbf\{f\}\_\{i\}^\{\\text\{tot\}\}\\sim\\mathcal\{N\}\(\\mathbf\{f\}\_\{i\}^\{\\text\{tot\}\},\\sigma^\{2\}\\mathbf\{I\}\_\{n\_\{i\}\}\)\.

### 3\.3Low\-Rank Inducing Variable Representation

To obtain a scalable inference scheme, we augment this hierarchical model with inducing variables for each GP component, yielding a low\-rank basis decomposition analogous to personalized dimensionality reduction methods\.

###### Proposition 2\(Low\-Rank Decomposition for Mean Representation\)

Let𝐙g∈ℝM×d\\mathbf\{Z\}\_\{g\}\\in\\mathbb\{R\}^\{M\\times d\}denote shared global inducing locations and𝐙i∈ℝMi×d\\mathbf\{Z\}\_\{i\}\\in\\mathbb\{R\}^\{M\_\{i\}\\times d\}denote client\-specific inducing locations\. Define the inducing variables

𝐮g=fg​\(𝐙g\)∈ℝM,𝜹i=fδ,i​\(𝐙g\)∈ℝM,𝐮i=fi​\(𝐙i\)∈ℝMi\.\\mathbf\{u\}\_\{g\}=f\_\{g\}\(\\mathbf\{Z\}\_\{g\}\)\\in\\mathbb\{R\}^\{M\},\\quad\\bm\{\\delta\}\_\{i\}=f\_\{\\delta,i\}\(\\mathbf\{Z\}\_\{g\}\)\\in\\mathbb\{R\}^\{M\},\\quad\\mathbf\{u\}\_\{i\}=f\_\{i\}\(\\mathbf\{Z\}\_\{i\}\)\\in\\mathbb\{R\}^\{M\_\{i\}\}\.\(5\)Under the GP priors specified in Section[3\.2](https://arxiv.org/html/2609.19337#S3.SS2), these inducing variables are jointly Gaussian:

𝐮g∼𝒩⁡\(𝟎,𝐊g​g\),𝜹i∼𝒩⁡\(𝟎,ϕ​𝐊g​g\),𝐮i∼𝒩⁡\(𝟎,𝐊i​i\),\\mathbf\{u\}\_\{g\}\\sim\\mathcal\{N\}\(\\mathbf\{0\},\\mathbf\{K\}\_\{gg\}\),\\quad\\bm\{\\delta\}\_\{i\}\\sim\\mathcal\{N\}\(\\mathbf\{0\},\\phi\\mathbf\{K\}\_\{gg\}\),\\quad\\mathbf\{u\}\_\{i\}\\sim\\mathcal\{N\}\(\\mathbf\{0\},\\mathbf\{K\}\_\{ii\}\),\(6\)where𝐊g​g=kg​\(𝐙g,𝐙g\)\\mathbf\{K\}\_\{gg\}=k\_\{g\}\(\\mathbf\{Z\}\_\{g\},\\mathbf\{Z\}\_\{g\}\)and𝐊i​i=ki​\(𝐙i,𝐙i\)\\mathbf\{K\}\_\{ii\}=k\_\{i\}\(\\mathbf\{Z\}\_\{i\},\\mathbf\{Z\}\_\{i\}\)\.

Conditioned on the inducing variables, the latent function values at clientii’s observations𝐗i\\mathbf\{X\}\_\{i\}have the following structure\. Let𝐟g,i=fg​\(𝐗i\)\\mathbf\{f\}\_\{g,i\}=f\_\{g\}\(\\mathbf\{X\}\_\{i\}\),𝐟δ,i=fδ,i​\(𝐗i\)\\mathbf\{f\}\_\{\\delta,i\}=f\_\{\\delta,i\}\(\\mathbf\{X\}\_\{i\}\), and𝐟i=fi​\(𝐗i\)\\mathbf\{f\}\_\{i\}=f\_\{i\}\(\\mathbf\{X\}\_\{i\}\)denote the latent function vectors evaluated at𝐗i\\mathbf\{X\}\_\{i\}\. Define the interpolation matrices

𝐀g,i=𝐊g​\(𝐗i,𝐙g\)​𝐊g​g−1,𝐀δ,i=𝐊g​\(𝐗i,𝐙g\)​𝐊g​g−1,𝐀i=𝐊i​\(𝐗i,𝐙i\)​𝐊i​i−1\.\\mathbf\{A\}\_\{g,i\}=\\mathbf\{K\}\_\{g\}\(\\mathbf\{X\}\_\{i\},\\mathbf\{Z\}\_\{g\}\)\\mathbf\{K\}\_\{gg\}^\{\-1\},\\quad\\mathbf\{A\}\_\{\\delta,i\}=\\mathbf\{K\}\_\{g\}\(\\mathbf\{X\}\_\{i\},\\mathbf\{Z\}\_\{g\}\)\\mathbf\{K\}\_\{gg\}^\{\-1\},\\quad\\mathbf\{A\}\_\{i\}=\\mathbf\{K\}\_\{i\}\(\\mathbf\{X\}\_\{i\},\\mathbf\{Z\}\_\{i\}\)\\mathbf\{K\}\_\{ii\}^\{\-1\}\.\(7\)The conditional mean of the total latent function admits the low\-rank basis expansion

𝔼\[𝐟itot∣𝐮g,𝜹i,𝐮i\]=𝐀g,i⏟shared basis\(𝐮g\+𝜹i\)\+𝐀i⏟local basis𝐮i,\\mathbb\{E\}\\big\[\\mathbf\{f\}\_\{i\}^\{\\mathrm\{tot\}\}\\mid\\mathbf\{u\}\_\{g\},\\bm\{\\delta\}\_\{i\},\\mathbf\{u\}\_\{i\}\\big\]=\\underbrace\{\\mathbf\{A\}\_\{g,i\}\}\_\{\\text\{shared basis\}\}\(\\mathbf\{u\}\_\{g\}\+\\bm\{\\delta\}\_\{i\}\)\+\\underbrace\{\\mathbf\{A\}\_\{i\}\}\_\{\\text\{local basis\}\}\\mathbf\{u\}\_\{i\},\(8\)where𝐟itot=𝐟g,i\+𝐟δ,i\+𝐟i\\mathbf\{f\}\_\{i\}^\{\\mathrm\{tot\}\}=\\mathbf\{f\}\_\{g,i\}\+\\mathbf\{f\}\_\{\\delta,i\}\+\\mathbf\{f\}\_\{i\}is the total latent function\. The proof follows from standard GP conditional mean formulas; see Appendix[D](https://arxiv.org/html/2609.19337#A4)for details\.

##### Connection to personalized PCA\.

The structure in \([8](https://arxiv.org/html/2609.19337#S3.E8)\) parallels personalized PCA \(perPCA\)\([Shi and Al Kontar, 2024](https://arxiv.org/html/2609.19337#bib.bib27)\), which decomposes client data into shared and unique features using linear bases\. In perPCA, each client’s data matrix is represented as𝐗i=𝐔𝐕g⊤\+𝐔i​𝐕i⊤\+𝐄i\\mathbf\{X\}\_\{i\}=\\mathbf\{U\}\\mathbf\{V\}\_\{g\}^\{\\top\}\+\\mathbf\{U\}\_\{i\}\\mathbf\{V\}\_\{i\}^\{\\top\}\+\\mathbf\{E\}\_\{i\}, where𝐔𝐕g⊤\\mathbf\{U\}\\mathbf\{V\}\_\{g\}^\{\\top\}is a shared low\-rank component and𝐔i​𝐕i⊤\\mathbf\{U\}\_\{i\}\\mathbf\{V\}\_\{i\}^\{\\top\}is a client\-specific low\-rank component\. Similarly, pFedHGP decomposes the latent function space into:

- •Ashared low\-rank subspacespanned by the global kernel basis𝐤g​\(𝐱,𝐙g\)⊤​𝐊g​g−1\\mathbf\{k\}\_\{g\}\(\\mathbf\{x\},\\mathbf\{Z\}\_\{g\}\)^\{\\top\}\\mathbf\{K\}\_\{gg\}^\{\-1\}with two components: the fixed global weights𝐮g\\mathbf\{u\}\_\{g\}\(common to all clients\) and client\-specific deviation weights𝜹i\\bm\{\\delta\}\_\{i\}\(capturing systematic offsets\)\.
- •Aclient\-specific low\-rank subspacespanned by the local kernel basis𝐤i​\(𝐱,𝐙i\)⊤​𝐊i​i−1\\mathbf\{k\}\_\{i\}\(\\mathbf\{x\},\\mathbf\{Z\}\_\{i\}\)^\{\\top\}\\mathbf\{K\}\_\{ii\}^\{\-1\}with local weights𝐮i\\mathbf\{u\}\_\{i\}\.

However, pFedHGP extends beyond perPCA in three critical ways: \(1\) it uses kernel\-induced nonlinear bases rather than linear bases, enabling flexible nonparametric function approximation; \(2\) it provides full probabilistic modeling with uncertainty quantification for all components via GP priors, including structured residual covariance, whereas perPCA yields point estimates with only simple \(often i\.i\.d\.\) noise; and \(3\) it naturally handles functional data and regression tasks through the GP framework, while perPCA is designed for matrix factorization\. The explicit regularization𝜹i∼𝒩⁡\(𝟎,ϕ​𝐊g​g\)\\bm\{\\delta\}\_\{i\}\\sim\\mathcal\{N\}\(\\mathbf\{0\},\\phi\\mathbf\{K\}\_\{gg\}\)ensures that client\-specific adjustments remain structurally consistent with the shared physical model encoded in the global kernel\.

### 3\.4Variational Sparse GP Framework

To enable scalable and tractable inference in the hierarchical model, we adopt the variational sparse Gaussian process framework based on the Variational Free Energy \(VFE\) method\([Titsias, 2009](https://arxiv.org/html/2609.19337#bib.bib29)\)\. This framework avoids restrictive assumptions about the observation covariance structure \(unlike FITC\) while maintaining computational efficiency\.

###### Assumption 2\(Mean Field Variational Inference\)

We assume a factorized variational posterior distribution over the inducing variables:

q⁡\(𝐮g,\{𝜹i\},\{𝐮i\}\)=q⁡\(𝐮g\)​∏i=1Tq⁡\(𝜹i\)​q​\(𝐮i\),q\(\\mathbf\{u\}\_\{g\},\\\{\\bm\{\\delta\}\_\{i\}\\\},\\\{\\mathbf\{u\}\_\{i\}\\\}\)=q\(\\mathbf\{u\}\_\{g\}\)\\prod\_\{i=1\}^\{T\}q\(\\bm\{\\delta\}\_\{i\}\)q\(\\mathbf\{u\}\_\{i\}\),\(9\)where each factor is a multivariate Gaussian distribution\.

This mean\-field assumption decouples the inference across the hierarchical components while preserving the full covariance structure of the inducing variables within each component\. Crucially, we do*not*assume that the conditional posterior over observations is diagonal \(the FITC assumption\)\. Instead, we leverage the property that the variational lower bound \(ELBO\) only depends on the marginal variances of the residual process\.

###### Proposition 3\(Low\-Rank Covariance Decomposition\)

Due to the independence of the hierarchical components \(fg⟂fδ,i⟂fif\_\{g\}\\perp f\_\{\\delta,i\}\\perp f\_\{i\}\) in the prior, the marginal variance of the total latent function at any input𝐱\\mathbf\{x\}is the sum of the component variances\. Within the VFE framework\([Titsias, 2009](https://arxiv.org/html/2609.19337#bib.bib29)\), the expected log\-likelihood term in the ELBO relies only on these marginal variances, allowing the residual termTr⁡\(𝐑i\)\\operatorname\{Tr\}\(\\mathbf\{R\}\_\{i\}\)to be computed additively as:

Tr⁡\(𝐑i\)=Tr⁡\(𝐊g,i​i−𝐐g,i​i\)\+Tr⁡\(𝐊δ,i​i−𝐐δ,i​i\)\+Tr⁡\(𝐊i,i​i−𝐐i,i​i\),\\operatorname\{Tr\}\(\\mathbf\{R\}\_\{i\}\)=\\operatorname\{Tr\}\(\\mathbf\{K\}\_\{g,ii\}\-\\mathbf\{Q\}\_\{g,ii\}\)\+\\operatorname\{Tr\}\(\\mathbf\{K\}\_\{\\delta,ii\}\-\\mathbf\{Q\}\_\{\\delta,ii\}\)\+\\operatorname\{Tr\}\(\\mathbf\{K\}\_\{i,ii\}\-\\mathbf\{Q\}\_\{i,ii\}\),\(10\)where𝐑i\\mathbf\{R\}\_\{i\}represents the dense residual covariance matrix of the total latent function\. The low\-rank approximation matrices𝐐⋅,i​i\\mathbf\{Q\}\_\{\\cdot,ii\}are defined as in Proposition[2](https://arxiv.org/html/2609.19337#Thmproposition2)\.

See Appendix[E](https://arxiv.org/html/2609.19337#A5)for the derivation showing that the expected log\-likelihood term in the ELBO involves the trace of the residual covariance𝔼q​\[log⁡p⁡\(𝐲\|𝐟\)\]∝−12​σ2​Tr⁡\(𝐊n​n−𝐐n​n\)\\mathbb\{E\}\_\{q\}\[\\log p\(\\mathbf\{y\}\|\\mathbf\{f\}\)\]\\propto\-\\frac\{1\}\{2\\sigma^\{2\}\}\\operatorname\{Tr\}\(\\mathbf\{K\}\_\{nn\}\-\\mathbf\{Q\}\_\{nn\}\), which decomposes additively due to the linearity of the trace operator and the independence of components\.

𝐊g\\mathbf\{K\}\_\{g\}Global Physics\+\+𝐊δ\\mathbf\{K\}\_\{\\delta\}Calibration\+\+𝐊i\\mathbf\{K\}\_\{i\}Local Noise==𝐊t​o​t​a​l\\mathbf\{K\}\_\{total\}Observation CovarianceLow\-Rank𝐐g\\mathbf\{Q\}\_\{g\}Low\-Rank𝐐i\\mathbf\{Q\}\_\{i\}Dense Matrix\(Approximated via VFE\)Figure 1:Structure of the additive covariance in pFedHGP\. While the components are independent \(block\-diagonal in latent space\), they sum to form the observation covariance\. The ELBO in the VFE framework allows us to efficiently compute this by utilizing only the diagonal elements of the residual𝐊−𝐐\\mathbf\{K\}\-\\mathbf\{Q\}\.This result is significant because it justifies the efficiency of pFedHGP without making the restrictive FITC assumption that off\-diagonal data correlations are zero\. We simply exploit the property that maximizing the ELBO only requires matching the marginal variances of the residual process\. Figure[1](https://arxiv.org/html/2609.19337#S3.F1)illustrates this additive structure\.

### 3\.5Federated Variational Inference

Given the hierarchical GP generative model in \([3](https://arxiv.org/html/2609.19337#S3.E3)\), the joint distribution over all latent variables and observations can be written as

p⁡\(\{𝐲i\},fg,\{fδ,i\},\{fi\},𝐮g,\{𝜹i\},\{𝐮i\}\)\\displaystyle p\\big\(\\\{\\mathbf\{y\}\_\{i\}\\\},f\_\{g\},\\\{f\_\{\\delta,i\}\\\},\\\{f\_\{i\}\\\},\\mathbf\{u\}\_\{g\},\\\{\\bm\{\\delta\}\_\{i\}\\\},\\\{\\mathbf\{u\}\_\{i\}\\\}\\big\)=p⁡\(𝐮g\)​∏i=1Tp⁡\(𝜹i\)​p​\(𝐮i\)\\displaystyle=p\(\\mathbf\{u\}\_\{g\}\)\\prod\_\{i=1\}^\{T\}p\(\\bm\{\\delta\}\_\{i\}\)\\,p\(\\mathbf\{u\}\_\{i\}\)\\;\(11\)×p⁡\(fg∣𝐮g\)​∏i=1Tp⁡\(fδ,i∣𝜹i\)​p​\(fi∣𝐮i\)\\displaystyle\\times p\(f\_\{g\}\\mid\\mathbf\{u\}\_\{g\}\)\\prod\_\{i=1\}^\{T\}p\(f\_\{\\delta,i\}\\mid\\bm\{\\delta\}\_\{i\}\)\\,p\(f\_\{i\}\\mid\\mathbf\{u\}\_\{i\}\)\\;×∏i=1Tp⁡\(𝐲i∣𝐟itot\),\\displaystyle\\times\\prod\_\{i=1\}^\{T\}p\(\\mathbf\{y\}\_\{i\}\\mid\\mathbf\{f\}\_\{i\}^\{\\mathrm\{tot\}\}\),where𝐟itot=fg​\(𝐗i\)\+fδ,i​\(𝐗i\)\+fi​\(𝐗i\)\\mathbf\{f\}\_\{i\}^\{\\mathrm\{tot\}\}=f\_\{g\}\(\\mathbf\{X\}\_\{i\}\)\+f\_\{\\delta,i\}\(\\mathbf\{X\}\_\{i\}\)\+f\_\{i\}\(\\mathbf\{X\}\_\{i\}\)andp⁡\(𝐲i∣𝐟itot\)=𝒩⁡\(𝐟itot,σ2​𝐈\)p\(\\mathbf\{y\}\_\{i\}\\mid\\mathbf\{f\}\_\{i\}^\{\\mathrm\{tot\}\}\)=\\mathcal\{N\}\(\\mathbf\{f\}\_\{i\}^\{\\mathrm\{tot\}\},\\sigma^\{2\}\\mathbf\{I\}\)\.

Exact posterior inference is analytically tractable in principle, but computationally prohibitive in the federated, large\-scale setting due to theO⁡\(N3\)O\(N^\{3\}\)coupling across all clients and components\. We therefore follow the variational sparse GP strategy and introduce a variational distribution of the form

q⁡\(fg,\{fδ,i\},\{fi\},𝐮g,\{𝜹i\},\{𝐮i\}\)=p⁡\(fg∣𝐮g\)​∏i=1Tp⁡\(fδ,i∣𝜹i\)​p​\(fi∣𝐮i\)​q​\(𝐮g\)​∏i=1Tq⁡\(𝜹i\)​q​\(𝐮i\),q\(f\_\{g\},\\\{f\_\{\\delta,i\}\\\},\\\{f\_\{i\}\\\},\\mathbf\{u\}\_\{g\},\\\{\\bm\{\\delta\}\_\{i\}\\\},\\\{\\mathbf\{u\}\_\{i\}\\\}\)=p\(f\_\{g\}\\mid\\mathbf\{u\}\_\{g\}\)\\prod\_\{i=1\}^\{T\}p\(f\_\{\\delta,i\}\\mid\\bm\{\\delta\}\_\{i\}\)\\,p\(f\_\{i\}\\mid\\mathbf\{u\}\_\{i\}\)\\;q\(\\mathbf\{u\}\_\{g\}\)\\prod\_\{i=1\}^\{T\}q\(\\bm\{\\delta\}\_\{i\}\)\\,q\(\\mathbf\{u\}\_\{i\}\),\(12\)where

q⁡\(𝐮g\)=𝒩⁡\(𝐦g,𝐒g\),q⁡\(𝜹i\)=𝒩⁡\(𝐦δi,𝐒δi\),q⁡\(𝐮i\)=𝒩⁡\(𝐦ui,𝐒ui\)q\(\\mathbf\{u\}\_\{g\}\)=\\mathcal\{N\}\(\\mathbf\{m\}\_\{g\},\\mathbf\{S\}\_\{g\}\),\\quad q\(\\bm\{\\delta\}\_\{i\}\)=\\mathcal\{N\}\(\\mathbf\{m\}\_\{\\delta\_\{i\}\},\\mathbf\{S\}\_\{\\delta\_\{i\}\}\),\\quad q\(\\mathbf\{u\}\_\{i\}\)=\\mathcal\{N\}\(\\mathbf\{m\}\_\{u\_\{i\}\},\\mathbf\{S\}\_\{u\_\{i\}\}\)are free Gaussian factors\. The Gaussian conditionalsp⁡\(fg∣𝐮g\)p\(f\_\{g\}\\mid\\mathbf\{u\}\_\{g\}\),p⁡\(fδ,i∣𝜹i\)p\(f\_\{\\delta,i\}\\mid\\bm\{\\delta\}\_\{i\}\), andp⁡\(fi∣𝐮i\)p\(f\_\{i\}\\mid\\mathbf\{u\}\_\{i\}\)are those implied by the GP priors and inducing locations, and the latent function values are analytically integrated out\.

##### Hierarchical ELBO\.

Following Proposition[1](https://arxiv.org/html/2609.19337#Thmproposition1), we integrate out the latent function values𝐟itot\\mathbf\{f\}\_\{i\}^\{\\mathrm\{tot\}\}at𝐗i\\mathbf\{X\}\_\{i\}conditional on\(𝐮g,𝜹i,𝐮i\)\(\\mathbf\{u\}\_\{g\},\\bm\{\\delta\}\_\{i\},\\mathbf\{u\}\_\{i\}\)and adopt the VFE variational distributionq⁡\(𝐟itot∣𝐮g,𝜹i,𝐮i\)​q​\(𝐮g\)​q​\(𝜹i\)​q​\(𝐮i\)q\(\\mathbf\{f\}\_\{i\}^\{\\mathrm\{tot\}\}\\mid\\mathbf\{u\}\_\{g\},\\bm\{\\delta\}\_\{i\},\\mathbf\{u\}\_\{i\}\)\\,q\(\\mathbf\{u\}\_\{g\}\)q\(\\bm\{\\delta\}\_\{i\}\)q\(\\mathbf\{u\}\_\{i\}\)with mean\-fieldqqfactors\. The inducing\-conditional model for clientiiis then

p⁡\(𝐲i∣𝐮g,𝜹i,𝐮i\):=𝒩⁡\(𝐀g,i​\(𝐮g\+𝜹i\)\+𝐀i​𝐮i,σ2​𝐈ni\),p\(\\mathbf\{y\}\_\{i\}\\mid\\mathbf\{u\}\_\{g\},\\bm\{\\delta\}\_\{i\},\\mathbf\{u\}\_\{i\}\):=\\mathcal\{N\}\\big\(\\mathbf\{A\}\_\{g,i\}\(\\mathbf\{u\}\_\{g\}\+\\bm\{\\delta\}\_\{i\}\)\+\\mathbf\{A\}\_\{i\}\\mathbf\{u\}\_\{i\},\\,\\sigma^\{2\}\\mathbf\{I\}\_\{n\_\{i\}\}\\big\),\(13\)
which is the hierarchical analogue ofp⁡\(𝐲∣𝐮\)p\(\\mathbf\{y\}\\mid\\mathbf\{u\}\)in \([2](https://arxiv.org/html/2609.19337#S3.E2)\) \(a low\-rank Gaussian in the inducing variables only\)\. It does*not*yet include the sparse GP residual covariance: under VFE, marginalizing𝐟itot\\mathbf\{f\}\_\{i\}^\{\\mathrm\{tot\}\}contributes the additive penalty−12​σ2​Tr⁡\(𝐑i\)\-\\frac\{1\}\{2\\sigma^\{2\}\}\\operatorname\{Tr\}\(\\mathbf\{R\}\_\{i\}\)withTr⁡\(𝐑i\)\\operatorname\{Tr\}\(\\mathbf\{R\}\_\{i\}\)given by \([10](https://arxiv.org/html/2609.19337#S3.E10)\)\. Moreover, becausep⁡\(𝐲i∣𝐮g,𝜹i,𝐮i\)p\(\\mathbf\{y\}\_\{i\}\\mid\\mathbf\{u\}\_\{g\},\\bm\{\\delta\}\_\{i\},\\mathbf\{u\}\_\{i\}\)is linear in\(𝐮g,𝜹i,𝐮i\)\(\\mathbf\{u\}\_\{g\},\\bm\{\\delta\}\_\{i\},\\mathbf\{u\}\_\{i\}\), the expectation𝔼q​\[log⁡p⁡\(𝐲i∣𝐮g,𝜹i,𝐮i\)\]\\mathbb\{E\}\_\{q\}\[\\log p\(\\mathbf\{y\}\_\{i\}\\mid\\mathbf\{u\}\_\{g\},\\bm\{\\delta\}\_\{i\},\\mathbf\{u\}\_\{i\}\)\]also propagates posterior covariance in each inducing block\. Proposition[5](https://arxiv.org/html/2609.19337#Thmproposition5)states the resulting closed form \(squared\-error, inducing\-covariance traces, andTr⁡\(𝐑i\)\\operatorname\{Tr\}\(\\mathbf\{R\}\_\{i\}\)\); Theorem[1](https://arxiv.org/html/2609.19337#Thmtheorem1)combines this with the KL decomposition below\.

The evidence lower bound \(ELBO\) for pFedHGP, expressed in terms of the global, client\-specific deviation, and local inducing variables, is

ℒ\\displaystyle\\mathcal\{L\}:=∑i=1T𝔼q⁡\(𝐮g\)​q​\(𝜹i\)​q​\(𝐮i\)​\[log⁡p⁡\(𝐲i∣𝐮g,𝜹i,𝐮i\)\]\\displaystyle:=\\sum\_\{i=1\}^\{T\}\\mathbb\{E\}\_\{q\(\\mathbf\{u\}\_\{g\}\)\\,q\(\\bm\{\\delta\}\_\{i\}\)\\,q\(\\mathbf\{u\}\_\{i\}\)\}\\bigl\[\\log p\(\\mathbf\{y\}\_\{i\}\\mid\\mathbf\{u\}\_\{g\},\\bm\{\\delta\}\_\{i\},\\mathbf\{u\}\_\{i\}\)\\bigr\]\(14\)−KL\(q\(𝐮g,\{𝜹i\},\{𝐮i\}\)∥p\(𝐮g\)∏i=1Tp\(𝜹i\)p\(𝐮i\)\),\\displaystyle\-\\mathrm\{KL\}\\bigl\(q\(\\mathbf\{u\}\_\{g\},\\\{\\bm\{\\delta\}\_\{i\}\\\},\\\{\\mathbf\{u\}\_\{i\}\\\}\)\\,\\\|\\,p\(\\mathbf\{u\}\_\{g\}\)\\prod\_\{i=1\}^\{T\}p\(\\bm\{\\delta\}\_\{i\}\)p\(\\mathbf\{u\}\_\{i\}\)\\bigr\),
where the expectation𝔼q​\[log⁡p⁡\(𝐲i∣𝐮g,𝜹i,𝐮i\)\]\\mathbb\{E\}\_\{q\}\[\\log p\(\\mathbf\{y\}\_\{i\}\\mid\\mathbf\{u\}\_\{g\},\\bm\{\\delta\}\_\{i\},\\mathbf\{u\}\_\{i\}\)\]is understood as the full VFE expected log\-likelihood in \([16](https://arxiv.org/html/2609.19337#S3.E16)\) \(not the plug\-in log\-density of \([13](https://arxiv.org/html/2609.19337#S3.E13)\) evaluated only at posterior means\)\. Proposition[4](https://arxiv.org/html/2609.19337#Thmproposition4)shows that the overall KL term decomposes into component\-wise KL divergences, enabling efficient federated optimization\.

###### Proposition 4\(Decomposability of KL Terms\)

Under the mean\-field variational posterior \(Assumption[2](https://arxiv.org/html/2609.19337#Thmassumption2)\), the KL divergence between the variational posterior and the prior factorizes across components:

KL\(q\(𝐮g,\{𝜹i\},\{𝐮i\}\)∥p\(𝐮g\)∏i=1Tp\(𝜹i\)p\(𝐮i\)\)=\\displaystyle\\mathrm\{KL\}\\bigl\(q\(\\mathbf\{u\}\_\{g\},\\\{\\bm\{\\delta\}\_\{i\}\\\},\\\{\\mathbf\{u\}\_\{i\}\\\}\)\\,\\\|\\,p\(\\mathbf\{u\}\_\{g\}\)\\prod\_\{i=1\}^\{T\}p\(\\bm\{\\delta\}\_\{i\}\)p\(\\mathbf\{u\}\_\{i\}\)\\bigr\)=\(15\)KL\(q\(𝐮g\)∥p\(𝐮g\)\)\+∑i=1TKL\(q\(𝜹i\)∥p\(𝜹i\)\)\+∑i=1TKL\(q\(𝐮i\)∥p\(𝐮i\)\)\.\\displaystyle\\mathrm\{KL\}\\bigl\(q\(\\mathbf\{u\}\_\{g\}\)\\\|p\(\\mathbf\{u\}\_\{g\}\)\\bigr\)\+\\sum\_\{i=1\}^\{T\}\\mathrm\{KL\}\\bigl\(q\(\\bm\{\\delta\}\_\{i\}\)\\\|p\(\\bm\{\\delta\}\_\{i\}\)\\bigr\)\+\\sum\_\{i=1\}^\{T\}\\mathrm\{KL\}\\bigl\(q\(\\mathbf\{u\}\_\{i\}\)\\\|p\(\\mathbf\{u\}\_\{i\}\)\\bigr\)\.

###### Proposition 5\(Decomposability of Expected Log\-Likelihood\)

Under the factorized variational posterior \(Assumption[2](https://arxiv.org/html/2609.19337#Thmassumption2)\), the VFE marginalization of𝐟itot\\mathbf\{f\}\_\{i\}^\{\\mathrm\{tot\}\}, and the additivity of marginal variances \(Proposition[3](https://arxiv.org/html/2609.19337#Thmproposition3)\), the expected log\-likelihood term in the ELBO \([14](https://arxiv.org/html/2609.19337#S3.E14)\) admits a decomposable closed\-form expression\. Equation \([13](https://arxiv.org/html/2609.19337#S3.E13)\) supplies the inducing\-conditional mean model; the expectation below adds posterior uncertainty in the inducing blocks and the sparse GP residual traces:

Define the interpolation matrices𝐀g,i\\mathbf\{A\}\_\{g,i\},𝐀δ,i\\mathbf\{A\}\_\{\\delta,i\}\(sincekδ=ϕ​kgk\_\{\\delta\}=\\phi k\_\{g\},𝐀δ,i=𝐀g,i\\mathbf\{A\}\_\{\\delta,i\}=\\mathbf\{A\}\_\{g,i\}\), and𝐀i\\mathbf\{A\}\_\{i\}from Proposition[2](https://arxiv.org/html/2609.19337#Thmproposition2)\. Then

𝔼q​\[log⁡p⁡\(𝐲i∣𝐮g,𝜹i,𝐮i\)\]\\displaystyle\\mathbb\{E\}\_\{q\}\\bigl\[\\log p\(\\mathbf\{y\}\_\{i\}\\mid\\mathbf\{u\}\_\{g\},\\bm\{\\delta\}\_\{i\},\\mathbf\{u\}\_\{i\}\)\\bigr\]=−ni2​log⁡\(2​π​σ2\)−12​σ2​‖𝐲i−𝝁𝐟i‖2\\displaystyle=\-\\frac\{n\_\{i\}\}\{2\}\\log\(2\\pi\\sigma^\{2\}\)\-\\frac\{1\}\{2\\sigma^\{2\}\}\\\|\\mathbf\{y\}\_\{i\}\-\\bm\{\\mu\}\_\{\\mathbf\{f\}\_\{i\}\}\\\|^\{2\}\(16\)−12​σ2​\(Tr⁡\(𝐀g,i​𝐒g​𝐀g,i⊤\)\+Tr⁡\(𝐀δ,i​𝐒δi​𝐀δ,i⊤\)\+Tr⁡\(𝐀i​𝐒ui​𝐀i⊤\)\)\\displaystyle\-\\frac\{1\}\{2\\sigma^\{2\}\}\\Bigl\(\\operatorname\{Tr\}\\bigl\(\\mathbf\{A\}\_\{g,i\}\\mathbf\{S\}\_\{g\}\\mathbf\{A\}\_\{g,i\}^\{\\top\}\\bigr\)\+\\operatorname\{Tr\}\\bigl\(\\mathbf\{A\}\_\{\\delta,i\}\\mathbf\{S\}\_\{\\delta\_\{i\}\}\\mathbf\{A\}\_\{\\delta,i\}^\{\\top\}\\bigr\)\+\\operatorname\{Tr\}\\bigl\(\\mathbf\{A\}\_\{i\}\\mathbf\{S\}\_\{u\_\{i\}\}\\mathbf\{A\}\_\{i\}^\{\\top\}\\bigr\)\\Bigr\)−12​σ2​Tr⁡\(𝐑i\),\\displaystyle\-\\frac\{1\}\{2\\sigma^\{2\}\}\\operatorname\{Tr\}\\bigl\(\\mathbf\{R\}\_\{i\}\\bigr\),
where the variational predictive mean decomposes as

𝝁𝐟i=𝐀g,i​𝐦g⏟global\+𝐀δ,i​𝐦δi⏟client​deviation\+𝐀i​𝐦ui⏟local​residual,\\bm\{\\mu\}\_\{\\mathbf\{f\}\_\{i\}\}=\\underbrace\{\\mathbf\{A\}\_\{g,i\}\\mathbf\{m\}\_\{g\}\}\_\{\\mathrm\{global\}\}\+\\underbrace\{\\mathbf\{A\}\_\{\\delta,i\}\\mathbf\{m\}\_\{\\delta\_\{i\}\}\}\_\{\\mathrm\{client\}\\ \\mathrm\{deviation\}\}\+\\underbrace\{\\mathbf\{A\}\_\{i\}\\mathbf\{m\}\_\{u\_\{i\}\}\}\_\{\\mathrm\{local\}\\ \\mathrm\{residual\}\},\(17\)and the residual trace decomposes additively via Proposition[3](https://arxiv.org/html/2609.19337#Thmproposition3):

Tr⁡\(𝐑i\)=Tr⁡\(𝐊g,i​i−𝐐g,i​i\)\+Tr⁡\(𝐊δ,i​i−𝐐δ,i​i\)\+Tr⁡\(𝐊i,i​i−𝐐i,i​i\)\.\\operatorname\{Tr\}\\big\(\\mathbf\{R\}\_\{i\}\\big\)=\\operatorname\{Tr\}\(\\mathbf\{K\}\_\{g,ii\}\-\\mathbf\{Q\}\_\{g,ii\}\)\+\\operatorname\{Tr\}\(\\mathbf\{K\}\_\{\\delta,ii\}\-\\mathbf\{Q\}\_\{\\delta,ii\}\)\+\\operatorname\{Tr\}\(\\mathbf\{K\}\_\{i,ii\}\-\\mathbf\{Q\}\_\{i,ii\}\)\.\(18\)

###### Theorem 1\(Explicit ELBO Decomposition\)

Under the mean\-field assumption \(Assumption[2](https://arxiv.org/html/2609.19337#Thmassumption2)\), substitute the closed\-form expected log\-likelihood from Proposition[5](https://arxiv.org/html/2609.19337#Thmproposition5)into the ELBO \([14](https://arxiv.org/html/2609.19337#S3.E14)\)\. The bound admits the following explicit decomposed form, combining the decompositions from Propositions[4](https://arxiv.org/html/2609.19337#Thmproposition4)and[5](https://arxiv.org/html/2609.19337#Thmproposition5):

ℒ=\\displaystyle\\mathcal\{L\}=\{\}∑i=1T\[−ni2log\(2πσ2\)−12​σ2∥𝐲i−𝝁𝐟i∥2\\displaystyle\\sum\_\{i=1\}^\{T\}\\Bigl\[\-\\tfrac\{n\_\{i\}\}\{2\}\\log\(2\\pi\\sigma^\{2\}\)\-\\tfrac\{1\}\{2\\sigma^\{2\}\}\\\|\\mathbf\{y\}\_\{i\}\-\\bm\{\\mu\}\_\{\\mathbf\{f\}\_\{i\}\}\\\|^\{2\}\(19\)−12​σ2​\(Tr⁡\(𝐀g,i​𝐒g​𝐀g,i⊤\)\+Tr⁡\(𝐀δ,i​𝐒δi​𝐀δ,i⊤\)\+Tr⁡\(𝐀i​𝐒ui​𝐀i⊤\)CLOSE\\displaystyle\-\\tfrac\{1\}\{2\\sigma^\{2\}\}\\bigl\(\\operatorname\{Tr\}\(\\mathbf\{A\}\_\{g,i\}\\mathbf\{S\}\_\{g\}\\mathbf\{A\}\_\{g,i\}^\{\\top\}\)\+\\operatorname\{Tr\}\(\\mathbf\{A\}\_\{\\delta,i\}\\mathbf\{S\}\_\{\\delta\_\{i\}\}\\mathbf\{A\}\_\{\\delta,i\}^\{\\top\}\)\+\\operatorname\{Tr\}\(\\mathbf\{A\}\_\{i\}\\mathbf\{S\}\_\{u\_\{i\}\}\\mathbf\{A\}\_\{i\}^\{\\top\}\)\+Tr\(𝐊g,i​i−𝐐g,i​i\)\+Tr\(𝐊δ,i​i−𝐐δ,i​i\)\+Tr\(𝐊i,i​i−𝐐i,i​i\)\)\]\\displaystyle\+\\operatorname\{Tr\}\(\\mathbf\{K\}\_\{g,ii\}\-\\mathbf\{Q\}\_\{g,ii\}\)\+\\operatorname\{Tr\}\(\\mathbf\{K\}\_\{\\delta,ii\}\-\\mathbf\{Q\}\_\{\\delta,ii\}\)\+\\operatorname\{Tr\}\(\\mathbf\{K\}\_\{i,ii\}\-\\mathbf\{Q\}\_\{i,ii\}\)\\bigr\)\\Bigr\]−KL\(q\(𝐮g\)∥p\(𝐮g\)\)−∑i=1TKL\(q\(𝜹i\)∥p\(𝜹i\)\)−∑i=1TKL\(q\(𝐮i\)∥p\(𝐮i\)\),\\displaystyle\-\\mathrm\{KL\}\\bigl\(q\(\\mathbf\{u\}\_\{g\}\)\\\|p\(\\mathbf\{u\}\_\{g\}\)\\bigr\)\-\\sum\_\{i=1\}^\{T\}\\mathrm\{KL\}\\bigl\(q\(\\bm\{\\delta\}\_\{i\}\)\\\|p\(\\bm\{\\delta\}\_\{i\}\)\\bigr\)\-\\sum\_\{i=1\}^\{T\}\\mathrm\{KL\}\\bigl\(q\(\\mathbf\{u\}\_\{i\}\)\\\|p\(\\mathbf\{u\}\_\{i\}\)\\bigr\),
where𝛍𝐟i\\bm\{\\mu\}\_\{\\mathbf\{f\}\_\{i\}\}is the decomposed variational predictive mean defined in \([17](https://arxiv.org/html/2609.19337#S3.E17)\)\.

This result explicitly connects the variational objective to the efficient trace computation derived in Proposition[3](https://arxiv.org/html/2609.19337#Thmproposition3)\. The least squares term fits the decomposed variational predictive mean, the𝐀⋅,i​𝐒⋅​𝐀⋅,i⊤\\mathbf\{A\}\_\{\\cdot,i\}\\mathbf\{S\}\_\{\\cdot\}\\mathbf\{A\}\_\{\\cdot,i\}^\{\\top\}traces propagate posterior uncertainty in each inducing block, and the residual traces penalize sparse\-approximation uncertainty not captured by inducing variables\. Optimization still requires only diagonal elements of𝐑i\\mathbf\{R\}\_\{i\}for the residual terms, preserving the efficiency of our framework\.

##### Summary and intuition \(ELBO\)\.

The explicit ELBO in Theorem[1](https://arxiv.org/html/2609.19337#Thmtheorem1)is the VFE bound in \([14](https://arxiv.org/html/2609.19337#S3.E14)\) with the expected log\-likelihood expanded as in Proposition[5](https://arxiv.org/html/2609.19337#Thmproposition5): a squared\-error term on the hierarchical predictive mean, inducing\-posterior covariance traces, additive sparse GP residual tracesTr⁡\(𝐑i\)\\operatorname\{Tr\}\(\\mathbf\{R\}\_\{i\}\), and KL penalties \(one global term on𝐮g\\mathbf\{u\}\_\{g\}versusTTclient terms on\{𝜹i\}\\\{\\bm\{\\delta\}\_\{i\}\\\}\), encouraging population\-level structure to reside infgf\_\{g\}\.

A detailed derivation is provided in Appendix[C](https://arxiv.org/html/2609.19337#A3)\.

##### Operational identifiability\.

Theorem[1](https://arxiv.org/html/2609.19337#Thmtheorem1)should be read as defining an*ELBO\-regularized allocation*of variation across layers rather than a proof of classical point identification of\(fg,fδ,i\)\(f\_\{g\},f\_\{\\delta,i\}\)from a single client’s marginal law\. Whenkδ=ϕ​kgk\_\{\\delta\}=\\phi k\_\{g\}, the sumfg\+fδ,if\_\{g\}\+f\_\{\\delta,i\}is itself a GP marginally; cross\-client sharing, distinct local kernelskik\_\{i\}, and asymmetric KL penalties on\(𝐮g,\{𝜹i\}\)\(\\mathbf\{u\}\_\{g\},\\\{\\bm\{\\delta\}\_\{i\}\\\}\)nonetheless provide a practical preference for placing multi\-client structure infgf\_\{g\}while reservingfδ,if\_\{\\delta,i\}for stable systematic shifts, as detailed in the paragraph*Identifiability and regularized attribution*in Section[3\.2](https://arxiv.org/html/2609.19337#S3.SS2)\.

### 3\.6pFedHGP Algorithm

Training proceeds in synchronous communication rounds indexed byr=0,1,2,…r=0,1,2,\\dots\. At the start of roundrr, the server holds the current global parameter block

Θglobal\(r\)=\(𝐦g\(r\),𝐒g\(r\),ℓg\(r\),σg2​\(r\),ηϕ\(r\),𝐙g\(r\)\),\\Theta\_\{\\mathrm\{global\}\}^\{\(r\)\}=\\bigl\(\\mathbf\{m\}\_\{g\}^\{\(r\)\},\\mathbf\{S\}\_\{g\}^\{\(r\)\},\\ell\_\{g\}^\{\(r\)\},\\sigma\_\{g\}^\{2\(r\)\},\\eta\_\{\\phi\}^\{\(r\)\},\\mathbf\{Z\}\_\{g\}^\{\(r\)\}\\bigr\),comprising both the global variational parameters\(𝐦g\(r\),𝐒g\(r\)\)\(\\mathbf\{m\}\_\{g\}^\{\(r\)\},\\mathbf\{S\}\_\{g\}^\{\(r\)\}\)and the global hyperparameters\(ℓg\(r\),σg2​\(r\),ηϕ\(r\),𝐙g\(r\)\)\(\\ell\_\{g\}^\{\(r\)\},\\sigma\_\{g\}^\{2\(r\)\},\\eta\_\{\\phi\}^\{\(r\)\},\\mathbf\{Z\}\_\{g\}^\{\(r\)\}\), with deviation scaleϕ\(r\)=softplus⁡\(ηϕ\(r\)\)\\phi^\{\(r\)\}=\\operatorname\{softplus\}\(\\eta\_\{\\phi\}^\{\(r\)\}\)\. Each clientiiholds its local parameter block

Θi,local\(r\)=\(𝐦δi\(r\),𝐒δi\(r\),𝐦ui\(r\),𝐒ui\(r\),ℓi\(r\),σi2​\(r\),𝐙i\(r\)\),\\Theta\_\{i,\\mathrm\{local\}\}^\{\(r\)\}=\\bigl\(\\mathbf\{m\}\_\{\\delta\_\{i\}\}^\{\(r\)\},\\mathbf\{S\}\_\{\\delta\_\{i\}\}^\{\(r\)\},\\mathbf\{m\}\_\{u\_\{i\}\}^\{\(r\)\},\\mathbf\{S\}\_\{u\_\{i\}\}^\{\(r\)\},\\ell\_\{i\}^\{\(r\)\},\\sigma\_\{i\}^\{2\(r\)\},\\mathbf\{Z\}\_\{i\}^\{\(r\)\}\\bigr\),containing local variational parameters\(𝐦δi\(r\),𝐒δi\(r\),𝐦ui\(r\),𝐒ui\(r\)\)\(\\mathbf\{m\}\_\{\\delta\_\{i\}\}^\{\(r\)\},\\mathbf\{S\}\_\{\\delta\_\{i\}\}^\{\(r\)\},\\mathbf\{m\}\_\{u\_\{i\}\}^\{\(r\)\},\\mathbf\{S\}\_\{u\_\{i\}\}^\{\(r\)\}\)and local hyperparameters\(ℓi\(r\),σi2​\(r\),𝐙i\(r\)\)\(\\ell\_\{i\}^\{\(r\)\},\\sigma\_\{i\}^\{2\(r\)\},\\mathbf\{Z\}\_\{i\}^\{\(r\)\}\)\. The algorithm alternates between client\-side local optimization and server\-side global aggregation; broadcasting is integrated at the end of the server update\.

#### 3\.6\.1Local Optimization

When clientiireceivesΘglobal\(r\)\\Theta\_\{\\mathrm\{global\}\}^\{\(r\)\}from the server at roundrr, it treats these global parameters as fixed and optimizes its own local block using its private data\{𝐗i,𝐲i\}\\\{\\mathbf\{X\}\_\{i\},\\mathbf\{y\}\_\{i\}\\\}\. For convenience we collect all client\-specific variables into

Θi,local\(r\)=\(𝐦δi\(r\),𝐒δi\(r\),𝐦ui\(r\),𝐒ui\(r\),ℓi\(r\),σi2​\(r\),𝐙i\(r\)\),\\Theta\_\{i,\\mathrm\{local\}\}^\{\(r\)\}=\\bigl\(\\mathbf\{m\}\_\{\\delta\_\{i\}\}^\{\(r\)\},\\mathbf\{S\}\_\{\\delta\_\{i\}\}^\{\(r\)\},\\mathbf\{m\}\_\{u\_\{i\}\}^\{\(r\)\},\\mathbf\{S\}\_\{u\_\{i\}\}^\{\(r\)\},\\ell\_\{i\}^\{\(r\)\},\\sigma\_\{i\}^\{2\(r\)\},\\mathbf\{Z\}\_\{i\}^\{\(r\)\}\\bigr\),whereq⁡\(𝜹i\)=𝒩⁡\(𝐦δi,𝐒δi\)q\(\\bm\{\\delta\}\_\{i\}\)=\\mathcal\{N\}\(\\mathbf\{m\}\_\{\\delta\_\{i\}\},\\mathbf\{S\}\_\{\\delta\_\{i\}\}\)andq⁡\(𝐮i\)=𝒩⁡\(𝐦ui,𝐒ui\)q\(\\mathbf\{u\}\_\{i\}\)=\\mathcal\{N\}\(\\mathbf\{m\}\_\{u\_\{i\}\},\\mathbf\{S\}\_\{u\_\{i\}\}\)are the local variational distributions, and\(ℓi,σi2,𝐙i\)\(\\ell\_\{i\},\\sigma\_\{i\}^\{2\},\\mathbf\{Z\}\_\{i\}\)are the kernel hyperparameters and inducing inputs of the idiosyncratic local kernel\. The priorp⁡\(𝜹i\)p\(\\bm\{\\delta\}\_\{i\}\)is specified by the deviation prior with hyperparameterϕ\(r\)=softplus⁡\(ηϕ\(r\)\)\\phi^\{\(r\)\}=\\operatorname\{softplus\}\(\\eta\_\{\\phi\}^\{\(r\)\}\), andp⁡\(𝐮i\)p\(\\mathbf\{u\}\_\{i\}\)by the local kernel with parameters\(ℓi,σi2,𝐙i\)\(\\ell\_\{i\},\\sigma\_\{i\}^\{2\},\\mathbf\{Z\}\_\{i\}\)\. In implementation each covariance matrix is represented by its Cholesky factor𝐒\(⋅\)=𝐋\(⋅\)​𝐋\(⋅\)⊤\\mathbf\{S\}\_\{\(\\cdot\)\}=\\mathbf\{L\}\_\{\(\\cdot\)\}\\mathbf\{L\}\_\{\(\\cdot\)\}^\{\\top\}; the optimization is carried out with respect to the entries of these Cholesky factors, although for notational simplicity we keep the covariance notation𝐒\(⋅\)\\mathbf\{S\}\_\{\(\\cdot\)\}\.

###### Proposition 6\(Local ELBO for Client Optimization\)

Based on the explicit ELBO decomposition in Theorem[1](https://arxiv.org/html/2609.19337#Thmtheorem1), clientiimaximizes the local ELBO, which is a function of the local parameters given the fixed global parameters:

ℒi​\(Θi,local∣Θglobal\(r\)\)\\displaystyle\\mathcal\{L\}\_\{i\}\\bigl\(\\Theta\_\{i,\\mathrm\{local\}\}\\mid\\Theta\_\{\\mathrm\{global\}\}^\{\(r\)\}\\bigr\)=−12​σ2​‖𝐲i−𝝁𝐟i‖2\\displaystyle=\-\\tfrac\{1\}\{2\\sigma^\{2\}\}\\bigl\\\|\\mathbf\{y\}\_\{i\}\-\\bm\{\\mu\}\_\{\\mathbf\{f\}\_\{i\}\}\\bigr\\\|^\{2\}\(20\)−12​σ2​\(Tr⁡\(𝐀δ,i​𝐒δi​𝐀δ,i⊤\)\+Tr⁡\(𝐀i​𝐒ui​𝐀i⊤\)\+Tr⁡\(𝐊i,i​i−𝐐i,i​i\)\)\\displaystyle\-\\tfrac\{1\}\{2\\sigma^\{2\}\}\\Bigl\(\\operatorname\{Tr\}\(\\mathbf\{A\}\_\{\\delta,i\}\\mathbf\{S\}\_\{\\delta\_\{i\}\}\\mathbf\{A\}\_\{\\delta,i\}^\{\\top\}\)\+\\operatorname\{Tr\}\(\\mathbf\{A\}\_\{i\}\\mathbf\{S\}\_\{u\_\{i\}\}\\mathbf\{A\}\_\{i\}^\{\\top\}\)\+\\operatorname\{Tr\}\(\\mathbf\{K\}\_\{i,ii\}\-\\mathbf\{Q\}\_\{i,ii\}\)\\Bigr\)−KL\(q\(𝜹i\)∥p\(𝜹i\)\)−KL\(q\(𝐮i\)∥p\(𝐮i\)\),\\displaystyle\-\\mathrm\{KL\}\\bigl\(q\(\\bm\{\\delta\}\_\{i\}\)\\,\\\|\\,p\(\\bm\{\\delta\}\_\{i\}\)\\bigr\)\-\\mathrm\{KL\}\\bigl\(q\(\\mathbf\{u\}\_\{i\}\)\\,\\\|\\,p\(\\mathbf\{u\}\_\{i\}\)\\bigr\),
where the decomposed variational predictive mean𝛍𝐟i\\bm\{\\mu\}\_\{\\mathbf\{f\}\_\{i\}\}is defined in \([17](https://arxiv.org/html/2609.19337#S3.E17)\) \(with𝐦g\\mathbf\{m\}\_\{g\}replaced by𝐦g\(r\)\\mathbf\{m\}\_\{g\}^\{\(r\)\}which is fixed fromΘglobal\(r\)\\Theta\_\{\\mathrm\{global\}\}^\{\(r\)\}\), and the low\-rank approximation matrix𝐐i,i​i=𝐊i​\(𝐗i,𝐙i\)​𝐊i​i−1​𝐊i​\(𝐙i,𝐗i\)\\mathbf\{Q\}\_\{i,ii\}=\\mathbf\{K\}\_\{i\}\(\\mathbf\{X\}\_\{i\},\\mathbf\{Z\}\_\{i\}\)\\mathbf\{K\}\_\{ii\}^\{\-1\}\\mathbf\{K\}\_\{i\}\(\\mathbf\{Z\}\_\{i\},\\mathbf\{X\}\_\{i\}\)depends on local kernel hyperparameters\(ℓi,σi2\)\(\\ell\_\{i\},\\sigma\_\{i\}^\{2\}\)and local inducing locations𝐙i\\mathbf\{Z\}\_\{i\}\. Constant terms with respect to local parameters \(includingTr⁡\(𝐀g,i​𝐒g\(r\)​𝐀g,i⊤\)\\operatorname\{Tr\}\(\\mathbf\{A\}\_\{g,i\}\\mathbf\{S\}\_\{g\}^\{\(r\)\}\\mathbf\{A\}\_\{g,i\}^\{\\top\}\),Tr⁡\(𝐊g,i​i−𝐐g,i​i\)\\operatorname\{Tr\}\(\\mathbf\{K\}\_\{g,ii\}\-\\mathbf\{Q\}\_\{g,ii\}\), andTr⁡\(𝐊δ,i​i−𝐐δ,i​i\)\\operatorname\{Tr\}\(\\mathbf\{K\}\_\{\\delta,ii\}\-\\mathbf\{Q\}\_\{\\delta,ii\}\), which depend only onΘglobal\(r\)\\Theta\_\{\\mathrm\{global\}\}^\{\(r\)\}\) have been omitted\. The notationℒi​\(Θi,local∣Θglobal\(r\)\)\\mathcal\{L\}\_\{i\}\(\\Theta\_\{i,\\mathrm\{local\}\}\\mid\\Theta\_\{\\mathrm\{global\}\}^\{\(r\)\}\)emphasizes that the global parameters are given/fixed, and the optimization is performed only over the local parametersΘi,local\\Theta\_\{i,\\mathrm\{local\}\}\.

The explicit decomposition in \([20](https://arxiv.org/html/2609.19337#S3.E20)\) reveals that the local objective separates into: \(1\) a data\-fitting least squares term−12​σ2​‖𝐲i−𝝁𝐟i‖2\-\\frac\{1\}\{2\\sigma^\{2\}\}\\\|\\mathbf\{y\}\_\{i\}\-\\bm\{\\mu\}\_\{\\mathbf\{f\}\_\{i\}\}\\\|^\{2\}involving the decomposed predictive mean𝝁𝐟i\\bm\{\\mu\}\_\{\\mathbf\{f\}\_\{i\}\}, \(2\) a trace regularization termTr⁡\(𝐊i,i​i−𝐐i,i​i\)\\operatorname\{Tr\}\(\\mathbf\{K\}\_\{i,ii\}\-\\mathbf\{Q\}\_\{i,ii\}\)for the local residual covariance, and \(3\) KL regularization terms for the client\-specific deviation and local components\. SinceΘglobal\(r\)\\Theta\_\{\\mathrm\{global\}\}^\{\(r\)\}is fixed, the global contributions𝐀g,i​𝐦g\(r\)\\mathbf\{A\}\_\{g,i\}\\mathbf\{m\}\_\{g\}^\{\(r\)\}in the predictive mean can be computed using only local data and the known global posterior approximation\.

Clientiiapplies a few steps of gradient ascent to \([20](https://arxiv.org/html/2609.19337#S3.E20)\) with respect toΘi,local\\Theta\_\{i,\\mathrm\{local\}\}, while holdingΘglobal\(r\)\\Theta\_\{\\mathrm\{global\}\}^\{\(r\)\}fixed\. This local optimization simultaneously updates the variational means𝐦δi\\mathbf\{m\}\_\{\\delta\_\{i\}\}and𝐦ui\\mathbf\{m\}\_\{u\_\{i\}\}, the Cholesky factors underlying𝐒δi\\mathbf\{S\}\_\{\\delta\_\{i\}\}and𝐒ui\\mathbf\{S\}\_\{u\_\{i\}\}, the local kernel hyperparametersℓi\\ell\_\{i\}\(lengthscale / bandwidth\) andσi2\\sigma\_\{i\}^\{2\}\(output scale\), and the local inducing locations𝐙i\\mathbf\{Z\}\_\{i\}\. We denote the updated local parameters byΘi,local\(r\+1\)\\Theta\_\{i,\\mathrm\{local\}\}^\{\(r\+1\)\}\.

#### 3\.6\.2Global Component Optimization

In this step we optimize the global parameter block

Θglobal=\(𝐦g,𝐒g,ℓg,σg2,ηϕ,𝐙g\),ϕ=softplus⁡\(ηϕ\),\\Theta\_\{\\mathrm\{global\}\}=\\bigl\(\\mathbf\{m\}\_\{g\},\\mathbf\{S\}\_\{g\},\\ell\_\{g\},\\sigma\_\{g\}^\{2\},\\eta\_\{\\phi\},\\mathbf\{Z\}\_\{g\}\\bigr\),\\qquad\\phi=\\operatorname\{softplus\}\(\\eta\_\{\\phi\}\),while keeping all local blocks\{Θi,local\(r\+1\)\}i=1T\\\{\\Theta\_\{i,\\mathrm\{local\}\}^\{\(r\+1\)\}\\\}\_\{i=1\}^\{T\}fixed at the values obtained from the client\-side updates in the current round\. We first isolate the part of the ELBO that depends onΘglobal\\Theta\_\{\\mathrm\{global\}\}and then show that its gradient decomposes into client\-wise contributions that can be computed locally and aggregated at the server\.

###### Proposition 7\(Global ELBO given fixed local blocks\)

Let\{Θi,local\(r\+1\)\}i=1T\\\{\\Theta\_\{i,\\mathrm\{local\}\}^\{\(r\+1\)\}\\\}\_\{i=1\}^\{T\}be fixed local parameters after the client\-side updates in roundrr, and let𝛍𝐟i\\bm\{\\mu\}\_\{\\mathbf\{f\}\_\{i\}\}be the decomposed predictive mean in \([17](https://arxiv.org/html/2609.19337#S3.E17)\)\. Up to an additive constant independent ofΘglobal\\Theta\_\{\\mathrm\{global\}\}, the ELBO \([19](https://arxiv.org/html/2609.19337#S3.E19)\) can be written as

ℒ⁡\(Θglobal∣\{Θi,local\(r\+1\)\}\)\\displaystyle\\mathcal\{L\}\\bigl\(\\Theta\_\{\\mathrm\{global\}\}\\mid\\\{\\Theta\_\{i,\\mathrm\{local\}\}^\{\(r\+1\)\}\\\}\\bigr\)=∑i=1Tℒiglobal​\(Θglobal∣Θi,local\(r\+1\)\)\\displaystyle=\\sum\_\{i=1\}^\{T\}\\mathcal\{L\}\_\{i\}^\{\\mathrm\{global\}\}\\bigl\(\\Theta\_\{\\mathrm\{global\}\}\\mid\\Theta\_\{i,\\mathrm\{local\}\}^\{\(r\+1\)\}\\bigr\)\(21\)−KL\(q\(𝐮g\)∥p\(𝐮g\)\)−∑i=1TKL\(q\(𝜹i\)∥p\(𝜹i;ϕ,𝐊g​g\)\)\+const,\\displaystyle\-\\mathrm\{KL\}\\bigl\(q\(\\mathbf\{u\}\_\{g\}\)\\,\\\|\\,p\(\\mathbf\{u\}\_\{g\}\)\\bigr\)\-\\sum\_\{i=1\}^\{T\}\\mathrm\{KL\}\\bigl\(q\(\\bm\{\\delta\}\_\{i\}\)\\,\\\|\\,p\(\\bm\{\\delta\}\_\{i\};\\phi,\\mathbf\{K\}\_\{gg\}\)\\bigr\)\+\\textnormal\{const\},where the client\-wise global contributions are

ℒiglobal​\(Θglobal∣Θi,local\(r\+1\)\)\\displaystyle\\mathcal\{L\}\_\{i\}^\{\\mathrm\{global\}\}\\bigl\(\\Theta\_\{\\mathrm\{global\}\}\\mid\\Theta\_\{i,\\mathrm\{local\}\}^\{\(r\+1\)\}\\bigr\)=−12​σ2​‖𝐲i−𝝁𝐟i‖2\\displaystyle=\-\\tfrac\{1\}\{2\\sigma^\{2\}\}\\bigl\\\|\\mathbf\{y\}\_\{i\}\-\\bm\{\\mu\}\_\{\\mathbf\{f\}\_\{i\}\}\\bigr\\\|^\{2\}\(22\)−12​σ2​\(Tr⁡\(𝐀g,i​𝐒g​𝐀g,i⊤\)\+Tr⁡\(𝐀δ,i​𝐒δi​𝐀δ,i⊤\)CLOSE\\displaystyle\-\\tfrac\{1\}\{2\\sigma^\{2\}\}\\Bigl\(\\operatorname\{Tr\}\(\\mathbf\{A\}\_\{g,i\}\\mathbf\{S\}\_\{g\}\\mathbf\{A\}\_\{g,i\}^\{\\top\}\)\+\\operatorname\{Tr\}\(\\mathbf\{A\}\_\{\\delta,i\}\\mathbf\{S\}\_\{\\delta\_\{i\}\}\\mathbf\{A\}\_\{\\delta,i\}^\{\\top\}\)OPEN\+Tr⁡\(𝐊g,i​i−𝐐g,i​i\)\+Tr⁡\(𝐊δ,i​i−𝐐δ,i​i\)\)\.\\displaystyle\+\\operatorname\{Tr\}\(\\mathbf\{K\}\_\{g,ii\}\-\\mathbf\{Q\}\_\{g,ii\}\)\+\\operatorname\{Tr\}\(\\mathbf\{K\}\_\{\\delta,ii\}\-\\mathbf\{Q\}\_\{\\delta,ii\}\)\\Bigr\)\.Here𝐊g,i​i\\mathbf\{K\}\_\{g,ii\},𝐐g,i​i\\mathbf\{Q\}\_\{g,ii\},𝐊δ,i​i\\mathbf\{K\}\_\{\\delta,ii\}, and𝐐δ,i​i\\mathbf\{Q\}\_\{\\delta,ii\}depend onΘglobal\\Theta\_\{\\mathrm\{global\}\}through the global kernel hyperparameters\(ℓg,σg2\)\(\\ell\_\{g\},\\sigma\_\{g\}^\{2\}\),ϕ=softplus⁡\(ηϕ\)\\phi=\\operatorname\{softplus\}\(\\eta\_\{\\phi\}\), and inducing locations𝐙g\\mathbf\{Z\}\_\{g\}; the predictive mean𝛍𝐟i\\bm\{\\mu\}\_\{\\mathbf\{f\}\_\{i\}\}depends onΘglobal\\Theta\_\{\\mathrm\{global\}\}through\(𝐦g,𝐒g\)\(\\mathbf\{m\}\_\{g\},\\mathbf\{S\}\_\{g\}\)and the global kernel; and𝐀g,i\\mathbf\{A\}\_\{g,i\},𝐀δ,i\\mathbf\{A\}\_\{\\delta,i\}inherit the same dependence through𝐙g\\mathbf\{Z\}\_\{g\}andkgk\_\{g\}\. The termsTr⁡\(𝐀i​𝐒ui​𝐀i⊤\)\\operatorname\{Tr\}\(\\mathbf\{A\}\_\{i\}\\mathbf\{S\}\_\{u\_\{i\}\}\\mathbf\{A\}\_\{i\}^\{\\top\}\),Tr⁡\(𝐊i,i​i−𝐐i,i​i\)\\operatorname\{Tr\}\(\\mathbf\{K\}\_\{i,ii\}\-\\mathbf\{Q\}\_\{i,ii\}\), and the local KL on𝐮i\\mathbf\{u\}\_\{i\}do not depend onΘglobal\\Theta\_\{\\mathrm\{global\}\}and are absorbed into the constant\. The deviation KL terms are retained in \([21](https://arxiv.org/html/2609.19337#S3.E21)\) becausep⁡\(𝛅i\)p\(\\bm\{\\delta\}\_\{i\}\)depends on\(ϕ,𝐊g​g\)\(\\phi,\\mathbf\{K\}\_\{gg\}\)\. A proof is provided in Appendix[F](https://arxiv.org/html/2609.19337#A6)\.

Proposition[7](https://arxiv.org/html/2609.19337#Thmproposition7)shows that optimizing the global block amounts to maximizing a sum of client\-wise objectivesℒiglobal\\mathcal\{L\}\_\{i\}^\{\\mathrm\{global\}\}minus client\-specific deviation KL penalties, together with the global inducing KL on𝐮g\\mathbf\{u\}\_\{g\}\. This structure enables a federated scheme in which each client transmits a gradient summary that already includes its private deviation\-KL contribution, while the server aggregates these summaries and subtracts only the shared global KL term \(which depends on\(𝐦g,𝐒g\)\(\\mathbf\{m\}\_\{g\},\\mathbf\{S\}\_\{g\}\)held at the server\)\.

For fixed local blocks\{Θi,local\(r\+1\)\}\\\{\\Theta\_\{i,\\mathrm\{local\}\}^\{\(r\+1\)\}\\\}, each clientiicomputes and transmits

𝒮~i\(r\+1\):=∇Θglobal\[ℒiglobal\(Θglobal∣Θi,local\(r\+1\)\)−KL\(q\(𝜹i\)∥p\(𝜹i;ϕ,𝐊g​g\)\)\]\|Θglobal=Θglobal\(r\),\\widetilde\{\\mathcal\{S\}\}\_\{i\}^\{\(r\+1\)\}:=\\nabla\_\{\\Theta\_\{\\mathrm\{global\}\}\}\\Bigl\[\\mathcal\{L\}\_\{i\}^\{\\mathrm\{global\}\}\\bigl\(\\Theta\_\{\\mathrm\{global\}\}\\mid\\Theta\_\{i,\\mathrm\{local\}\}^\{\(r\+1\)\}\\bigr\)\-\\mathrm\{KL\}\\bigl\(q\(\\bm\{\\delta\}\_\{i\}\)\\,\\\|\\,p\(\\bm\{\\delta\}\_\{i\};\\phi,\\mathbf\{K\}\_\{gg\}\)\\bigr\)\\Bigr\]\\Big\|\_\{\\Theta\_\{\\mathrm\{global\}\}=\\Theta\_\{\\mathrm\{global\}\}^\{\(r\)\}\},\(23\)using its private observations\(𝐗i,𝐲i\)\(\\mathbf\{X\}\_\{i\},\\mathbf\{y\}\_\{i\}\), updated local blockΘi,local\(r\+1\)\\Theta\_\{i,\\mathrm\{local\}\}^\{\(r\+1\)\}\(which containsq⁡\(𝜹i\)q\(\\bm\{\\delta\}\_\{i\}\)\), and broadcastΘglobal\(r\)\\Theta\_\{\\mathrm\{global\}\}^\{\(r\)\}\. Only the vector𝒮~i\(r\+1\)\\widetilde\{\\mathcal\{S\}\}\_\{i\}^\{\(r\+1\)\}is sent to the server;𝐦δi\\mathbf\{m\}\_\{\\delta\_\{i\}\},𝐒δi\\mathbf\{S\}\_\{\\delta\_\{i\}\}, and raw data remain local\. The full global gradient is then

∇Θglobalℒ\(Θglobal∣\{Θi,local\(r\+1\)\}\)=∑i=1T𝒮~i\(r\+1\)−∇ΘglobalKL\(q\(𝐮g\)∥p\(𝐮g\)\),\\nabla\_\{\\Theta\_\{\\mathrm\{global\}\}\}\\mathcal\{L\}\\bigl\(\\Theta\_\{\\mathrm\{global\}\}\\mid\\\{\\Theta\_\{i,\\mathrm\{local\}\}^\{\(r\+1\)\}\\\}\\bigr\)=\\sum\_\{i=1\}^\{T\}\\widetilde\{\\mathcal\{S\}\}\_\{i\}^\{\(r\+1\)\}\-\\nabla\_\{\\Theta\_\{\\mathrm\{global\}\}\}\\,\\mathrm\{KL\}\\bigl\(q\(\\mathbf\{u\}\_\{g\}\)\\,\\\|\\,p\(\\mathbf\{u\}\_\{g\}\)\\bigr\),\(24\)
evaluated atΘglobal=Θglobal\(r\)\\Theta\_\{\\mathrm\{global\}\}=\\Theta\_\{\\mathrm\{global\}\}^\{\(r\)\}\.

Upon collecting\{𝒮~i\(r\+1\)\}i=1T\\\{\\widetilde\{\\mathcal\{S\}\}\_\{i\}^\{\(r\+1\)\}\\\}\_\{i=1\}^\{T\}, the server aggregates

𝒮~\(r\+1\):=∑i=1T𝒮~i\(r\+1\),\\widetilde\{\\mathcal\{S\}\}^\{\(r\+1\)\}:=\\sum\_\{i=1\}^\{T\}\\widetilde\{\\mathcal\{S\}\}\_\{i\}^\{\(r\+1\)\},forms∇Θglobalℒ\\nabla\_\{\\Theta\_\{\\mathrm\{global\}\}\}\\mathcal\{L\}using \([24](https://arxiv.org/html/2609.19337#S3.E24)\) \(subtracting only the global inducing KL, which is evaluated from\(𝐦g,𝐒g\)\(\\mathbf\{m\}\_\{g\},\\mathbf\{S\}\_\{g\}\)at the server\), and updates the entire global block by a gradient step, e\.g\.

Θglobal\(r\+1\)=Θglobal\(r\)\+ηglobal​∇Θglobalℒ,\\Theta\_\{\\mathrm\{global\}\}^\{\(r\+1\)\}=\\Theta\_\{\\mathrm\{global\}\}^\{\(r\)\}\+\\eta\_\{\\mathrm\{global\}\}\\,\\nabla\_\{\\Theta\_\{\\mathrm\{global\}\}\}\\mathcal\{L\},\(25\)where the gradient in \([25](https://arxiv.org/html/2609.19337#S3.E25)\) is taken with respect to the unconstrained global coordinates inΘglobal\(r\)\\Theta\_\{\\mathrm\{global\}\}^\{\(r\)\}\(includingηϕ\(r\)\\eta\_\{\\phi\}^\{\(r\)\}rather thanϕ\\phidirectly\)\. Positive hyperparameters such asϕ\\phienter the ELBO only through reparameterizationsϕ=softplus⁡\(ηϕ\)\\phi=\\operatorname\{softplus\}\(\\eta\_\{\\phi\}\)and analogous softplus maps for lengthscales and variances \(Section[3\.7](https://arxiv.org/html/2609.19337#S3.SS7)\)\. The updatedΘglobal\(r\+1\)\\Theta\_\{\\mathrm\{global\}\}^\{\(r\+1\)\}is then broadcast to clients for the next communication round\. This procedure shows explicitly that the global component optimization is decomposable across clients and implementable in a federated way via local gradients and server\-side aggregation\.

##### Deviation scaleϕ\\phias a federated hyperparameter\.

The scale factorϕ\>0\\phi\>0is a globally shared hyperparameter that controls the prior variance of the systematic deviation layer,kδ=ϕ​kgk\_\{\\delta\}=\\phi k\_\{g\}\. It should be interpreted as a population\-level regularization scale rather than as a constraint that all clients have identical realized deviation magnitudes\. Client\-specific amplitudes and shapes are still learned through the local variational posteriorsq⁡\(𝜹i\)=𝒩⁡\(𝐦δi,𝐒δi\)q\(\\bm\{\\delta\}\_\{i\}\)=\\mathcal\{N\}\(\\mathbf\{m\}\_\{\\delta\_\{i\}\},\\mathbf\{S\}\_\{\\delta\_\{i\}\}\)\. To enforceϕ\>0\\phi\>0throughout optimization, we include an unconstrained coordinateηϕ∈ℝ\\eta\_\{\\phi\}\\in\\mathbb\{R\}inΘglobal\\Theta\_\{\\mathrm\{global\}\}and setϕ=softplus⁡\(ηϕ\)\\phi=\\operatorname\{softplus\}\(\\eta\_\{\\phi\}\)\. After the client\-side updates in roundr\+1r\+1, the server aggregates client\-wise contributions and updates

ηϕ\(r\+1\)=ηϕ\(r\)\+ηglobal​∑i=1T𝒮i,ηϕ\(r\+1\),ϕ\(r\+1\)=softplus⁡\(ηϕ\(r\+1\)\),\\eta\_\{\\phi\}^\{\(r\+1\)\}=\\eta\_\{\\phi\}^\{\(r\)\}\+\\eta\_\{\\mathrm\{global\}\}\\sum\_\{i=1\}^\{T\}\\mathcal\{S\}\_\{i,\\eta\_\{\\phi\}\}^\{\(r\+1\)\},\\qquad\\phi^\{\(r\+1\)\}=\\operatorname\{softplus\}\\\!\\bigl\(\\eta\_\{\\phi\}^\{\(r\+1\)\}\\bigr\),\(26\)where clientiicontributes

𝒮i,ηϕ\(r\+1\)=∂ℒiglobal∂ηϕ−∂∂ηϕKL\(q\(𝜹i\)∥p\(𝜹i;ϕ,𝐊g​g\)\)\|Θglobal=Θglobal\(r\),\\mathcal\{S\}\_\{i,\\eta\_\{\\phi\}\}^\{\(r\+1\)\}=\\left\.\\frac\{\\partial\\mathcal\{L\}\_\{i\}^\{\\mathrm\{global\}\}\}\{\\partial\\eta\_\{\\phi\}\}\-\\frac\{\\partial\}\{\\partial\\eta\_\{\\phi\}\}\\mathrm\{KL\}\\bigl\(q\(\\bm\{\\delta\}\_\{i\}\)\\,\\\|\\,p\(\\bm\{\\delta\}\_\{i\};\\phi,\\mathbf\{K\}\_\{gg\}\)\\bigr\)\\right\|\_\{\\Theta\_\{\\mathrm\{global\}\}=\\Theta\_\{\\mathrm\{global\}\}^\{\(r\)\}\},\(27\)which is theηϕ\\eta\_\{\\phi\}\-component of𝒮~i\(r\+1\)\\widetilde\{\\mathcal\{S\}\}\_\{i\}^\{\(r\+1\)\}in \([23](https://arxiv.org/html/2609.19337#S3.E23)\)\. Under the shared\-basis parameterizationkδ=ϕ​kgk\_\{\\delta\}=\\phi k\_\{g\}, the deviation interpolation matrix is𝐀δ,i=ϕ​𝐊g​\(𝐗i,𝐙g\)​\(ϕ​𝐊g​g\)−1=𝐀g,i\\mathbf\{A\}\_\{\\delta,i\}=\\phi\\mathbf\{K\}\_\{g\}\(\\mathbf\{X\}\_\{i\},\\mathbf\{Z\}\_\{g\}\)\(\\phi\\mathbf\{K\}\_\{gg\}\)^\{\-1\}=\\mathbf\{A\}\_\{g,i\}\(Proposition[2](https://arxiv.org/html/2609.19337#Thmproposition2)\), so withq⁡\(𝜹i\)q\(\\bm\{\\delta\}\_\{i\}\)fixed at the global step the deviation contribution𝐀δ,i​𝐦δi\\mathbf\{A\}\_\{\\delta,i\}\\mathbf\{m\}\_\{\\delta\_\{i\}\}to𝝁𝐟i\\bm\{\\mu\}\_\{\\mathbf\{f\}\_\{i\}\}has no directϕ\\phidependence\. The squared\-error term−12​σ2​‖𝐲i−𝝁𝐟i‖2\-\\frac\{1\}\{2\\sigma^\{2\}\}\\\|\\mathbf\{y\}\_\{i\}\-\\bm\{\\mu\}\_\{\\mathbf\{f\}\_\{i\}\}\\\|^\{2\}and the inducing traceTr⁡\(𝐀δ,i​𝐒δi​𝐀δ,i⊤\)\\operatorname\{Tr\}\(\\mathbf\{A\}\_\{\\delta,i\}\\mathbf\{S\}\_\{\\delta\_\{i\}\}\\mathbf\{A\}\_\{\\delta,i\}^\{\\top\}\)inℒiglobal\\mathcal\{L\}\_\{i\}^\{\\mathrm\{global\}\}therefore contribute zero directϕ\\phi\-derivatives in the standard scalar model unless the implementation parameterizes means or loadings differently\. By the chain rule,

∂ℒiglobal∂ϕ−∂∂ϕKL\(q\(𝜹i\)∥p\(𝜹i;ϕ,𝐊g​g\)\)\|Θglobal=Θglobal\(r\)\\displaystyle\\left\.\\frac\{\\partial\\mathcal\{L\}\_\{i\}^\{\\mathrm\{global\}\}\}\{\\partial\\phi\}\-\\frac\{\\partial\}\{\\partial\\phi\}\\mathrm\{KL\}\\bigl\(q\(\\bm\{\\delta\}\_\{i\}\)\\,\\\|\\,p\(\\bm\{\\delta\}\_\{i\};\\phi,\\mathbf\{K\}\_\{gg\}\)\\bigr\)\\right\|\_\{\\Theta\_\{\\mathrm\{global\}\}=\\Theta\_\{\\mathrm\{global\}\}^\{\(r\)\}\}\(28\)=∂∂ϕ\[−12​σ2Tr\(𝐊δ,i​i−𝐐δ,i​i\)\\displaystyle=\\left\.\\frac\{\\partial\}\{\\partial\\phi\}\\Bigl\[\-\\tfrac\{1\}\{2\\sigma^\{2\}\}\\operatorname\{Tr\}\(\\mathbf\{K\}\_\{\\delta,ii\}\-\\mathbf\{Q\}\_\{\\delta,ii\}\)\\right\.−KL\(q\(𝜹i\)∥p\(𝜹i;ϕ,𝐊g​g\)\)\]\|Θglobal=Θglobal\(r\),\\displaystyle\\left\.\-\\mathrm\{KL\}\\bigl\(q\(\\bm\{\\delta\}\_\{i\}\)\\,\\\|\\,p\(\\bm\{\\delta\}\_\{i\};\\phi,\\mathbf\{K\}\_\{gg\}\)\\bigr\)\\Bigr\]\\right\|\_\{\\Theta\_\{\\mathrm\{global\}\}=\\Theta\_\{\\mathrm\{global\}\}^\{\(r\)\}\},so that

𝒮i,ηϕ\(r\+1\)\\displaystyle\\mathcal\{S\}\_\{i,\\eta\_\{\\phi\}\}^\{\(r\+1\)\}=\(∂ℒiglobal∂ϕ−∂∂ϕKL\(q\(𝜹i\)∥p\(𝜹i;ϕ,𝐊g​g\)\)\|Θglobal=Θglobal\(r\)\)∂ϕ∂ηϕ,\\displaystyle=\\left\(\\left\.\\frac\{\\partial\\mathcal\{L\}\_\{i\}^\{\\mathrm\{global\}\}\}\{\\partial\\phi\}\-\\frac\{\\partial\}\{\\partial\\phi\}\\mathrm\{KL\}\\bigl\(q\(\\bm\{\\delta\}\_\{i\}\)\\,\\\|\\,p\(\\bm\{\\delta\}\_\{i\};\\phi,\\mathbf\{K\}\_\{gg\}\)\\bigr\)\\right\|\_\{\\Theta\_\{\\mathrm\{global\}\}=\\Theta\_\{\\mathrm\{global\}\}^\{\(r\)\}\}\\right\)\\frac\{\\partial\\phi\}\{\\partial\\eta\_\{\\phi\}\},\(29\)∂ϕ∂ηϕ\\displaystyle\\frac\{\\partial\\phi\}\{\\partial\\eta\_\{\\phi\}\}=σ⁡\(ηϕ\)∈\(0,1\),\\displaystyle=\\sigma\(\\eta\_\{\\phi\}\)\\in\(0,1\),withσ⁡\(⋅\)\\sigma\(\\cdot\)the logistic function\. Under the shared\-basis parameterizationkδ=ϕ​kgk\_\{\\delta\}=\\phi k\_\{g\},ϕ\\phiaffects the global step primarily through the deviation residual traceTr⁡\(𝐊δ,i​i−𝐐δ,i​i\)\\operatorname\{Tr\}\(\\mathbf\{K\}\_\{\\delta,ii\}\-\\mathbf\{Q\}\_\{\\delta,ii\}\)and the deviation prior/KL term in \([28](https://arxiv.org/html/2609.19337#S3.E28)\)\. It can influence the fitted deviation mean only indirectly through subsequent local variational updates, becauseq⁡\(𝜹i\)q\(\\bm\{\\delta\}\_\{i\}\)is reoptimized under the updated prior scale afterϕ\\phiis broadcast\. In practice these gradients are computed by automatic differentiation within the federated block\-coordinate loop\. A more flexible extension would replace the sharedϕ\\phiwith client\-specific scalesϕi\\phi\_\{i\}or a hierarchical prior over\{ϕi\}\\\{\\phi\_\{i\}\\\}; we leave this extension for future work to avoid increasing communication and server\-side hyperparameter complexity\.

Algorithm 1pFedHGP federated variational optimization \(one communication round\)\.1:Server input:

Θglobal\(r\)\\Theta\_\{\\mathrm\{global\}\}^\{\(r\)\};Clientsi=1,…,Ti=1,\\ldots,T:private data

\(𝐗i,𝐲i\)\(\\mathbf\{X\}\_\{i\},\\mathbf\{y\}\_\{i\}\), local blocks

Θi,local\(r\)\\Theta\_\{i,\\mathrm\{local\}\}^\{\(r\)\}
2:Broadcast

Θglobal\(r\)\\Theta\_\{\\mathrm\{global\}\}^\{\(r\)\}from server to all clients

3:foreach client

iiin paralleldo

4:Maximize local ELBO \([20](https://arxiv.org/html/2609.19337#S3.E20)\) over

Θi,local\\Theta\_\{i,\\mathrm\{local\}\}with

Θglobal\(r\)\\Theta\_\{\\mathrm\{global\}\}^\{\(r\)\}held fixed; obtain

Θi,local\(r\+1\)\\Theta\_\{i,\\mathrm\{local\}\}^\{\(r\+1\)\}
5:Compute

𝒮~i\(r\+1\)\\widetilde\{\\mathcal\{S\}\}\_\{i\}^\{\(r\+1\)\}from \([23](https://arxiv.org/html/2609.19337#S3.E23)\) using

\(𝐗i,𝐲i\)\(\\mathbf\{X\}\_\{i\},\\mathbf\{y\}\_\{i\}\),

Θi,local\(r\+1\)\\Theta\_\{i,\\mathrm\{local\}\}^\{\(r\+1\)\}, and

Θglobal\(r\)\\Theta\_\{\\mathrm\{global\}\}^\{\(r\)\}
6:Transmit

𝒮~i\(r\+1\)\\widetilde\{\\mathcal\{S\}\}\_\{i\}^\{\(r\+1\)\}to server \(no raw data or local variational parameters\)

7:endfor

8:Server aggregation:

𝒮~\(r\+1\)←∑i=1T𝒮~i\(r\+1\)\\widetilde\{\\mathcal\{S\}\}^\{\(r\+1\)\}\\leftarrow\\sum\_\{i=1\}^\{T\}\\widetilde\{\\mathcal\{S\}\}\_\{i\}^\{\(r\+1\)\}; form

∇Θglobalℒ\\nabla\_\{\\Theta\_\{\\mathrm\{global\}\}\}\\mathcal\{L\}via \([24](https://arxiv.org/html/2609.19337#S3.E24)\)

9:Global update:

Θglobal\(r\+1\)←Θglobal\(r\)\+ηglobal​∇Θglobalℒ\\Theta\_\{\\mathrm\{global\}\}^\{\(r\+1\)\}\\leftarrow\\Theta\_\{\\mathrm\{global\}\}^\{\(r\)\}\+\\eta\_\{\\mathrm\{global\}\}\\nabla\_\{\\Theta\_\{\\mathrm\{global\}\}\}\\mathcal\{L\}\(cf\. \([25](https://arxiv.org/html/2609.19337#S3.E25)\)\)

10:Broadcast

Θglobal\(r\+1\)\\Theta\_\{\\mathrm\{global\}\}^\{\(r\+1\)\}to clients for round

r\+1r\+1

Each round therefore alternates a*local*block update with a*global*aggregation step; within the local step the global parameters act as fixed conditioning values for the client sub\-problem \(Proposition[6](https://arxiv.org/html/2609.19337#Thmproposition6)\), while the transmitted gradient𝒮~i\(r\+1\)\\widetilde\{\\mathcal\{S\}\}\_\{i\}^\{\(r\+1\)\}couples updated local states back to the global block for the subsequent server update\.

##### Summary and intuition \(federated optimization\)\.

Algorithm[1](https://arxiv.org/html/2609.19337#alg1)implements block coordinate descent: clients first maximize their local ELBO blocks withΘglobal\\Theta\_\{\\mathrm\{global\}\}fixed, then transmit low\-dimensional gradient summaries𝒮~i\\widetilde\{\\mathcal\{S\}\}\_\{i\}that already include each client’s deviation\-KL contribution, so the server updates the shared backbone—includingηϕ\\eta\_\{\\phi\}\(henceϕ=softplus⁡\(ηϕ\)\\phi=\\operatorname\{softplus\}\(\\eta\_\{\\phi\}\)\) and𝐙g\\mathbf\{Z\}\_\{g\}—using only aggregated gradients and the global inducing variational parameters, without pooling raw observations or local deviation posteriors\.

### 3\.7Implementation Details

##### Covariance factorization\.

We use Cholesky decomposition to factorize covariance matrices for numerical stability and to avoid explicit matrix inversion during variational updates\. Specifically, for any covariance matrix𝐒\\mathbf\{S\}, we compute the lower triangular Cholesky factor𝐋\\mathbf\{L\}such that𝐒=𝐋𝐋⊤\\mathbf\{S\}=\\mathbf\{L\}\\mathbf\{L\}^\{\\top\}\. All covariance matrices are regularized with a fixed jitterϵ=10−6\\epsilon=10^\{\-6\}before decomposition, i\.e\., we compute the Cholesky factorization of𝐊\+ϵ​𝐈\\mathbf\{K\}\+\\epsilon\\mathbf\{I\}where𝐊\\mathbf\{K\}is the original covariance matrix\. This ensures numerical stability while maintaining the positive definiteness of the matrices\.

##### Whitened variational parameterization\.

For numerical stability, we utilize the whitened variational strategy parameterization\([Matthews et al\., 2017](https://arxiv.org/html/2609.19337#bib.bib21)\)\. Given a prior distributionp⁡\(𝐮\)=𝒩⁡\(𝟎,𝐊\)p\(\\mathbf\{u\}\)=\\mathcal\{N\}\(\\mathbf\{0\},\\mathbf\{K\}\)with Cholesky factorization𝐊=𝐋𝐋⊤\\mathbf\{K\}=\\mathbf\{L\}\\mathbf\{L\}^\{\\top\}, we parameterize the inducing variables as𝐮=𝐋𝐯\\mathbf\{u\}=\\mathbf\{L\}\\mathbf\{v\}, where𝐯\\mathbf\{v\}are the whitened variational parameters\. The variational distribution is then expressed asq⁡\(𝐯\)=𝒩⁡\(𝐦~,𝐒~\)q\(\\mathbf\{v\}\)=\\mathcal\{N\}\(\\tilde\{\\mathbf\{m\}\},\\tilde\{\\mathbf\{S\}\}\)with whitened mean𝐦~\\tilde\{\\mathbf\{m\}\}and covariance𝐒~\\tilde\{\\mathbf\{S\}\}\. The original variational parameters are recovered via𝐦=𝐋​𝐦~\\mathbf\{m\}=\\mathbf\{L\}\\tilde\{\\mathbf\{m\}\}and𝐒=𝐋​𝐒~​𝐋⊤\\mathbf\{S\}=\\mathbf\{L\}\\tilde\{\\mathbf\{S\}\}\\mathbf\{L\}^\{\\top\}\. This parameterization improves optimization by removing correlations in the variational posterior, leading to better\-conditioned optimization problems and more stable gradient computations\.

##### Hyperparameter optimization\.

All positive\-constrained hyperparameters \(lengthscales, variances\) are optimized in an unconstrained space via a Softplus transformation\. Specifically, for a positive hyperparameterθ\>0\\theta\>0, we introduce an unconstrained parameterθ~∈ℝ\\tilde\{\\theta\}\\in\\mathbb\{R\}and defineθ=log⁡\(1\+exp⁡\(θ~\)\)=softplus​\(θ~\)\\theta=\\log\(1\+\\exp\(\\tilde\{\\theta\}\)\)=\\text\{softplus\}\(\\tilde\{\\theta\}\)\. The gradient with respect toθ~\\tilde\{\\theta\}is computed using the chain rule:∂ℒ∂θ~=∂ℒ∂θ⋅σ⁡\(θ~\)\\frac\{\\partial\\mathcal\{L\}\}\{\\partial\\tilde\{\\theta\}\}=\\frac\{\\partial\\mathcal\{L\}\}\{\\partial\\theta\}\\cdot\\sigma\(\\tilde\{\\theta\}\), whereσ⁡\(⋅\)\\sigma\(\\cdot\)is the sigmoid function\. This transformation ensures that the hyperparameters remain positive throughout optimization while allowing unconstrained gradient\-based optimization algorithms \(e\.g\., Adam, SGD\) to be applied directly\.

##### Inducing locations and kernel hyperparameters\.

We jointly optimize inducing point locations and kernel parameters via gradient ascent on the ELBO \(marginal likelihood\)\. The global component and the client\-specific deviations share the same inducing locations𝐙g\\mathbf\{Z\}\_\{g\}to preserve a common basis, while each local component uses its own𝐙i\\mathbf\{Z\}\_\{i\}to capture client\-specific residual structure\. For the global kernelkgk\_\{g\}with hyperparameters\(ℓg,σg2\)\(\\ell\_\{g\},\\sigma\_\{g\}^\{2\}\)and the local kernelkik\_\{i\}with hyperparameters\(ℓi,σi2\)\(\\ell\_\{i\},\\sigma\_\{i\}^\{2\}\), we optimize\{𝐙g,ℓg,σg2,ηϕ,\{𝐙i,ℓi,σi2\}i=1T\}\\\{\\mathbf\{Z\}\_\{g\},\\ell\_\{g\},\\sigma\_\{g\}^\{2\},\\eta\_\{\\phi\},\\\{\\mathbf\{Z\}\_\{i\},\\ell\_\{i\},\\sigma\_\{i\}^\{2\}\\\}\_\{i=1\}^\{T\}\\\}jointly with the variational parameters by maximizing the ELBO, withϕ=softplus⁡\(ηϕ\)\\phi=\\operatorname\{softplus\}\(\\eta\_\{\\phi\}\)enforced by the transformation above\. The gradients∇𝐙gℒ\\nabla\_\{\\mathbf\{Z\}\_\{g\}\}\\mathcal\{L\},∇ℓgℒ\\nabla\_\{\\ell\_\{g\}\}\\mathcal\{L\},∇σg2ℒ\\nabla\_\{\\sigma\_\{g\}^\{2\}\}\\mathcal\{L\},∇𝐙iℒ\\nabla\_\{\\mathbf\{Z\}\_\{i\}\}\\mathcal\{L\},∇ℓiℒ\\nabla\_\{\\ell\_\{i\}\}\\mathcal\{L\}, and∇σi2ℒ\\nabla\_\{\\sigma\_\{i\}^\{2\}\}\\mathcal\{L\}are computed via automatic differentiation, enabling automatic adaptation of both global and local kernels\.

##### VFE\-style diagonal computations\.

For efficiency we employ diagonal computations for the trace terms in the VFE objective\. Specifically, the trace regularization termsTr⁡\(𝐊g,i​i−𝐐g,i​i\)\\operatorname\{Tr\}\(\\mathbf\{K\}\_\{g,ii\}\-\\mathbf\{Q\}\_\{g,ii\}\),Tr⁡\(𝐊δ,i​i−𝐐δ,i​i\)\\operatorname\{Tr\}\(\\mathbf\{K\}\_\{\\delta,ii\}\-\\mathbf\{Q\}\_\{\\delta,ii\}\), andTr⁡\(𝐊i,i​i−𝐐i,i​i\)\\operatorname\{Tr\}\(\\mathbf\{K\}\_\{i,ii\}\-\\mathbf\{Q\}\_\{i,ii\}\)are computed by summing only the diagonal elements:∑j=1ni\[𝐊⋅,i​i−𝐐⋅,i​i\]j​j\\sum\_\{j=1\}^\{n\_\{i\}\}\[\\mathbf\{K\}\_\{\\cdot,ii\}\-\\mathbf\{Q\}\_\{\\cdot,ii\}\]\_\{jj\}, where\[⋅\]j​j\[\\cdot\]\_\{jj\}denotes thejj\-th diagonal element\. Our inference is framed in the VFE objective: we retain the exact GP prior and rely on the trace of the residual covariance \(which depends only on diagonal elements\) rather than assuming a diagonal observation covariance structure\. This approach provides computational efficiency while maintaining the exact prior and avoiding restrictive assumptions on the covariance structure\.

### 3\.8Practical Usage of pFedHGP

Beyond regression and calibration tasks, pFedHGP naturally extends to classification and clustering applications by leveraging its hierarchical predictive distributions\. Section[3\.8\.1](https://arxiv.org/html/2609.19337#S3.SS8.SSS1)presents a maximum\-evidence classification approach that uses the full Gaussian predictive log\-likelihood under \([30](https://arxiv.org/html/2609.19337#S3.E30)\) to assign test samples to client models, enabling federated fault diagnosis and anomaly detection while preserving privacy\. Section[3\.8\.2](https://arxiv.org/html/2609.19337#S3.SS8.SSS2)introduces a structural comparison framework that computes client\-specific projection operators in the inducing feature space, enabling unsupervised discovery of client groups with similar patterns without sharing raw observations\.

#### 3\.8\.1Classification via pFedHGP

For each clientii, the variational posterior induces a Gaussian predictive law for observations at arbitrary inputs\. At the training inputs𝐗i\\mathbf\{X\}\_\{i\}, we have already expressed the expected log likelihood in closed form in \([16](https://arxiv.org/html/2609.19337#S3.E16)\) in terms of the hierarchical predictive mean𝝁𝐟i\\bm\{\\mu\}\_\{\\mathbf\{f\}\_\{i\}\}, posterior covariance in each inducing block \(𝐒g\\mathbf\{S\}\_\{g\},𝐒δi\\mathbf\{S\}\_\{\\delta\_\{i\}\},𝐒ui\\mathbf\{S\}\_\{u\_\{i\}\}\), and the residual covariance𝐑i\\mathbf\{R\}\_\{i\}\. The mean decomposes into global, client\-specific deviation, and local contributions as in \([17](https://arxiv.org/html/2609.19337#S3.E17)\), and the residual traceTr⁡\(𝐑i\)\\operatorname\{Tr\}\(\\mathbf\{R\}\_\{i\}\)decomposes additively across the same three components as in \([10](https://arxiv.org/html/2609.19337#S3.E10)\)\.

For a batch of test inputs𝐗⋆\\mathbf\{X\}\_\{\\star\}and observations𝐲⋆\\mathbf\{y\}\_\{\\star\}, we reuse exactly this structure\. Let𝐀g,i⋆\\mathbf\{A\}\_\{g,i\}^\{\\star\},𝐀δ,i⋆\\mathbf\{A\}\_\{\\delta,i\}^\{\\star\}, and𝐀i⋆\\mathbf\{A\}\_\{i\}^\{\\star\}denote the interpolation matrices from Proposition[2](https://arxiv.org/html/2609.19337#Thmproposition2)evaluated at𝐗⋆\\mathbf\{X\}\_\{\\star\}\(so𝐀g,i⋆=𝐊g​\(𝐗⋆,𝐙g\)​𝐊g​g−1\\mathbf\{A\}\_\{g,i\}^\{\\star\}=\\mathbf\{K\}\_\{g\}\(\\mathbf\{X\}\_\{\\star\},\\mathbf\{Z\}\_\{g\}\)\\mathbf\{K\}\_\{gg\}^\{\-1\},𝐀δ,i⋆=𝐀g,i⋆\\mathbf\{A\}\_\{\\delta,i\}^\{\\star\}=\\mathbf\{A\}\_\{g,i\}^\{\\star\}sincekδ=ϕ​kgk\_\{\\delta\}=\\phi k\_\{g\}, and𝐀i⋆=𝐊i​\(𝐗⋆,𝐙i\)​𝐊i​i−1\\mathbf\{A\}\_\{i\}^\{\\star\}=\\mathbf\{K\}\_\{i\}\(\\mathbf\{X\}\_\{\\star\},\\mathbf\{Z\}\_\{i\}\)\\mathbf\{K\}\_\{ii\}^\{\-1\}\)\. Under the trained variational posterior, clientiiassigns the Gaussian predictive law

𝐲⋆\|𝐗⋆,i∼𝒩⁡\(𝝁𝐟i​\(𝐗⋆\),𝚺i​\(𝐗⋆\)\),\\mathbf\{y\}\_\{\\star\}\\mid\\mathbf\{X\}\_\{\\star\},i~\\sim~\\mathcal\{N\}\\\!\\bigl\(\\bm\{\\mu\}\_\{\\mathbf\{f\}\_\{i\}\}\(\\mathbf\{X\}\_\{\\star\}\),\\,\\bm\{\\Sigma\}\_\{i\}\(\\mathbf\{X\}\_\{\\star\}\)\\bigr\),with hierarchical predictive mean𝝁𝐟i​\(𝐗⋆\)\\bm\{\\mu\}\_\{\\mathbf\{f\}\_\{i\}\}\(\\mathbf\{X\}\_\{\\star\}\)as in \([17](https://arxiv.org/html/2609.19337#S3.E17)\) and predictive covariance

𝚺i​\(𝐗⋆\)=𝐀g,i⋆​𝐒g​\(𝐀g,i⋆\)⊤\+𝐀δ,i⋆​𝐒δi​\(𝐀δ,i⋆\)⊤\+𝐀i⋆​𝐒ui​\(𝐀i⋆\)⊤\+𝐑i​\(𝐗⋆\)\+σ2​𝐈,\\bm\{\\Sigma\}\_\{i\}\(\\mathbf\{X\}\_\{\\star\}\)=\\mathbf\{A\}\_\{g,i\}^\{\\star\}\\mathbf\{S\}\_\{g\}\(\\mathbf\{A\}\_\{g,i\}^\{\\star\}\)^\{\\top\}\+\\mathbf\{A\}\_\{\\delta,i\}^\{\\star\}\\mathbf\{S\}\_\{\\delta\_\{i\}\}\(\\mathbf\{A\}\_\{\\delta,i\}^\{\\star\}\)^\{\\top\}\+\\mathbf\{A\}\_\{i\}^\{\\star\}\\mathbf\{S\}\_\{u\_\{i\}\}\(\\mathbf\{A\}\_\{i\}^\{\\star\}\)^\{\\top\}\+\\mathbf\{R\}\_\{i\}\(\\mathbf\{X\}\_\{\\star\}\)\+\\sigma^\{2\}\\mathbf\{I\},\(30\)where the first three terms are inducing\-posterior covariance contributions and𝐑i​\(𝐗⋆\)\\mathbf\{R\}\_\{i\}\(\\mathbf\{X\}\_\{\\star\}\)is the VFE residual covariance at𝐗⋆\\mathbf\{X\}\_\{\\star\}not explained by the inducing variables \(with trace decomposition as in \([10](https://arxiv.org/html/2609.19337#S3.E10)\)\); the sum is the full predictive covariance and does not double\-count uncertainty\. The blocks in \([30](https://arxiv.org/html/2609.19337#S3.E30)\) are the same uncertainty components that enter the training\-time VFE objective \([16](https://arxiv.org/html/2609.19337#S3.E16)\), but classification uses the*marginal*Gaussian log\-density at test inputs, not the trace\-weighted expected log\-likelihood from the ELBO\. In particular, \([16](https://arxiv.org/html/2609.19337#S3.E16)\) scales squared error and all covariance contributions by1/\(2​σ2\)1/\(2\\sigma^\{2\}\)and does not include the𝚺i−1\\bm\{\\Sigma\}\_\{i\}^\{\-1\}quadratic form orlog⁡\|𝚺i\|\\log\|\\bm\{\\Sigma\}\_\{i\}\|term of a full multivariate predictive law\.

Letn⋆:=dim\(𝐲⋆\)n\_\{\\star\}:=\\dim\(\\mathbf\{y\}\_\{\\star\}\)\. The classification score is the Gaussian predictive log\-likelihood

ℓi​\(𝐗⋆,𝐲⋆\)\\displaystyle\\ell\_\{i\}\(\\mathbf\{X\}\_\{\\star\},\\mathbf\{y\}\_\{\\star\}\):=log⁡𝒩⁡\(𝐲⋆,𝝁𝐟i​\(𝐗⋆\),𝚺i​\(𝐗⋆\)\)\\displaystyle:=\\log\\mathcal\{N\}\\\!\\bigl\(\\mathbf\{y\}\_\{\\star\};\\,\\bm\{\\mu\}\_\{\\mathbf\{f\}\_\{i\}\}\(\\mathbf\{X\}\_\{\\star\}\),\\,\\bm\{\\Sigma\}\_\{i\}\(\\mathbf\{X\}\_\{\\star\}\)\\bigr\)\(31\)=−12​\[\(𝐲⋆−𝝁𝐟i​\(𝐗⋆\)\)⊤​𝚺i​\(𝐗⋆\)−1​\(𝐲⋆−𝝁𝐟i​\(𝐗⋆\)\)\+log⁡\|𝚺i​\(𝐗⋆\)\|\+n⋆​log⁡\(2​π\)\]\.\\displaystyle=\-\\frac\{1\}\{2\}\\Bigl\[\\bigl\(\\mathbf\{y\}\_\{\\star\}\-\\bm\{\\mu\}\_\{\\mathbf\{f\}\_\{i\}\}\(\\mathbf\{X\}\_\{\\star\}\)\\bigr\)^\{\\top\}\\bm\{\\Sigma\}\_\{i\}\(\\mathbf\{X\}\_\{\\star\}\)^\{\-1\}\\bigl\(\\mathbf\{y\}\_\{\\star\}\-\\bm\{\\mu\}\_\{\\mathbf\{f\}\_\{i\}\}\(\\mathbf\{X\}\_\{\\star\}\)\\bigr\)\+\\log\\bigl\|\\bm\{\\Sigma\}\_\{i\}\(\\mathbf\{X\}\_\{\\star\}\)\\bigr\|\+n\_\{\\star\}\\log\(2\\pi\)\\Bigr\]\.We evaluate \([31](https://arxiv.org/html/2609.19337#S3.E31)\) in practice via Cholesky factorization of𝚺i​\(𝐗⋆\)\\bm\{\\Sigma\}\_\{i\}\(\\mathbf\{X\}\_\{\\star\}\)\.

We then use a standard model\-evidence decision rule, treating each client model as a class\-conditional generator\. The predicted label for\(𝐗⋆,𝐲⋆\)\(\\mathbf\{X\}\_\{\\star\},\\mathbf\{y\}\_\{\\star\}\)is chosen by maximum predictive log likelihood:

c^=arg⁡maxi∈\{1,…,T\}​ℓi​\(𝐗⋆,𝐲⋆\)\.\\hat\{c\}\\;=\\;\\arg\\max\_\{i\\in\\\{1,\\dots,T\\\}\}~\\ell\_\{i\}\(\\mathbf\{X\}\_\{\\star\},\\mathbf\{y\}\_\{\\star\}\)\.\(32\)Because𝝁𝐟i​\(𝐗⋆\)\\bm\{\\mu\}\_\{\\mathbf\{f\}\_\{i\}\}\(\\mathbf\{X\}\_\{\\star\}\),𝚺i​\(𝐗⋆\)\\bm\{\\Sigma\}\_\{i\}\(\\mathbf\{X\}\_\{\\star\}\), and𝐑i​\(𝐗⋆\)\\mathbf\{R\}\_\{i\}\(\\mathbf\{X\}\_\{\\star\}\)inherit the hierarchical decompositions in \([17](https://arxiv.org/html/2609.19337#S3.E17)\) and \([10](https://arxiv.org/html/2609.19337#S3.E10)\), the scoreℓi\\ell\_\{i\}automatically discounts patterns shared through the global component while emphasizing client\-specific offsets and local structure, which is particularly valuable in heterogeneous classification settings\.

#### 3\.8\.2Clustering via pFedHGP

For structural comparison across clients, each client is represented by a projection operator in the inducing feature space that captures the subspaces effectively spanned by its offset and local components\. Using only𝐗i\\mathbf\{X\}\_\{i\}and the trained kernels \(privacy preserved\), define the offset operatorAδ,i=Kg​\(𝐙g,𝐗i\)​\(Kg​\(𝐗i,𝐗i\)\+σ2​Ini\)−1​Kg​\(𝐗i,𝐙g\)∈ℝM×MA\_\{\\delta,i\}=K\_\{g\}\(\\mathbf\{Z\}\_\{g\},\\mathbf\{X\}\_\{i\}\)\\,\(K\_\{g\}\(\\mathbf\{X\}\_\{i\},\\mathbf\{X\}\_\{i\}\)\+\\sigma^\{2\}I\_\{n\_\{i\}\}\)^\{\-1\}\\,K\_\{g\}\(\\mathbf\{X\}\_\{i\},\\mathbf\{Z\}\_\{g\}\)\\in\\mathbb\{R\}^\{M\\times M\}, which leverages the global kernelKgK\_\{g\}since the client\-specific offsetδi\\delta\_\{i\}does not possess an independent kernel\. The local operator isAloc,i=Ki​\(𝐙i,𝐗i\)​\(Ki​\(𝐗i,𝐗i\)\+σ2​Ini\)−1​Ki​\(𝐗i,𝐙i\)∈ℝMi×MiA\_\{\\mathrm\{loc\},i\}=K\_\{i\}\(\\mathbf\{Z\}\_\{i\},\\mathbf\{X\}\_\{i\}\)\\,\(K\_\{i\}\(\\mathbf\{X\}\_\{i\},\\mathbf\{X\}\_\{i\}\)\+\\sigma^\{2\}I\_\{n\_\{i\}\}\)^\{\-1\}\\,K\_\{i\}\(\\mathbf\{X\}\_\{i\},\\mathbf\{Z\}\_\{i\}\)\\in\\mathbb\{R\}^\{M\_\{i\}\\times M\_\{i\}\}\. The joint personalized operator is given byAi=blockdiag⁡\(Aδ,i,Aloc,i\)A\_\{i\}=\\mathrm\{blockdiag\}\(A\_\{\\delta,i\},A\_\{\\mathrm\{loc\},i\}\), and its trace\-normalized version isA^i=Ai/tr⁡\(Ai\)\\widehat\{A\}\_\{i\}=A\_\{i\}/\\mathrm\{tr\}\(A\_\{i\}\)\. Dissimilarity between clients is quantified by the Frobenius distanceρ⁡\(i,j\)=‖A^i−A^j‖F2\\rho\(i,j\)=\\\|\\widehat\{A\}\_\{i\}\-\\widehat\{A\}\_\{j\}\\\|\_\{F\}^\{2\}, and the dissimilarity matrix\[ρ⁡\(i,j\)\]i,j=1T\\big\[\\rho\(i,j\)\\big\]\_\{i,j=1\}^\{T\}can be fed directly into standard clustering routines\. SinceA^i\\widehat\{A\}\_\{i\}is computed locally from𝐗i\\mathbf\{X\}\_\{i\}and kernel evaluations, clustering requires no exposure of raw observations and remains coherent with the federated training pipeline\.

### 3\.9ILMM Formulation of pFedHGP

We connect pFedHGP with the*Instantaneous Linear Mixing Model*\(ILMM\)\([Bruinsma et al\., 2020](https://arxiv.org/html/2609.19337#bib.bib2)\), a framework for constructing multi\-output Gaussian processes by linearly mixing latent Gaussian processes\. In function space, the hierarchical decomposition from Assumption[1](https://arxiv.org/html/2609.19337#Thmassumption1)\(wherefitot​\(x\)=fg​\(x\)\+fδ,i​\(x\)\+fi​\(x\)f\_\{i\}^\{\\mathrm\{tot\}\}\(x\)=f\_\{g\}\(x\)\+f\_\{\\delta,i\}\(x\)\+f\_\{i\}\(x\)\) can be recast as an ILMM by stacking all clients’ latent functions into aTT\-dimensional output vector𝐟tot​\(x\)\\mathbf\{f\}^\{\\mathrm\{tot\}\}\(x\)and a\(2​T\+1\)\(2T\+1\)\-dimensional latent vector𝐱⁡\(x\)=\[fg​\(x\),fδ,1​\(x\),…,fδ,T​\(x\),f1​\(x\),…,fT​\(x\)\]⊤\\mathbf\{x\}\(x\)=\[f\_\{g\}\(x\),f\_\{\\delta,1\}\(x\),\\dots,f\_\{\\delta,T\}\(x\),f\_\{1\}\(x\),\\dots,f\_\{T\}\(x\)\]^\{\\top\}\. The mapping from latent to output processes is then linear:

𝐟tot​\(x\)=H​𝐱​\(x\),H=\[1T​\|IT\|​IT\]∈ℝT×\(2​T\+1\),\\mathbf\{f\}^\{\\mathrm\{tot\}\}\(x\)=H\\,\\mathbf\{x\}\(x\),\\quad H=\\bigl\[\\,\\mathbf\{1\}\_\{T\}\\;\\big\|\\;I\_\{T\}\\;\\big\|\\;I\_\{T\}\\,\\bigr\]\\in\\mathbb\{R\}^\{T\\times\(2T\+1\)\},\(33\)where𝟏T\\mathbf\{1\}\_\{T\}is theTT\-vector of ones andITI\_\{T\}is theT×TT\\times Tidentity matrix\. With observation noiseΣ=σ2​IT\\Sigma=\\sigma^\{2\}I\_\{T\}, this yields an ILMM withp=Tp=Toutputs andm=2​T\+1m=2T\+1latent processes\. The choicekδ=ϕ​kgk\_\{\\delta\}=\\phi k\_\{g\}\(from Assumption[1](https://arxiv.org/html/2609.19337#Thmassumption1)\) ensures that client\-specific deviationsfδ,if\_\{\\delta,i\}live in the same latent subspace as the shared global processfgf\_\{g\}\.

In inducing\-variable space, Proposition[2](https://arxiv.org/html/2609.19337#Thmproposition2)reveals a similar linear mixing structure through the low\-rank expansion \([8](https://arxiv.org/html/2609.19337#S3.E8)\), where the shared basis𝐀g,i\\mathbf\{A\}\_\{g,i\}\(common to𝐮g\\mathbf\{u\}\_\{g\}and𝜹i\\bm\{\\delta\}\_\{i\}\) and client\-specific local bases𝐀i\\mathbf\{A\}\_\{i\}organize the hierarchical decomposition\. Under our VFE implementation, the residual covariance𝐑i\\mathbf\{R\}\_\{i\}\(Proposition[3](https://arxiv.org/html/2609.19337#Thmproposition3)\) is incorporated through trace regularization terms in the ELBO \(Theorem[1](https://arxiv.org/html/2609.19337#Thmtheorem1)\), rather than directly modifying the observation model\. This ILMM perspective embeds pFedHGP within the multi\-output GP mixing\-model hierarchy while maintaining the privacy\-preserving federated structure\.

##### Scalar\-to\-multi\-sensor vector extension\.

The same ILMM viewpoint extends pFedHGP from scalar client responses to synchronized multi\-sensor observations\. Suppose clientiiobserves𝐲i​\(x\)=\(yi​1​\(x\),…,yi​Q​\(x\)\)⊤∈ℝQ\\mathbf\{y\}\_\{i\}\(x\)=\(y\_\{i1\}\(x\),\\ldots,y\_\{iQ\}\(x\)\)^\{\\top\}\\in\\mathbb\{R\}^\{Q\}\. A vector\-valued pFedHGP preserves the same hierarchical decomposition by representing theQQ\-dimensional global, deviation, and local components through latent GP factors\. We denote ILMM*loading*matrices byB⋅B\_\{\\cdot\}to distinguish them from sparse\-GP*interpolation*matrices𝐀⋅,i\\mathbf\{A\}\_\{\\cdot,i\}in Section[3\.3](https://arxiv.org/html/2609.19337#S3.SS3):

𝐲i​\(x\)=Bg​𝐡g​\(x\)\+Bδ,i​𝐡δ,i​\(x\)\+Bl,i​𝐡l,i​\(x\)\+ϵi​\(x\),ϵi​\(x\)∼𝒩⁡\(𝟎,Σϵ\),\\mathbf\{y\}\_\{i\}\(x\)=B\_\{g\}\\mathbf\{h\}\_\{g\}\(x\)\+B\_\{\\delta,i\}\\mathbf\{h\}\_\{\\delta,i\}\(x\)\+B\_\{l,i\}\\mathbf\{h\}\_\{l,i\}\(x\)\+\\bm\{\\epsilon\}\_\{i\}\(x\),\\qquad\\bm\{\\epsilon\}\_\{i\}\(x\)\\sim\\mathcal\{N\}\(\\mathbf\{0\},\\Sigma\_\{\\epsilon\}\),\(34\)whereBg∈ℝQ×RgB\_\{g\}\\in\\mathbb\{R\}^\{Q\\times R\_\{g\}\},Bδ,i∈ℝQ×RδB\_\{\\delta,i\}\\in\\mathbb\{R\}^\{Q\\times R\_\{\\delta\}\}, andBl,i∈ℝQ×RlB\_\{l,i\}\\in\\mathbb\{R\}^\{Q\\times R\_\{l\}\}map latent global, deviation, and local GP factors to theQQsensor channels\. The latent vectors𝐡g​\(x\)\\mathbf\{h\}\_\{g\}\(x\),𝐡δ,i​\(x\)\\mathbf\{h\}\_\{\\delta,i\}\(x\), and𝐡l,i​\(x\)\\mathbf\{h\}\_\{l,i\}\(x\)represent shared population\-level factors, client\-specific systematic calibration factors, and local residual factors, respectively\.

This construction induces the cross\-client and cross\-sensor covariance

Cov⁡\{𝐲i​\(x\),𝐲j​\(x′\)\}\\displaystyle\\operatorname\{Cov\}\\\{\\mathbf\{y\}\_\{i\}\(x\),\\mathbf\{y\}\_\{j\}\(x^\{\\prime\}\)\\\}=Bg​Kg​\(x,x′\)​Bg⊤\\displaystyle=B\_\{g\}K\_\{g\}\(x,x^\{\\prime\}\)B\_\{g\}^\{\\top\}\(35\)\+𝟏\[i=j\]Bδ,iKδ\(x,x′\)Bδ,i⊤\\displaystyle\+\\mathbf\{1\}\[i=j\]B\_\{\\delta,i\}K\_\{\\delta\}\(x,x^\{\\prime\}\)B\_\{\\delta,i\}^\{\\top\}\+𝟏\[i=j\]Bl,iKl\(x,x′\)Bl,i⊤\+𝟏\[i=j\]Σϵ𝟏\[x=x′\],\\displaystyle\+\\mathbf\{1\}\[i=j\]B\_\{l,i\}K\_\{l\}\(x,x^\{\\prime\}\)B\_\{l,i\}^\{\\top\}\+\\mathbf\{1\}\[i=j\]\\Sigma\_\{\\epsilon\}\\mathbf\{1\}\[x=x^\{\\prime\}\],whereKg​\(x,x′\)=Cov⁡\{𝐡g​\(x\),𝐡g​\(x′\)\}K\_\{g\}\(x,x^\{\\prime\}\)=\\operatorname\{Cov\}\\\{\\mathbf\{h\}\_\{g\}\(x\),\\mathbf\{h\}\_\{g\}\(x^\{\\prime\}\)\\\},Kδ​\(x,x′\)=Cov⁡\{𝐡δ,i​\(x\),𝐡δ,i​\(x′\)\}K\_\{\\delta\}\(x,x^\{\\prime\}\)=\\operatorname\{Cov\}\\\{\\mathbf\{h\}\_\{\\delta,i\}\(x\),\\mathbf\{h\}\_\{\\delta,i\}\(x^\{\\prime\}\)\\\}, andKl​\(x,x′\)=Cov⁡\{𝐡l,i​\(x\),𝐡l,i​\(x′\)\}K\_\{l\}\(x,x^\{\\prime\}\)=\\operatorname\{Cov\}\\\{\\mathbf\{h\}\_\{l,i\}\(x\),\\mathbf\{h\}\_\{l,i\}\(x^\{\\prime\}\)\\\}\. Equation \([35](https://arxiv.org/html/2609.19337#S3.E35)\) shows that cross\-sensor dependence is generated by latent trajectories projected through the loading matrices\. The global term is shared across clients, whereas the deviation, local, and noise terms are client\-specific\.

##### Summary and intuition \(ILMM\)\.

The ILMM viewpoint supplies a positive\-semidefiniteQ×QQ\\times Qcovariance construction for correlated sensor channels; it does not introduce a separate scalar computational shortcut\. Scalability follows from sparse variational inducing variables, trace regularization, and federated aggregation of global statistics\.

## 4Numerical Study

This section evaluates pFedHGP in a controlled multi\-output federated regression problem\. The experiment is designed to test four properties that are central to the proposed framework: predictive accuracy, recovery of cross\-output dependence, uncertainty calibration, and robustness to violations of the assumed additive hierarchy\. Each client observes multiple correlated output channels over a common input domain\. The ground\-truth data\-generating process is known, which permits direct evaluation of both prediction quality and the learned output covariance structure\. All reported multi\-output results are based on30 independent runs\. In each run we regenerate the latent component trajectories, draw new component loading matrices, resample client training inputs, reinitialize every fitted model, and draw fresh observation noise\. Metrics are averaged over theT=6T=6clients within each run; tables report run\-level means with95%95\\%confidence intervals computed across the 30 runs\.

### 4\.1Multi\-output Synthetic Benchmark

We considerT=6T=6clients andQ=4Q=4output channels over the input domain\[0,10\]\[0,10\]\. For each independent run, clienti∈\{1,…,T\}i\\in\\\{1,\\ldots,T\\\}receivesntrain=50n\_\{\\mathrm\{train\}\}=50training inputs sampled uniformly from\[0,10\]\[0,10\], and all clients are evaluated on a common dense test grid withntest=150n\_\{\\mathrm\{test\}\}=150inputs\. The observation at inputxxis a vector𝐲i​\(x\)=\(yi​1​\(x\),…,yi​Q​\(x\)\)⊤∈ℝQ\\mathbf\{y\}\_\{i\}\(x\)=\(y\_\{i1\}\(x\),\\ldots,y\_\{iQ\}\(x\)\)^\{\\top\}\\in\\mathbb\{R\}^\{Q\}, generated by the additive multi\-output hierarchy

𝐲i​\(x\)=𝐟g​\(x\)\+𝐟δ,i​\(x\)\+𝐟i​\(x\)⏟𝐲iclean​\(x\)\+ϵi​\(x\),ϵi​\(x\)∼N⁡\(𝟎,σ2​IQ\),\\mathbf\{y\}\_\{i\}\(x\)=\\underbrace\{\\mathbf\{f\}\_\{g\}\(x\)\+\\mathbf\{f\}\_\{\\delta,i\}\(x\)\+\\mathbf\{f\}\_\{i\}\(x\)\}\_\{\\mathbf\{y\}^\{\\mathrm\{clean\}\}\_\{i\}\(x\)\}\+\\bm\{\\epsilon\}\_\{i\}\(x\),\\qquad\\bm\{\\epsilon\}\_\{i\}\(x\)\\sim N\(\\mathbf\{0\},\\sigma^\{2\}I\_\{Q\}\),\(36\)whereσ=0\.05\\sigma=0\.05\. The components𝐟g\\mathbf\{f\}\_\{g\},𝐟δ,i\\mathbf\{f\}\_\{\\delta,i\}, and𝐟i\\mathbf\{f\}\_\{i\}represent, respectively, the shared global function, the systematic client\-specific deviation, and the local site\-specific variation\. Cross\-output dependence is introduced within these vector\-valued component processes\. All component processes are generated from RBF covariance functions of the form

k⁡\(x,x′\)=τ2​exp⁡\{−\(x−x′\)22​ℓ2\}\.k\(x,x^\{\\prime\}\)=\\tau^\{2\}\\exp\\left\\\{\-\\frac\{\(x\-x^\{\\prime\}\)^\{2\}\}\{2\\ell^\{2\}\}\\right\\\}\.\(37\)To induce correlated four\-dimensional outputs, each vector\-valued component is formed by mixingR=2R=2independent scalar GP draws into theQ=4Q=4output channels\. Concretely, for componentc∈\{g,δ,l\}c\\in\\\{g,\\delta,l\\\}we draw𝐡c​\(x\)=\(hc,1​\(x\),hc,2​\(x\)\)⊤\\mathbf\{h\}\_\{c\}\(x\)=\(h\_\{c,1\}\(x\),h\_\{c,2\}\(x\)\)^\{\\top\}and set

𝐟c​\(x\)=Bc​𝐡c​\(x\),Bc∈ℝQ×R,\\mathbf\{f\}\_\{c\}\(x\)=B\_\{c\}\\,\\mathbf\{h\}\_\{c\}\(x\),\\qquad B\_\{c\}\\in\\mathbb\{R\}^\{Q\\times R\},with independent standard\-normal entries in eachBcB\_\{c\}resampled at every run \(and client\-specificBδ,iB\_\{\\delta,i\},Bl,iB\_\{l,i\}drawn independently in the deviation and local components\)\. The two global latent draws use lengthscalesℓg,1=1\.8\\ell\_\{g,1\}=1\.8andℓg,2=2\.5\\ell\_\{g,2\}=2\.5, with output scalesτg,1=1\.0\\tau\_\{g,1\}=1\.0andτg,2=0\.7\\tau\_\{g,2\}=0\.7\. The deviation latents use the same smooth RBF lengthscales as the global component, with simulator deviation scaleϕsim=1\.20\\phi\_\{\\mathrm\{sim\}\}=1\.20\. The local latents use shorter lengthscalesℓl,1=0\.25\\ell\_\{l,1\}=0\.25andℓl,2=0\.45\\ell\_\{l,2\}=0\.45, with output scalesτl,1=1\.0\\tau\_\{l,1\}=1\.0andτl,2=0\.8\\tau\_\{l,2\}=0\.8\.

### 4\.2Fitted ILMM–pFedHGP configuration and training

All fitted methods target the ILMM observation model in \([34](https://arxiv.org/html/2609.19337#S3.E34)\) withQ=4Q=4and latent ranksRg=Rδ=Rl=2R\_\{g\}=R\_\{\\delta\}=R\_\{l\}=2, matching the two\-latent\-draw simulator above\. For each scalar latent dimension we use the sparse hierarchical pFedHGP parameterization from Sections[3\.3](https://arxiv.org/html/2609.19337#S3.SS3)–[3\.6](https://arxiv.org/html/2609.19337#S3.SS6): shared global inducing locations𝐙g∈ℝM×1\\mathbf\{Z\}\_\{g\}\\in\\mathbb\{R\}^\{M\\times 1\}withM=25M=25points initialized on a uniform grid over\[0,10\]\[0,10\], client\-specific local inducing grids𝐙i∈ℝMi×1\\mathbf\{Z\}\_\{i\}\\in\\mathbb\{R\}^\{M\_\{i\}\\times 1\}withMi=25M\_\{i\}=25, squared\-exponential kernelskgk\_\{g\}and\{ki\}\\\{k\_\{i\}\\\}, and deviation layerkδ=ϕ​kgk\_\{\\delta\}=\\phi k\_\{g\}with population scaleϕ\\philearned on the server\. The loading matrices are*not*fixed to the simulator draws: we treatBg∈ℝ4×2B\_\{g\}\\in\\mathbb\{R\}^\{4\\times 2\}as a globally shared learnable block and\{Bδ,i,Bl,i\}\\\{B\_\{\\delta,i\},B\_\{l,i\}\\\}as client\-specific learnable blocks, all optimized jointly with the variational parameters and kernel hyperparameters by gradient ascent on the federated ELBO\. Observation noise variance is learned through a softplus reparameterization with initializationσ2=σsim2=0\.052\\sigma^\{2\}=\\sigma\_\{\\mathrm\{sim\}\}^\{2\}=0\.05^\{2\}\. Initial kernel lengthscales are set toℓg=2\.0\\ell\_\{g\}=2\.0,ℓi=0\.35\\ell\_\{i\}=0\.35, and output variances to1\.01\.0; inducing locations remain learnable as in Section[3\.7](https://arxiv.org/html/2609.19337#S3.SS7)\.

Training follows Algorithm[1](https://arxiv.org/html/2609.19337#alg1)forRcomm=20R\_\{\\mathrm\{comm\}\}=20synchronous communication rounds\. In roundrr, each client first runsSloc=80S\_\{\\mathrm\{loc\}\}=80Adam steps \(learning rateηloc=0\.1\\eta\_\{\\mathrm\{loc\}\}=0\.1\) on its local ELBO block withΘglobal\(r\)\\Theta\_\{\\mathrm\{global\}\}^\{\(r\)\}held fixed, then transmits𝒮~i\(r\+1\)\\widetilde\{\\mathcal\{S\}\}\_\{i\}^\{\(r\+1\)\}from \([23](https://arxiv.org/html/2609.19337#S3.E23)\)\. The server forms∇Θglobalℒ\\nabla\_\{\\Theta\_\{\\mathrm\{global\}\}\}\\mathcal\{L\}via \([24](https://arxiv.org/html/2609.19337#S3.E24)\) \(theηϕ\\eta\_\{\\phi\}block is the corresponding component of∑i𝒮~i\(r\+1\)\\sum\_\{i\}\\widetilde\{\\mathcal\{S\}\}\_\{i\}^\{\(r\+1\)\}, cf\. \([26](https://arxiv.org/html/2609.19337#S3.E26)\)\) and applies one Adam step with learning rateηglobal=0\.1\\eta\_\{\\mathrm\{global\}\}=0\.1\. Reduced models use the same protocol:No\-deviationremoves𝐟δ,i\\mathbf\{f\}\_\{\\delta,i\}andBδ,iB\_\{\\delta,i\};No\-localremoves𝐟i\\mathbf\{f\}\_\{i\}andBl,iB\_\{l,i\};Global\-onlyandLocal\-onlyretain only the corresponding hierarchical blocks;pFedGPkeeps client\-specific multi\-output blocks with FedAvg\-style synchronization of shared input\-kernel hyperparameters but without the explicit global/deviation/local ILMM decomposition\.

At test time, clientiireturns the joint predictive Gaussian𝒩⁡\(𝝁^i​\(x\),Σ^i​\(x\)\)\\mathcal\{N\}\(\\widehat\{\\bm\{\\mu\}\}\_\{i\}\(x\),\\widehat\{\\Sigma\}\_\{i\}\(x\)\)obtained by composing the hierarchical predictive mean and covariance in Section[3\.8\.1](https://arxiv.org/html/2609.19337#S3.SS8.SSS1)with the learned ILMM loadings\. Let𝐀g⋆​\(x\)∈ℝRg×M\\mathbf\{A\}\_\{g\}^\{\\star\}\(x\)\\in\\mathbb\{R\}^\{R\_\{g\}\\times M\},𝐀δ,i⋆​\(x\)∈ℝRδ×M\\mathbf\{A\}\_\{\\delta,i\}^\{\\star\}\(x\)\\in\\mathbb\{R\}^\{R\_\{\\delta\}\\times M\}, and𝐀l,i⋆​\(x\)∈ℝRl×Mi\\mathbf\{A\}\_\{l,i\}^\{\\star\}\(x\)\\in\\mathbb\{R\}^\{R\_\{l\}\\times M\_\{i\}\}denote the rows of the interpolation matrices in Proposition[2](https://arxiv.org/html/2609.19337#Thmproposition2)evaluated at the scalar inputxx\(one row per latent GP coordinate, as in \([30](https://arxiv.org/html/2609.19337#S3.E30)\) at a single test location\)\. Then

𝝁^i​\(x\)=Bg​𝐡^g​\(x\)\+Bδ,i​𝐡^δ,i​\(x\)\+Bl,i​𝐡^l,i​\(x\),\\widehat\{\\bm\{\\mu\}\}\_\{i\}\(x\)=B\_\{g\}\\,\\widehat\{\\mathbf\{h\}\}\_\{g\}\(x\)\+B\_\{\\delta,i\}\\,\\widehat\{\\mathbf\{h\}\}\_\{\\delta,i\}\(x\)\+B\_\{l,i\}\\,\\widehat\{\\mathbf\{h\}\}\_\{l,i\}\(x\),\(38\)Σ^i​\(x\)=Bg​𝐒^g​\(x\)​Bg⊤\+Bδ,i​𝐒^δ,i​\(x\)​Bδ,i⊤\+Bl,i​𝐒^l,i​\(x\)​Bl,i⊤\+𝐑^i​\(x\)\+σ2​IQ,\\widehat\{\\Sigma\}\_\{i\}\(x\)=B\_\{g\}\\,\\widehat\{\\mathbf\{S\}\}\_\{g\}\(x\)\\,B\_\{g\}^\{\\top\}\+B\_\{\\delta,i\}\\,\\widehat\{\\mathbf\{S\}\}\_\{\\delta,i\}\(x\)\\,B\_\{\\delta,i\}^\{\\top\}\+B\_\{l,i\}\\,\\widehat\{\\mathbf\{S\}\}\_\{l,i\}\(x\)\\,B\_\{l,i\}^\{\\top\}\+\\widehat\{\\mathbf\{R\}\}\_\{i\}\(x\)\+\\sigma^\{2\}I\_\{Q\},\(39\)with inducing\-posterior latent summaries𝐡^g​\(x\)=𝐀g⋆​\(x\)​𝐦g\\widehat\{\\mathbf\{h\}\}\_\{g\}\(x\)=\\mathbf\{A\}\_\{g\}^\{\\star\}\(x\)\\mathbf\{m\}\_\{g\},𝐡^δ,i​\(x\)=𝐀δ,i⋆​\(x\)​𝐦δi\\widehat\{\\mathbf\{h\}\}\_\{\\delta,i\}\(x\)=\\mathbf\{A\}\_\{\\delta,i\}^\{\\star\}\(x\)\\mathbf\{m\}\_\{\\delta\_\{i\}\},𝐡^l,i​\(x\)=𝐀l,i⋆​\(x\)​𝐦ui\\widehat\{\\mathbf\{h\}\}\_\{l,i\}\(x\)=\\mathbf\{A\}\_\{l,i\}^\{\\star\}\(x\)\\mathbf\{m\}\_\{u\_\{i\}\}, and𝐒^g​\(x\)=𝐀g⋆​\(x\)​𝐒g​\{𝐀g⋆​\(x\)\}⊤\\widehat\{\\mathbf\{S\}\}\_\{g\}\(x\)=\\mathbf\{A\}\_\{g\}^\{\\star\}\(x\)\\mathbf\{S\}\_\{g\}\\\{\\mathbf\{A\}\_\{g\}^\{\\star\}\(x\)\\\}^\{\\top\},𝐒^δ,i​\(x\)=𝐀δ,i⋆​\(x\)​𝐒δi​\{𝐀δ,i⋆​\(x\)\}⊤\\widehat\{\\mathbf\{S\}\}\_\{\\delta,i\}\(x\)=\\mathbf\{A\}\_\{\\delta,i\}^\{\\star\}\(x\)\\mathbf\{S\}\_\{\\delta\_\{i\}\}\\\{\\mathbf\{A\}\_\{\\delta,i\}^\{\\star\}\(x\)\\\}^\{\\top\},𝐒^l,i​\(x\)=𝐀l,i⋆​\(x\)​𝐒ui​\{𝐀l,i⋆​\(x\)\}⊤\\widehat\{\\mathbf\{S\}\}\_\{l,i\}\(x\)=\\mathbf\{A\}\_\{l,i\}^\{\\star\}\(x\)\\mathbf\{S\}\_\{u\_\{i\}\}\\\{\\mathbf\{A\}\_\{l,i\}^\{\\star\}\(x\)\\\}^\{\\top\}\. Here𝐒^⋅​\(x\)\\widehat\{\\mathbf\{S\}\}\_\{\\cdot\}\(x\)is the*inducing\-posterior*covariance contribution in latent space \(analogous to the𝐀⋅,i⋆​𝐒⋅​\(𝐀⋅,i⋆\)⊤\\mathbf\{A\}\_\{\\cdot,i\}^\{\\star\}\\mathbf\{S\}\_\{\\cdot\}\(\\mathbf\{A\}\_\{\\cdot,i\}^\{\\star\}\)^\{\\top\}blocks in \([30](https://arxiv.org/html/2609.19337#S3.E30)\)\), not the full latent predictive covariance\. The term𝐑^i​\(x\)∈ℝQ×Q\\widehat\{\\mathbf\{R\}\}\_\{i\}\(x\)\\in\\mathbb\{R\}^\{Q\\times Q\}is the VFE residual covariance atxx—uncertainty not explained by the inducing variables, with the same additive global/deviation/local decomposition as in \([10](https://arxiv.org/html/2609.19337#S3.E10)\)—mapped to the output channels through\{Bg,Bδ,i,Bl,i\}\\\{B\_\{g\},B\_\{\\delta,i\},B\_\{l,i\}\\\}\. Adding𝐑^i​\(x\)\\widehat\{\\mathbf\{R\}\}\_\{i\}\(x\)therefore does not double\-count the inducing\-posterior terms carried by𝐒^⋅​\(x\)\\widehat\{\\mathbf\{S\}\}\_\{\\cdot\}\(x\)\.

### 4\.3Assumption\-mismatch Scenarios

We evaluate one matched scenario and three deliberately misspecified scenarios\. These scenarios isolate different ways in which the fitted pFedHGP assumptions may fail in practice\.

##### Scenario A: matched hierarchy\.

The data follow the additive structure in \([36](https://arxiv.org/html/2609.19337#S4.E36)\)\. The global, deviation, and local component processes are mutually independent; the deviation process is smooth; and the local process has short lengthscale\. This scenario tests whether the fitted hierarchy can recover the intended multi\-output decomposition when the modeling assumptions are approximately correct\.

##### Scenario B: deviation\-kernel mismatch\.

The fitted model retains an RBF deviation component, but the true deviation process contains an additional periodic structure\. Specifically, the deviation kernel is

kδtrue​\(x,x′\)=0\.65​τδ2​exp⁡\{−\(x−x′\)22​ℓδ2\}\+0\.35​τδ2​exp⁡\{−2​sin2⁡\{π​\|x−x′\|/3\}0\.82\}\.k\_\{\\delta\}^\{\\mathrm\{true\}\}\(x,x^\{\\prime\}\)=0\.65\\,\\tau\_\{\\delta\}^\{2\}\\exp\\left\\\{\-\\frac\{\(x\-x^\{\\prime\}\)^\{2\}\}\{2\\ell\_\{\\delta\}^\{2\}\}\\right\\\}\+0\.35\\,\\tau\_\{\\delta\}^\{2\}\\exp\\left\\\{\-\\frac\{2\\sin^\{2\}\\\{\\pi\|x\-x^\{\\prime\}\|/3\\\}\}\{0\.8^\{2\}\}\\right\\\}\.\(40\)This tests robustness to an incorrect functional form for systematic client drift\.

##### Scenario C: non\-additive interaction\.

The true clean signal includes a multiplicative interaction between the shared global component and the client\-deviation component:

𝐲iclean​\(x\)=𝐟g​\(x\)\+𝐟δ,i​\(x\)\+𝐟i​\(x\)\+γi​\{𝐟g​\(x\)⊙𝐟δ,i​\(x\)\},\\mathbf\{y\}^\{\\mathrm\{clean\}\}\_\{i\}\(x\)=\\mathbf\{f\}\_\{g\}\(x\)\+\\mathbf\{f\}\_\{\\delta,i\}\(x\)\+\\mathbf\{f\}\_\{i\}\(x\)\+\\gamma\_\{i\}\\\{\\mathbf\{f\}\_\{g\}\(x\)\\odot\\mathbf\{f\}\_\{\\delta,i\}\(x\)\\\},\(41\)where⊙\\odotdenotes elementwise multiplication andγi∼Unif⁡\(0\.35,0\.60\)\\gamma\_\{i\}\\sim\\operatorname\{Unif\}\(0\.35,0\.60\)\. This setting violates the additive pFedHGP decomposition\.

##### Scenario D: correlated hierarchical components\.

The fitted model assumes independent hierarchical components, whereas the true local residual is partially correlated with the deviation process:

𝐟i​\(x\)=ρ​𝐟δ,i​\(x\)\+1−ρ2​𝐟~i​\(x\),ρ=0\.5,\\mathbf\{f\}\_\{i\}\(x\)=\\rho\\mathbf\{f\}\_\{\\delta,i\}\(x\)\+\\sqrt\{1\-\\rho^\{2\}\}\\,\\widetilde\{\\mathbf\{f\}\}\_\{i\}\(x\),\\qquad\\rho=0\.5,\(42\)where𝐟~i\\widetilde\{\\mathbf\{f\}\}\_\{i\}is an independent short\-lengthscale local GP\. This scenario evaluates sensitivity to dependence between the deviation and local residual levels\.

### 4\.4Compared Models

We compare Full pFedHGP with five alternatives that represent different modeling assumptions and degrees of hierarchical structure\.

##### No\-deviation\.

This ablated variant removes the client\-specific deviation component𝐟δ,i\\mathbf\{f\}\_\{\\delta,i\}and fits𝐟g\+𝐟i\\mathbf\{f\}\_\{g\}\+\\mathbf\{f\}\_\{i\}\. It tests whether systematic client\-specific calibration provides additional value beyond a simpler global\-plus\-local decomposition\.

##### No\-local\.

This ablated variant removes the local residual component𝐟i\\mathbf\{f\}\_\{i\}and fits𝐟g\+𝐟δ,i\\mathbf\{f\}\_\{g\}\+\\mathbf\{f\}\_\{\\delta,i\}\. It evaluates the effect of excluding short\-scale client\-specific residual variation from the hierarchy\.

##### Global\-only\.

This model fits only the shared global component𝐟g\\mathbf\{f\}\_\{g\}, pooling all clients without personalization\. It represents the assumption that all clients follow a common population\-level process\.

##### Local\-only\.

This model fits client\-specific private multi\-output GPs without shared global or systematic deviation structure\. It represents the opposite assumption that each client should be modeled independently, without borrowing information across clients\.

##### pFedGP\.

This baseline uses the personalized federated Gaussian process method proposed by[Achituve et al\. \(2021\)](https://arxiv.org/html/2609.19337#bib.bib1)\. In our multi\-output benchmark, each client maintains a personal multi\-output GP block, and model parameters are synchronized by FedAvg\-style averaging\. Unlike pFedHGP, this baseline does not decompose the response into explicit global, deviation, and local components\.

### 4\.5Evaluation Protocol and Metrics

Each scenario is run for3030independent random seeds\. Each client has5050training observations and is evaluated on a common grid of150150test inputs\. Point prediction is measured by the client\-averaged root mean squared error \(RMSE\):

RMSE=1T​∑i=1T\{1ntest​Q​∑j=1ntest‖𝐲i​\(xj\)−𝝁^i​\(xj\)‖22\}1/2\.\\operatorname\{RMSE\}=\\frac\{1\}\{T\}\\sum\_\{i=1\}^\{T\}\\left\\\{\\frac\{1\}\{n\_\{\\mathrm\{test\}\}Q\}\\sum\_\{j=1\}^\{n\_\{\\mathrm\{test\}\}\}\\\|\\mathbf\{y\}\_\{i\}\(x\_\{j\}\)\-\\widehat\{\\bm\{\\mu\}\}\_\{i\}\(x\_\{j\}\)\\\|\_\{2\}^\{2\}\\right\\\}^\{1/2\}\.\(43\)
Probabilistic fit is measured by the client\-averaged multivariate negative log\-likelihood \(NLL\)\. If𝝁^i​j\\widehat\{\\bm\{\\mu\}\}\_\{ij\}andΣ^i​j\\widehat\{\\Sigma\}\_\{ij\}are the predictive mean and fullQ×QQ\\times Qcovariance block for clientiiat test inputxjx\_\{j\}, then

NLL=1T​∑i=1T1ntest​∑j=1ntest12​\[Q​log⁡\(2​π\)\+log⁡\|Σ^i​j\|\+\{𝐲i​\(xj\)−𝝁^i​j\}⊤​Σ^i​j−1​\{𝐲i​\(xj\)−𝝁^i​j\}\]\.\\operatorname\{NLL\}=\\frac\{1\}\{T\}\\sum\_\{i=1\}^\{T\}\\frac\{1\}\{n\_\{\\mathrm\{test\}\}\}\\sum\_\{j=1\}^\{n\_\{\\mathrm\{test\}\}\}\\frac\{1\}\{2\}\\left\[Q\\log\(2\\pi\)\+\\log\|\\widehat\{\\Sigma\}\_\{ij\}\|\+\\\{\\mathbf\{y\}\_\{i\}\(x\_\{j\}\)\-\\widehat\{\\bm\{\\mu\}\}\_\{ij\}\\\}^\{\\top\}\\widehat\{\\Sigma\}\_\{ij\}^\{\-1\}\\\{\\mathbf\{y\}\_\{i\}\(x\_\{j\}\)\-\\widehat\{\\bm\{\\mu\}\}\_\{ij\}\\\}\\right\]\.\(44\)
We also report the Gaussian closed\-form continuous ranked probability score \(CRPS\)\. For outputqq, letμ^i​j​q\\widehat\{\\mu\}\_\{ijq\}andσ^i​j​q2\\widehat\{\\sigma\}\_\{ijq\}^\{2\}denote the marginal predictive mean and variance, and definezi​j​q=\(yi​q​\(xj\)−μ^i​j​q\)/σ^i​j​qz\_\{ijq\}=\(y\_\{iq\}\(x\_\{j\}\)\-\\widehat\{\\mu\}\_\{ijq\}\)/\\widehat\{\\sigma\}\_\{ijq\}\. Then

CRPS=1T​ntest​Q​∑i=1T∑j=1ntest∑q=1Qσ^i​j​q​\[zi​j​q​\{2​Φ​\(zi​j​q\)−1\}\+2​φ​\(zi​j​q\)−1π\],\\operatorname\{CRPS\}=\\frac\{1\}\{Tn\_\{\\mathrm\{test\}\}Q\}\\sum\_\{i=1\}^\{T\}\\sum\_\{j=1\}^\{n\_\{\\mathrm\{test\}\}\}\\sum\_\{q=1\}^\{Q\}\\widehat\{\\sigma\}\_\{ijq\}\\left\[z\_\{ijq\}\\\{2\\Phi\(z\_\{ijq\}\)\-1\\\}\+2\\varphi\(z\_\{ijq\}\)\-\\frac\{1\}\{\\sqrt\{\\pi\}\}\\right\],\(45\)whereΦ⁡\(⋅\)\\Phi\(\\cdot\)andφ⁡\(⋅\)\\varphi\(\\cdot\)denote the standard normal cumulative distribution function and probability density function, respectively \(distinct from the deviation scaleϕ\\phi\)\.

We further report empirical95%95\\%interval coverage and95%95\\%interval width\. Letz0\.95=Φ−1​\(0\.975\)z\_\{0\.95\}=\\Phi^\{\-1\}\(0\.975\)denote the standard normal critical value for a two\-sided95%95\\%interval:

Cov0\.95\\displaystyle\\operatorname\{Cov\}\_\{0\.95\}=1T​ntest​Q∑i=1T∑j=1ntest∑q=1Q𝕀\{\|yi​q\(xj\)−μ^i​j​q\|≤z0\.95σ^i​j​q\},\\displaystyle=\\frac\{1\}\{Tn\_\{\\mathrm\{test\}\}Q\}\\sum\_\{i=1\}^\{T\}\\sum\_\{j=1\}^\{n\_\{\\mathrm\{test\}\}\}\\sum\_\{q=1\}^\{Q\}\\mathbb\{I\}\\left\\\{\|y\_\{iq\}\(x\_\{j\}\)\-\\widehat\{\\mu\}\_\{ijq\}\|\\leq z\_\{0\.95\}\\widehat\{\\sigma\}\_\{ijq\}\\right\\\},\(46\)Width0\.95\\displaystyle\\operatorname\{Width\}\_\{0\.95\}=1T​ntest​Q​∑i=1T∑j=1ntest∑q=1Q2​z0\.95​σ^i​j​q\.\\displaystyle=\\frac\{1\}\{Tn\_\{\\mathrm\{test\}\}Q\}\\sum\_\{i=1\}^\{T\}\\sum\_\{j=1\}^\{n\_\{\\mathrm\{test\}\}\}\\sum\_\{q=1\}^\{Q\}2z\_\{0\.95\}\\widehat\{\\sigma\}\_\{ijq\}\.\(47\)
To evaluate recovery of cross\-output dependence, we first form the ground\-truth output covariance on the common test grid𝒳test\\mathcal\{X\}\_\{\\mathrm\{test\}\}:

Σy,itrue=Covx∈𝒳test⁡\{𝐲iclean​\(x\)\}\+σsim2​IQ\.\\Sigma^\{\\mathrm\{true\}\}\_\{y,i\}=\\operatorname\{Cov\}\_\{x\\in\\mathcal\{X\}\_\{\\mathrm\{test\}\}\}\\bigl\\\{\\mathbf\{y\}^\{\\mathrm\{clean\}\}\_\{i\}\(x\)\\bigr\\\}\+\\sigma\_\{\\mathrm\{sim\}\}^\{2\}I\_\{Q\}\.\(48\)For each fitted method, predictive distributions are joint Gaussians over theQ=4Q=4channels at every test input\. Let𝝁^i​j∈ℝQ\\widehat\{\\bm\{\\mu\}\}\_\{ij\}\\in\\mathbb\{R\}^\{Q\}andΣ^i​j∈ℝQ×Q\\widehat\{\\Sigma\}\_\{ij\}\\in\\mathbb\{R\}^\{Q\\times Q\}denote the predictive mean and full covariance block for clientiiat test inputxj∈𝒳testx\_\{j\}\\in\\mathcal\{X\}\_\{\\mathrm\{test\}\}\. We aggregate these blocks into a model\-implied output covariance

Σ^y,i=Covx∈𝒳test⁡\{𝝁^i​\(x\)\}\+1ntest​∑j=1ntestΣ^i​j,\\widehat\{\\Sigma\}\_\{y,i\}=\\operatorname\{Cov\}\_\{x\\in\\mathcal\{X\}\_\{\\mathrm\{test\}\}\}\\bigl\\\{\\widehat\{\\bm\{\\mu\}\}\_\{i\}\(x\)\\bigr\\\}\+\\frac\{1\}\{n\_\{\\mathrm\{test\}\}\}\\sum\_\{j=1\}^\{n\_\{\\mathrm\{test\}\}\}\\widehat\{\\Sigma\}\_\{ij\},\(49\)and report the relative Frobenius error

CovErr=1T​∑i=1T‖Σ^y,i−Σy,itrue‖F‖Σy,itrue‖F\.\\operatorname\{CovErr\}=\\frac\{1\}\{T\}\\sum\_\{i=1\}^\{T\}\\frac\{\\\|\\widehat\{\\Sigma\}\_\{y,i\}\-\\Sigma^\{\\mathrm\{true\}\}\_\{y,i\}\\\|\_\{F\}\}\{\\\|\\Sigma^\{\\mathrm\{true\}\}\_\{y,i\}\\\|\_\{F\}\}\.\(50\)LowerCovErr\\operatorname\{CovErr\}indicates more accurate recovery of cross\-output covariance structure\. Lower NLL, CRPS, RMSE, and width are also better, while coverage should be close to the nominal level\.

### 4\.6Results

Table[2](https://arxiv.org/html/2609.19337#S4.T2)reports the quantitative comparison, and Figure[2](https://arxiv.org/html/2609.19337#S4.F2)summarizes the main accuracy, probabilistic, and covariance\-recovery metrics graphically\. Lower values are better for RMSE, NLL, CRPS, CovErr, and 95% Width, while 95% coverage should be close to the nominal level\. We organize the discussion around three empirical patterns: the value of the full hierarchy, robustness and covariance recovery under misspecification, and the coverage–sharpness behavior of the predictive uncertainty\.

Table 2:Multi\-output performance across matched and assumption\-mismatch scenarios\. Entries report run\-level means with 95% confidence intervals in brackets\.![Refer to caption](https://arxiv.org/html/2609.19337v1/fig_robustness_summary_ci95_1.png)Figure 2:Robustness summary across matched, deviation\-kernel mismatch, non\-additive, and correlated\-component scenarios\. Error bars denote 95% confidence intervals computed from the results of 30 independent runs\. Bars exceeding the plotting range are clipped for readability, with annotations showing their true values\.The first pattern is that the full pFedHGP is consistently the strongest specification\. In the matched scenario, Full pFedHGP gives the best overall combination of point prediction, probabilistic fit, interval sharpness, and covariance recovery\. The ablations show why this improvement is not attributable to a single pooled multi\-output GP\. Removing the deviation component produces a smaller but systematic degradation, indicating that𝐟δ,i\\mathbf\{f\}\_\{\\delta,i\}captures smooth client\-specific calibration beyond a global\-plus\-local model\. Removing the local component produces a much larger loss in point prediction and probabilistic sharpness, indicating that𝐟i\\mathbf\{f\}\_\{i\}is needed to represent short\-scale variation\. Global\-only underperforms because it cannot personalize to heterogeneous clients, while Local\-only and pFedGP improve on Global\-only but do not recover the same balance between shared structure, client drift, local variation behavior, and output covariance\.

![Refer to caption](https://arxiv.org/html/2609.19337v1/fig_multioutput_prediction_seed0_client0_1.png)Figure 3:Representative four\-channel prediction for one client in the matched benchmark\. The model jointly predicts correlated outputs and reports 95% predictive intervals\.Figure[3](https://arxiv.org/html/2609.19337#S4.F3)provides a representative client visualization for the matched multi\-output setting\. The four channels exhibit different amplitudes and local trajectory patterns, reflecting the cross\-output mixing used in the data\-generating process\. This figure is therefore used as a visual check of the multi\-output regression behavior; the recovery of cross\-output covariance is evaluated quantitatively byCovErr\\operatorname\{CovErr\}in Table[2](https://arxiv.org/html/2609.19337#S4.T2), using the ground\-truth and model\-implied covariances in \([48](https://arxiv.org/html/2609.19337#S4.E48)\)–\([50](https://arxiv.org/html/2609.19337#S4.E50)\) and the fitted ILMM–pFedHGP protocol in Section[4\.2](https://arxiv.org/html/2609.19337#S4.SS2)\.

The second pattern is robustness of prediction under misspecification, which we interpret as a representation\-efficiency issue rather than a failure of the function class\. The fitted layers use squared\-exponential kernels, which are universal on compact input domains: continuous targets on the experimental interval can be approximated by an RBF GP even when the true covariance is periodic, the true map is multiplicative, or the latent layers are correlated\. Universality means the misspecified targets remain approximable\. It does not make those representations equally cheap, and it is not a claim that the true kernel, the multiplicative mechanism, or the independent\-layer split is recovered\. The hierarchy is an efficient parameterization when its assumptions hold\. When they fail, the same universal class can still fit the observable mean and covariance, but it uses extra residual variance, shorter effective local lengthscales, and a finite inducing budget less economically\. Point prediction can still be approximated, so RMSE need not collapse; leftover structure is pushed into residual variance, which is why NLL worsens and intervals widen\.

This view matches Table[2](https://arxiv.org/html/2609.19337#S4.T2)\. Full pFedHGP remains the best or near\-best method in Scenarios B–D, yet the mismatch is not free\. In Scenario B, RMSE rises from0\.2230\.223to0\.2560\.256, NLL from−2\.275\-2\.275to−0\.996\-0\.996, and the average 95% interval width from0\.7650\.765to0\.9190\.919; coverage stays near nominal because the intervals widen\. The periodic deviation is approximable by RBF layers, but a scaled copy ofkgk\_\{g\}is an inefficient code for periodicity, so the local residual must absorb the leftover oscillation\. Reduced models have fewer universal components and therefore pay a larger efficiency cost: No\-local and Global\-only fold the unmatched structure into a single kernel\. In Scenario C the product𝐟g⊙𝐟δ,i\\mathbf\{f\}\_\{g\}\\odot\\mathbf\{f\}\_\{\\delta,i\}is a continuous function ofxxand can be approximated by the remaining RBF layers, especiallyfif\_\{i\}; the additive split is simply a less efficient code for a multiplicative interaction, which shows up mainly as a worse NLL \(−2\.275\-2\.275to−1\.923\-1\.923\) rather than a large RMSE gap\. In Scenario D,𝐟i=ρ​𝐟δ,i\+1−ρ2​𝐟~i\\mathbf\{f\}\_\{i\}=\\rho\\mathbf\{f\}\_\{\\delta,i\}\+\\sqrt\{1\-\\rho^\{2\}\}\\widetilde\{\\mathbf\{f\}\}\_\{i\}is already a linear reparameterization of two GPs, so predictive scores need not degrade; the efficiency cost is replaced by an interpretability cost, because the fitted independent\-layer split is not the true correlated attribution\.

Thus, what remains robust is the ability of a richer universal\-kernel hierarchy to approximate the observable predictive distribution\. The price of mismatch is representation efficiency: weaker NLL and wider intervals when the deviation kernel is wrong, and a weaker literal reading offgf\_\{g\},fδ,if\_\{\\delta,i\}, andfif\_\{i\}when additivity or independence fails\.

![Refer to caption](https://arxiv.org/html/2609.19337v1/fig_reliability_diagram_1.png)Figure 4:Uncertainty calibration and sharpness diagnostics for the multi\-output synthetic benchmark\. The left panel plots empirical coverage minus nominal coverage, with the dashed horizontal line denoting ideal calibration\. The middle panel reports expected calibration error \(ECE\), and the right panel reports the average 95% interval width\.The third pattern concerns uncertainty calibration and sharpness\. Across scenarios, the empirical 95% coverages are close to the nominal level for all methods\. Full pFedHGP is not the closest method to nominal coverage, nor does it have the smallest calibration error\. We do not interpret coverage alone as a ranking criterion\. Instead, the relevant comparison is the coverage–sharpness tradeoff together with proper probabilistic scores\. Figure[4](https://arxiv.org/html/2609.19337#S4.F4)makes this tradeoff explicit: the residual calibration curves show that all methods have small deviations from nominal coverage, the ECE panel summarizes the remaining calibration error across nominal levels, and the interval\-width panel compares sharpness at comparable calibration\. Reduced models such as No\-local and Global\-only achieve smaller calibration error by producing wider, more conservative intervals, reflecting their need to absorb unmodeled client heterogeneity into predictive variance\. In contrast, the full pFedHGP explicitly represents shared global structure, client\-specific deviation, and local site\-specific variation, leading to sharper intervals while remaining close to nominal coverage\. Together with the RMSE, NLL, CRPS, CovErr, and width columns in Table[2](https://arxiv.org/html/2609.19337#S4.T2), these diagnostics show that pFedHGP provides a more favorable probabilistic prediction tradeoff\.

## 5Application Case Studies

### 5\.1Manufacturing Monitoring: Press Tonnage and Fault Diagnosis

This case study evaluates pFedHGP as a federated personalized Gaussian process model for real\-time monitoring and fault diagnosis of a multi\-station metal stamping press\. The physical system is a multi\-operation forming line where billets pass through sequential stamping stations\. Four strain\-gauge tonnage sensors are mounted on the press uprights, as illustrated in Figure[5](https://arxiv.org/html/2609.19337#S5.F5)\. Each production cycle generates a four\-channel force–time waveform that serves as the observed monitoring trajectory\([LEI et al\., 2010](https://arxiv.org/html/2609.19337#bib.bib18)\)\. The prediction task is to distinguish normal operation from four distinct fault modes, corresponding to missing parts at different stations, while maintaining calibration to the shared physical press dynamics\. The manufacturing data provide a high\-dimensional waveform benchmark with shared process structure and condition\-specific variation\. Since all cycles are generated by the same press process, the waveforms contain common mechanical patterns induced by the material, tooling, and forming sequence\. Missing\-part conditions alter these patterns in different local regions of the waveform, making the task suitable for evaluating whether pFedHGP can learn shared structure while adapting to condition\-specific behavior\.

The dataset comprises305305normal cycles and6969cycles for each of four fault classes, with12011201time points per channel\. Thus, the normal class has shape4×1201×3054\\times 1201\\times 305, and each fault class has shape4×1201×694\\times 1201\\times 69\. Representative waveforms for the normal and four fault conditions are shown in Figure[6](https://arxiv.org/html/2609.19337#S5.F6)\. The five clients in this federated benchmark are therefore the normal class and the four fault classes\. These clients correspond to operating conditions in the available single\-press dataset, rather than physically isolated factories or production lines\.

##### Why federation is relevant here\.

Although the available tonnage dataset is collected from a single press system, it provides a prototype for distributed manufacturing monitoring settings in which waveform repositories are generated and owned locally by different machines, lines, plants, suppliers, or production units\. Raw tonnage trajectories may reveal machine condition, tooling configuration, process settings, and fault signatures, and therefore may be difficult to pool centrally\. The high dimensionality of each cycle, with four channels and12011201time points, also creates communication and storage costs if raw trajectories are continuously transmitted to a central server\. pFedHGP addresses this deployment setting by performing posterior updates locally and communicating model\-level updates tied to the shared representation, while retaining client\-specific deviation and local residual structure on private data\.

![Refer to caption](https://arxiv.org/html/2609.19337v1/tonnage_create.png)Figure 5:Schematic of tonnage monitoring system with four strain\-gauge sensors on the press uprights\.![Refer to caption](https://arxiv.org/html/2609.19337v1/tonnage_curve.png)Figure 6:Average tonnage waveforms under normal operation and four fault conditions, with shaded bands denoting standard deviation across samples\.For data preparation, each cycle is represented directly by its raw four\-channel time indices without filtering or resampling\. Channels are concatenated at the cycle level to form the input representation\. The training set includes the first4040normal samples and1010samples per fault class, while the test set includes the remaining265265normal and5959fault samples per class\. This corresponds to using only13\.77%13\.77\\%of the available data for training\.

The pFedHGP model is configured with a shared inducing grid𝐙g\\mathbf\{Z\}\_\{g\}ofM=81M=81points, coupled with a squared\-exponential kernel parameterized by learnable lengthscale and variance\. Deviationsδi\\delta\_\{i\}for each client are defined on the same global basis, with scaling parameterϕ\\phimodulating their amplitude relative to the global latent vector𝐮g\\mathbf\{u\}\_\{g\}\. Each client also maintains a local inducing grid𝐙i\\mathbf\{Z\}\_\{i\}ofMi=41M\_\{i\}=41points, with its own squared\-exponential kernel to capture higher\-frequency, class\-specific variations\. Initial hyperparameters reflect the observed temporal scales: global lengthscale5050, local lengthscale1010, and output variances1\.01\.0and0\.50\.5, respectively\. Training proceeds for 10 federated rounds\. In each round, clients perform 100 local ELBO steps using Adam with learning rate0\.10\.1, after which the server aggregates client\-wise global gradient contributions and broadcasts the updated global block\. Inducing locations are learnable and updated jointly with variational parameters\.

At inference, given a test cycle\(𝐗⋆,𝐲⋆\)\(\\mathbf\{X\}\_\{\\star\},\\mathbf\{y\}\_\{\\star\}\), clientiiproduces the Gaussian predictive law in Section[3\.8\.1](https://arxiv.org/html/2609.19337#S3.SS8.SSS1), with mean𝝁𝐟i​\(𝐗⋆\)\\bm\{\\mu\}\_\{\\mathbf\{f\}\_\{i\}\}\(\\mathbf\{X\}\_\{\\star\}\)and full covariance𝚺i​\(𝐗⋆\)\\bm\{\\Sigma\}\_\{i\}\(\\mathbf\{X\}\_\{\\star\}\)given by \([30](https://arxiv.org/html/2609.19337#S3.E30)\)\. Classification uses the predictive log\-likelihood \([31](https://arxiv.org/html/2609.19337#S3.E31)\) and the maximum\-evidence rule \([32](https://arxiv.org/html/2609.19337#S3.E32)\) over client models\{normal,fault 1,…,fault 4\}\\\{\\text\{normal\},\\text\{fault 1\},\\dots,\\text\{fault 4\}\\\}\.

We compare pFedHGP with established manufacturing monitoring approaches, including personalized Tucker decomposition\([Hu et al\., 2025](https://arxiv.org/html/2609.19337#bib.bib12)\), hierarchical PCA pipelines\([LEI et al\., 2010](https://arxiv.org/html/2609.19337#bib.bib18)\), tensor discriminant learning\([Ye et al\., 2019](https://arxiv.org/html/2609.19337#bib.bib35)\), CNNs with multilinear PCA features\([Guo et al\., 2022](https://arxiv.org/html/2609.19337#bib.bib11)\), and recurrence\-plot imaging\([Zhou et al\., 2015](https://arxiv.org/html/2609.19337#bib.bib42)\)\. For pFedHGP, we repeat the experiment five times with different random seeds, and report the resulting mean accuracy and standard deviation\. Because the competing methods were reported in separate studies and their full implementations are not available, we treat them as published\-reference baselines\. Accordingly, Table[3](https://arxiv.org/html/2609.19337#S5.T3)reports the training percentages and classification accuracies from the corresponding papers\. Standard deviations are included if they are reported in original paper; otherwise, they are omitted because the original papers did not report the necessary variability information\.

Table[3](https://arxiv.org/html/2609.19337#S5.T3)demonstrates that the proposed federated personalized GP model achieves perfect classification accuracy with only13\.77%13\.77\\%of available data, outperforming even perTucker \(15\.49%15\.49\\%data requirement\) and substantially exceeding traditional centralized methods that require7575–90%90\\%of data\. This data efficiency is critical for distributed manufacturing monitoring environments where: \(1\) fault conditions are rare and labeled training data is expensive, \(2\) privacy or data\-governance constraints may prevent sharing production data across facilities, and \(3\) edge devices have limited storage and computational resources\. The hierarchical decomposition enables the model to rapidly adapt to new fault signatures through the deviation component while maintaining the shared physical press model, supporting proactive maintenance planning and real\-time decision support\. The probabilistic predictive scores produced by pFedHGP can support risk\-aware intervention strategies, although a dedicated calibration analysis for the tonnage case is left for future work\.

Table 3:Classification accuracy on the tonnage dataset\. pFedHGP achieves perfect accuracy with the smallest training fraction\. Baseline results are published\-reference results reported in the corresponding papers\. Standard deviations are included when they can be obtained from the reported accuracy variance\.

### 5\.2Urban Air Quality: Federated Monitoring and Zone Discovery

We next evaluate pFedHGP as a federated personalized Gaussian process model for distributed environmental monitoring and urban zone discovery\. The study comprises twelve air\-quality monitoring stations across Beijing, spanning urban core sites and suburban periphery zones\. Each station collects hourly meteorology and pollutant concentrations and is treated as one client maintaining its own local sensor records\. The prediction and discovery task is to learn shared regional air\-quality structure, capture station\-specific variation, identify distinct pollution behaviors across urban zones, and avoid centralized aggregation of station\-level time series\. To construct a temporal series that captures both diurnal cycles and seasonal trends while maintaining computational efficiency, we subsample measurements at 02:00 and 14:00 daily, yielding approximately 730 observations per year per station\. The observed PM2\.5 concentration time series is modeled using pFedHGP following the federated training procedure in Section[3](https://arxiv.org/html/2609.19337#S3)\.

##### Why federation is essential here\.

This air\-quality monitoring study directly reflects a distributed sensing setting\. Monitoring stations are geographically separated sites, and station\-level meteorological and pollutant histories may be owned, managed, or governed by different agencies, municipal units, or local operators\. Centralizing all raw station records would require continuously transmitting and storing large volumes of temporal sensor data in a single repository, which creates substantial communication and storage costs\. pFedHGP addresses this setting by keeping observations at each station, performing local posterior updates on private data, and communicating model\-level updates needed to learn the shared regional representation\.

Federated model configuration\.The shared component usesM=30M=30inducing points placed uniformly along the annual temporal grid to capture Beijing’s regional seasonal pollution dynamics driven by meteorology and emission patterns\. Each station maintains a local inducing grid of sizeMi=12M\_\{i\}=12to model site\-specific transient variations from local traffic, industrial sources, and microclimate effects\. During federated model synchronization, each station updates its local variational parameters using private sensor streams, while the central server aggregates global parameters representing the shared regional air\-quality model\. This hierarchical decomposition separates city\-wide pollution trends from station\- and zone\-specific behavior, supporting interpretable urban environmental modeling under data locality\.

Urban zone discovery\.After training, each station is characterized by its learned projection operator combining deviation and local components \(Section[3\.8](https://arxiv.org/html/2609.19337#S3.SS8)\), which encodes the station’s pollution behavior relative to the shared regional baseline\. These operators enable unsupervised discovery of urban zones with similar environmental characteristics\. Pairwise dissimilarities computed via Frobenius distance produce a station similarity matrix, which is input to spectral clustering with three groups matching the expected geographical structure: urban core, western suburbs, and eastern suburbs\.

Results\.Figure[7](https://arxiv.org/html/2609.19337#S5.F7)shows the ground\-truth geographical distribution: urban core \(blue\), western suburbs including Changping and Dingling \(orange\), and eastern suburbs including Huairou and Shunyi \(green\)\. Figure[8](https://arxiv.org/html/2609.19337#S5.F8)shows the pFedHGP clustering result using 2013 PM2\.5 measurements\. pFedHGP recovers the three geographical zones: urban stations cluster together despite their physical separation, while western and eastern suburban stations form distinct groups\. This result shows that the learned hierarchical representation captures shared regional variation and station\-specific pollution behavior without using explicit geographical coordinates or centralizing raw station histories\. Such automatic zone discovery supports distributed environmental monitoring tasks including targeted pollution analysis, sensor\-network assessment, and air\-quality modeling that accounts for both regional patterns and local urban effects\. The federated architecture keeps station\-level measurements local while still enabling city\-wide coordination through shared model updates\.

![Refer to caption](https://arxiv.org/html/2609.19337v1/map1.png)Figure 7:True geographical distribution of Beijing monitoring stations\. Colors indicate ground\-truth categories: Urban \(blue\), West Suburb \(orange\), East Suburb \(green\)\.![Refer to caption](https://arxiv.org/html/2609.19337v1/clustering_result1.png)Figure 8:pFedHGP clustering result on PM2\.5 measurements \(2013, 02:00 & 14:00 only\)\. Colors denote spectral clustering output: Urban \(blue\), West Suburb \(orange\), East Suburb \(green\)\. The grouping matches the true geographical distribution\.

## 6Limitations and Future Work

Despite its advantages, pFedHGP has several limitations\. First, the current federated protocol keeps raw observations local but does not by itself provide formal differential privacy or secure aggregation guarantees, so model\-update leakage remains a possible risk\. Second, robustness to malicious or corrupted clients is not explicitly handled\. Third, becausekδ=ϕ​kgk\_\{\\delta\}=\\phi k\_\{g\}, the split betweenfgf\_\{g\}andfδ,if\_\{\\delta,i\}is not uniquely identifiable from a single client’s marginal process, and the hierarchical attribution can become less stable when the number of clients is small, when clients have highly unbalanced sample sizes, or when cross\-client heterogeneity is weak\. Fourth, performance depends on kernel and inducing\-budget choices; overly restrictive kernels may underfit systematic deviations, whereas overly flexible local kernels may absorb shared structure\. Fifth, communication cost grows with the number of inducing variables and output factors, so adaptive inducing\-budget selection is an important direction for larger deployments\. Future work will study the conditions under which this hierarchical attribution is stable, including the number of clients, sample\-size balance, and degree of cross\-client heterogeneity, as well as differentially private updates, online updating, robust aggregation under adversarial participation, and scalable inducing\-budget selection in multi\-output federated settings\.

## 7Conclusion

We introduced pFedHGP, a federated personalized Gaussian\-process framework for probabilistic modeling across heterogeneous clients with data that remain locally held\. Its hierarchical construction combines a shared global process, a structured client\-specific deviation defined on the shared basis, and a flexible local residual process\. This organization supplies a regularized representation of common and client\-specific variation while sparse variational inference and federated aggregation allow the shared representation to be learned without centralizing raw observations\.

The ILMM formulation places pFedHGP within a standard multi\-output GP construction and provides a principled representation of cross\-sensor dependence\. Our multi\-output numerical study evaluates the resulting predictive distributions through point\-prediction, proper probabilistic, calibration, sharpness, and output\-covariance metrics\. Across the matched setting, reduced hierarchical variants, and the specified departures from the fitted assumptions, the full hierarchy provides the most favorable overall predictive trade\-off\. These results concern recovery of the observable predictive mean and covariance, rather than exact identification of every latent component\.

The application studies illustrate complementary distributed prediction tasks\. In the manufacturing study, pFedHGP supports fault diagnosis from high\-dimensional multi\-sensor tonnage waveforms using a limited labeled training set\. In the air\-quality study, it learns a shared regional representation while retaining station\-level records locally and supports the discovery of geographically meaningful station groups without explicit coordinates\. Together, these studies demonstrate how federated hierarchical GP modeling can combine information across heterogeneous clients while preserving local adaptation and uncertainty\-aware prediction\.

## Appendix ADerivation of pFedHGP Generative Model

We detail the construction of the interpolation matrices𝐀g,i\\mathbf\{A\}\_\{g,i\}and𝐀i\\mathbf\{A\}\_\{i\}appearing in Proposition[2](https://arxiv.org/html/2609.19337#Thmproposition2)\. The global functionfg​\(⋅\)f\_\{g\}\(\\cdot\)and client\-specific deviationfδ,i​\(⋅\)f\_\{\\delta,i\}\(\\cdot\)share the same kernelkgk\_\{g\}and inducing inputs𝐙g\\mathbf\{Z\}\_\{g\}\. They are approximated via standard kernel interpolation \(sparse GP\) formulas:

fg​\(𝐱\)=𝐊g​\(𝐱,𝐙g\)​𝐊g​\(𝐙g,𝐙g\)−1​𝐮g,fδ,i​\(𝐱\)=𝐊g​\(𝐱,𝐙g\)​𝐊g​\(𝐙g,𝐙g\)−1​𝜹i\.f\_\{g\}\(\\mathbf\{x\}\)=\\mathbf\{K\}\_\{g\}\(\\mathbf\{x\},\\mathbf\{Z\}\_\{g\}\)\\,\\mathbf\{K\}\_\{g\}\(\\mathbf\{Z\}\_\{g\},\\mathbf\{Z\}\_\{g\}\)^\{\-1\}\\mathbf\{u\}\_\{g\},\\quad f\_\{\\delta,i\}\(\\mathbf\{x\}\)=\\mathbf\{K\}\_\{g\}\(\\mathbf\{x\},\\mathbf\{Z\}\_\{g\}\)\\,\\mathbf\{K\}\_\{g\}\(\\mathbf\{Z\}\_\{g\},\\mathbf\{Z\}\_\{g\}\)^\{\-1\}\\bm\{\\delta\}\_\{i\}\.Summing these yields the total global contribution for clientii:fg,i​\(𝐱\)=𝐊g​\(𝐱,𝐙g\)​𝐊g​g−1​\(𝐮g\+𝜹i\)f\_\{g,i\}\(\\mathbf\{x\}\)=\\mathbf\{K\}\_\{g\}\(\\mathbf\{x\},\\mathbf\{Z\}\_\{g\}\)\\,\\mathbf\{K\}\_\{gg\}^\{\-1\}\(\\mathbf\{u\}\_\{g\}\+\\bm\{\\delta\}\_\{i\}\), where𝐊g​g=𝐊g​\(𝐙g,𝐙g\)\\mathbf\{K\}\_\{gg\}=\\mathbf\{K\}\_\{g\}\(\\mathbf\{Z\}\_\{g\},\\mathbf\{Z\}\_\{g\}\)\.

Similarly, the local residualfi​\(⋅\)f\_\{i\}\(\\cdot\)uses a client\-specific kernelkik\_\{i\}and inducing inputs𝐙i\\mathbf\{Z\}\_\{i\}:fi​\(𝐱\)=𝐊i​\(𝐱,𝐙i\)​𝐊i​i−1​𝐮if\_\{i\}\(\\mathbf\{x\}\)=\\mathbf\{K\}\_\{i\}\(\\mathbf\{x\},\\mathbf\{Z\}\_\{i\}\)\\,\\mathbf\{K\}\_\{ii\}^\{\-1\}\\mathbf\{u\}\_\{i\}, where𝐊i​i=𝐊i​\(𝐙i,𝐙i\)\\mathbf\{K\}\_\{ii\}=\\mathbf\{K\}\_\{i\}\(\\mathbf\{Z\}\_\{i\},\\mathbf\{Z\}\_\{i\}\)\.

For a dataset𝐗i\\mathbf\{X\}\_\{i\}, the projection matrices are:

𝐀g,i=𝐊g​\(𝐗i,𝐙g\)​𝐊g​g−1,𝐀i=𝐊i​\(𝐗i,𝐙i\)​𝐊i​i−1\.\\mathbf\{A\}\_\{g,i\}=\\mathbf\{K\}\_\{g\}\(\\mathbf\{X\}\_\{i\},\\mathbf\{Z\}\_\{g\}\)\\,\\mathbf\{K\}\_\{gg\}^\{\-1\},\\quad\\mathbf\{A\}\_\{i\}=\\mathbf\{K\}\_\{i\}\(\\mathbf\{X\}\_\{i\},\\mathbf\{Z\}\_\{i\}\)\\,\\mathbf\{K\}\_\{ii\}^\{\-1\}\.Substituting these into the additive modely=fg\+fδ,i\+fi\+εy=f\_\{g\}\+f\_\{\\delta,i\}\+f\_\{i\}\+\\varepsilonyields Eq\. \([3](https://arxiv.org/html/2609.19337#S3.E3)\)\.

## Appendix BProof of Variational Free Energy ELBO \(Proposition[1](https://arxiv.org/html/2609.19337#Thmproposition1)\)

We derive the Variational Free Energy \(VFE\) ELBO following[Titsias \(2009\)](https://arxiv.org/html/2609.19337#bib.bib29)\. Starting from the standard GP model with observations𝐲=𝐟\+𝜺\\mathbf\{y\}=\\mathbf\{f\}\+\\bm\{\\varepsilon\}where𝜺∼𝒩⁡\(𝟎,σ2​𝐈\)\\bm\{\\varepsilon\}\\sim\\mathcal\{N\}\(\\mathbf\{0\},\\sigma^\{2\}\\mathbf\{I\}\), we introduce inducing variables𝐮\\mathbf\{u\}at locations𝐙\\mathbf\{Z\}\.

The variational distribution factorizes as:

q⁡\(𝐟,𝐮\)=p⁡\(𝐟∣𝐮\)​q​\(𝐮\),q⁡\(𝐮\)=𝒩⁡\(𝐦,𝐒\),q\(\\mathbf\{f\},\\mathbf\{u\}\)=p\(\\mathbf\{f\}\\mid\\mathbf\{u\}\)q\(\\mathbf\{u\}\),\\quad q\(\\mathbf\{u\}\)=\\mathcal\{N\}\(\\mathbf\{m\},\\mathbf\{S\}\),wherep⁡\(𝐟∣𝐮\)=𝒩⁡\(𝐊n​m​𝐊m​m−1​𝐮,𝐊n​n−𝐐n​n\)p\(\\mathbf\{f\}\\mid\\mathbf\{u\}\)=\\mathcal\{N\}\(\\mathbf\{K\}\_\{nm\}\\mathbf\{K\}\_\{mm\}^\{\-1\}\\mathbf\{u\},\\mathbf\{K\}\_\{nn\}\-\\mathbf\{Q\}\_\{nn\}\)is the exact conditional distribution\.

The ELBO is defined as:

ℒVFE=𝔼q⁡\(𝐟,𝐮\)\[logp\(𝐲∣𝐟\)\]−KL\(q\(𝐟,𝐮\)∥p\(𝐟,𝐮\)\)\.\\mathcal\{L\}\_\{\\mathrm\{VFE\}\}=\\mathbb\{E\}\_\{q\(\\mathbf\{f\},\\mathbf\{u\}\)\}\[\\log p\(\\mathbf\{y\}\\mid\\mathbf\{f\}\)\]\-\\mathrm\{KL\}\(q\(\\mathbf\{f\},\\mathbf\{u\}\)\\\|p\(\\mathbf\{f\},\\mathbf\{u\}\)\)\.
Expanding the KL divergence term:

KL\(q\(𝐟,𝐮\)∥p\(𝐟,𝐮\)\)\\displaystyle\\mathrm\{KL\}\(q\(\\mathbf\{f\},\\mathbf\{u\}\)\\\|p\(\\mathbf\{f\},\\mathbf\{u\}\)\)=KL\(p\(𝐟∣𝐮\)q\(𝐮\)∥p\(𝐟∣𝐮\)p\(𝐮\)\)\\displaystyle=\\mathrm\{KL\}\(p\(\\mathbf\{f\}\\mid\\mathbf\{u\}\)q\(\\mathbf\{u\}\)\\\|p\(\\mathbf\{f\}\\mid\\mathbf\{u\}\)p\(\\mathbf\{u\}\)\)=𝔼q⁡\(𝐮\)\[KL\(p\(𝐟∣𝐮\)∥p\(𝐟∣𝐮\)\)\]\+KL\(q\(𝐮\)∥p\(𝐮\)\)\\displaystyle=\\mathbb\{E\}\_\{q\(\\mathbf\{u\}\)\}\[\\mathrm\{KL\}\(p\(\\mathbf\{f\}\\mid\\mathbf\{u\}\)\\\|p\(\\mathbf\{f\}\\mid\\mathbf\{u\}\)\)\]\+\\mathrm\{KL\}\(q\(\\mathbf\{u\}\)\\\|p\(\\mathbf\{u\}\)\)=KL\(q\(𝐮\)∥p\(𝐮\)\),\\displaystyle=\\mathrm\{KL\}\(q\(\\mathbf\{u\}\)\\\|p\(\\mathbf\{u\}\)\),where the first term vanishes sincep⁡\(𝐟∣𝐮\)p\(\\mathbf\{f\}\\mid\\mathbf\{u\}\)appears in both distributions\.

For the expected log\-likelihood term, we have:

𝔼q⁡\(𝐟,𝐮\)​\[log⁡p⁡\(𝐲∣𝐟\)\]\\displaystyle\\mathbb\{E\}\_\{q\(\\mathbf\{f\},\\mathbf\{u\}\)\}\[\\log p\(\\mathbf\{y\}\\mid\\mathbf\{f\}\)\]=𝔼q⁡\(𝐮\)​\[𝔼p⁡\(𝐟∣𝐮\)​\[log⁡p⁡\(𝐲∣𝐟\)\]\]\\displaystyle=\\mathbb\{E\}\_\{q\(\\mathbf\{u\}\)\}\\left\[\\mathbb\{E\}\_\{p\(\\mathbf\{f\}\\mid\\mathbf\{u\}\)\}\[\\log p\(\\mathbf\{y\}\\mid\\mathbf\{f\}\)\]\\right\]=𝔼q⁡\(𝐮\)​\[log⁡p⁡\(𝐲∣𝐮\)\]−12​σ2​Tr⁡\(𝐊n​n−𝐐n​n\),\\displaystyle=\\mathbb\{E\}\_\{q\(\\mathbf\{u\}\)\}\[\\log p\(\\mathbf\{y\}\\mid\\mathbf\{u\}\)\]\-\\frac\{1\}\{2\\sigma^\{2\}\}\\operatorname\{Tr\}\(\\mathbf\{K\}\_\{nn\}\-\\mathbf\{Q\}\_\{nn\}\),where the trace term arises from the difference between the full covariance𝐊n​n\\mathbf\{K\}\_\{nn\}and the low\-rank approximation𝐐n​n=𝐊n​m​𝐊m​m−1​𝐊m​n\\mathbf\{Q\}\_\{nn\}=\\mathbf\{K\}\_\{nm\}\\mathbf\{K\}\_\{mm\}^\{\-1\}\\mathbf\{K\}\_\{mn\}in the conditional distributionp⁡\(𝐟∣𝐮\)p\(\\mathbf\{f\}\\mid\\mathbf\{u\}\)\.

Combining these terms yields the VFE ELBO:

ℒVFE=𝔼q⁡\(𝐮\)\[logp\(𝐲∣𝐮\)\]−KL\(q\(𝐮\)∥p\(𝐮\)\)−12​σ2Tr\(𝐊n​n−𝐐n​n\),\\mathcal\{L\}\_\{\\mathrm\{VFE\}\}=\\mathbb\{E\}\_\{q\(\\mathbf\{u\}\)\}\[\\log p\(\\mathbf\{y\}\\mid\\mathbf\{u\}\)\]\-\\mathrm\{KL\}\(q\(\\mathbf\{u\}\)\\\|p\(\\mathbf\{u\}\)\)\-\\frac\{1\}\{2\\sigma^\{2\}\}\\operatorname\{Tr\}\(\\mathbf\{K\}\_\{nn\}\-\\mathbf\{Q\}\_\{nn\}\),which completes the proof\.□\\square

## Appendix CProof of ELBO Decomposition \(Theorem[1](https://arxiv.org/html/2609.19337#Thmtheorem1)\)

We provide the derivation of the ELBO decomposition\. The variational lower bound on the marginal likelihood is defined as:

ℒ=𝔼q​\[log⁡p⁡\(\{𝐲i\},𝐮g,\{𝜹i\},\{𝐮i\}\)\]−𝔼q​\[log⁡q⁡\(𝐮g,\{𝜹i\},\{𝐮i\}\)\],\\mathcal\{L\}=\\mathbb\{E\}\_\{q\}\[\\log p\(\\\{\\mathbf\{y\}\_\{i\}\\\},\\mathbf\{u\}\_\{g\},\\\{\\bm\{\\delta\}\_\{i\}\\\},\\\{\\mathbf\{u\}\_\{i\}\\\}\)\]\-\\mathbb\{E\}\_\{q\}\[\\log q\(\\mathbf\{u\}\_\{g\},\\\{\\bm\{\\delta\}\_\{i\}\\\},\\\{\\mathbf\{u\}\_\{i\}\\\}\)\],where the expectation is taken with respect to the variational posteriorq⁡\(𝐮g,\{𝜹i\},\{𝐮i\}\)=q⁡\(𝐮g\)​∏i=1Tq⁡\(𝜹i\)​q​\(𝐮i\)q\(\\mathbf\{u\}\_\{g\},\\\{\\bm\{\\delta\}\_\{i\}\\\},\\\{\\mathbf\{u\}\_\{i\}\\\}\)=q\(\\mathbf\{u\}\_\{g\}\)\\prod\_\{i=1\}^\{T\}q\(\\bm\{\\delta\}\_\{i\}\)q\(\\mathbf\{u\}\_\{i\}\)\.

Substituting the factorized joint distributionp⁡\(\{𝐲i\},𝐮g,\{𝜹i\},\{𝐮i\}\)=p⁡\(𝐮g\)​∏i=1Tp⁡\(𝐲i∣𝐮g,𝜹i,𝐮i\)​p​\(𝜹i\)​p​\(𝐮i\)p\(\\\{\\mathbf\{y\}\_\{i\}\\\},\\mathbf\{u\}\_\{g\},\\\{\\bm\{\\delta\}\_\{i\}\\\},\\\{\\mathbf\{u\}\_\{i\}\\\}\)=p\(\\mathbf\{u\}\_\{g\}\)\\prod\_\{i=1\}^\{T\}p\(\\mathbf\{y\}\_\{i\}\\mid\\mathbf\{u\}\_\{g\},\\bm\{\\delta\}\_\{i\},\\mathbf\{u\}\_\{i\}\)p\(\\bm\{\\delta\}\_\{i\}\)p\(\\mathbf\{u\}\_\{i\}\), we can expand the first term:

𝔼q​\[log⁡p\]\\displaystyle\\mathbb\{E\}\_\{q\}\[\\log p\]=𝔼q​\[log⁡p⁡\(𝐮g\)\+∑i=1T\(log⁡p⁡\(𝐲i∣𝐮g,𝜹i,𝐮i\)\+log⁡p⁡\(𝜹i\)\+log⁡p⁡\(𝐮i\)\)\]\\displaystyle=\\mathbb\{E\}\_\{q\}\\left\[\\log p\(\\mathbf\{u\}\_\{g\}\)\+\\sum\_\{i=1\}^\{T\}\\left\(\\log p\(\\mathbf\{y\}\_\{i\}\\mid\\mathbf\{u\}\_\{g\},\\bm\{\\delta\}\_\{i\},\\mathbf\{u\}\_\{i\}\)\+\\log p\(\\bm\{\\delta\}\_\{i\}\)\+\\log p\(\\mathbf\{u\}\_\{i\}\)\\right\)\\right\]=𝔼q⁡\(𝐮g\)​\[log⁡p⁡\(𝐮g\)\]\+∑i=1T\(𝔼⁡\[log⁡p⁡\(𝐲i∣𝐮g,𝜹i,𝐮i\)\]\+𝔼⁡\[log⁡p⁡\(𝜹i\)\]\+𝔼⁡\[log⁡p⁡\(𝐮i\)\]\)\.\\displaystyle=\\mathbb\{E\}\_\{q\(\\mathbf\{u\}\_\{g\}\)\}\[\\log p\(\\mathbf\{u\}\_\{g\}\)\]\+\\sum\_\{i=1\}^\{T\}\\left\(\\mathbb\{E\}\[\\log p\(\\mathbf\{y\}\_\{i\}\\mid\\mathbf\{u\}\_\{g\},\\bm\{\\delta\}\_\{i\},\\mathbf\{u\}\_\{i\}\)\]\+\\mathbb\{E\}\[\\log p\(\\bm\{\\delta\}\_\{i\}\)\]\+\\mathbb\{E\}\[\\log p\(\\mathbf\{u\}\_\{i\}\)\]\\right\)\.
Similarly, the entropy term expands as:

𝔼q​\[log⁡q\]\\displaystyle\\mathbb\{E\}\_\{q\}\[\\log q\]=𝔼q​\[log⁡q⁡\(𝐮g\)\+∑i=1T\(log⁡q⁡\(𝜹i\)\+log⁡q⁡\(𝐮i\)\)\]\\displaystyle=\\mathbb\{E\}\_\{q\}\\left\[\\log q\(\\mathbf\{u\}\_\{g\}\)\+\\sum\_\{i=1\}^\{T\}\\left\(\\log q\(\\bm\{\\delta\}\_\{i\}\)\+\\log q\(\\mathbf\{u\}\_\{i\}\)\\right\)\\right\]=𝔼q⁡\(𝐮g\)​\[log⁡q⁡\(𝐮g\)\]\+∑i=1T\(𝔼⁡\[log⁡q⁡\(𝜹i\)\]\+𝔼⁡\[log⁡q⁡\(𝐮i\)\]\)\.\\displaystyle=\\mathbb\{E\}\_\{q\(\\mathbf\{u\}\_\{g\}\)\}\[\\log q\(\\mathbf\{u\}\_\{g\}\)\]\+\\sum\_\{i=1\}^\{T\}\\left\(\\mathbb\{E\}\[\\log q\(\\bm\{\\delta\}\_\{i\}\)\]\+\\mathbb\{E\}\[\\log q\(\\mathbf\{u\}\_\{i\}\)\]\\right\)\.
Combining terms by variable groups:

- •Global terms:𝔼q⁡\(𝐮g\)\[logp\(𝐮g\)\]−𝔼q⁡\(𝐮g\)\[logq\(𝐮g\)\]=−KL\(q\(𝐮g\)∥p\(𝐮g\)\)\\mathbb\{E\}\_\{q\(\\mathbf\{u\}\_\{g\}\)\}\[\\log p\(\\mathbf\{u\}\_\{g\}\)\]\-\\mathbb\{E\}\_\{q\(\\mathbf\{u\}\_\{g\}\)\}\[\\log q\(\\mathbf\{u\}\_\{g\}\)\]=\-\\mathrm\{KL\}\(q\(\\mathbf\{u\}\_\{g\}\)\\\|p\(\\mathbf\{u\}\_\{g\}\)\)\.
- •Local terms for clientii: 𝔼⁡\[log⁡p⁡\(𝐲i∣𝐮g,𝜹i,𝐮i\)\]\+\(𝔼⁡\[log⁡p⁡\(𝜹i\)\]−𝔼⁡\[log⁡q⁡\(𝜹i\)\]\)\+\(𝔼⁡\[log⁡p⁡\(𝐮i\)\]−𝔼⁡\[log⁡q⁡\(𝐮i\)\]\)\\displaystyle\\mathbb\{E\}\[\\log p\(\\mathbf\{y\}\_\{i\}\\mid\\mathbf\{u\}\_\{g\},\\bm\{\\delta\}\_\{i\},\\mathbf\{u\}\_\{i\}\)\]\+\(\\mathbb\{E\}\[\\log p\(\\bm\{\\delta\}\_\{i\}\)\]\-\\mathbb\{E\}\[\\log q\(\\bm\{\\delta\}\_\{i\}\)\]\)\+\(\\mathbb\{E\}\[\\log p\(\\mathbf\{u\}\_\{i\}\)\]\-\\mathbb\{E\}\[\\log q\(\\mathbf\{u\}\_\{i\}\)\]\)=𝔼\[logp\(𝐲i∣𝐮g,𝜹i,𝐮i\)\]−KL\(q\(𝜹i\)∥p\(𝜹i\)\)−KL\(q\(𝐮i\)∥p\(𝐮i\)\)\.\\displaystyle=\\mathbb\{E\}\[\\log p\(\\mathbf\{y\}\_\{i\}\\mid\\mathbf\{u\}\_\{g\},\\bm\{\\delta\}\_\{i\},\\mathbf\{u\}\_\{i\}\)\]\-\\mathrm\{KL\}\(q\(\\bm\{\\delta\}\_\{i\}\)\\\|p\(\\bm\{\\delta\}\_\{i\}\)\)\-\\mathrm\{KL\}\(q\(\\mathbf\{u\}\_\{i\}\)\\\|p\(\\mathbf\{u\}\_\{i\}\)\)\.

Summing these components yields the result in Theorem[1](https://arxiv.org/html/2609.19337#Thmtheorem1):

ℒ=∑i=1T\(𝔼\[logp\(𝐲i∣𝐮g,𝜹i,𝐮i\)\]−KL\(q\(𝜹i\)∥p\(𝜹i\)\)−KL\(q\(𝐮i\)∥p\(𝐮i\)\)\)⏟ℒi−KL\(q\(𝐮g\)∥p\(𝐮g\)\)\.\\mathcal\{L\}=\\sum\_\{i=1\}^\{T\}\\underbrace\{\\left\(\\mathbb\{E\}\[\\log p\(\\mathbf\{y\}\_\{i\}\\mid\\mathbf\{u\}\_\{g\},\\bm\{\\delta\}\_\{i\},\\mathbf\{u\}\_\{i\}\)\]\-\\mathrm\{KL\}\(q\(\\bm\{\\delta\}\_\{i\}\)\\\|p\(\\bm\{\\delta\}\_\{i\}\)\)\-\\mathrm\{KL\}\(q\(\\mathbf\{u\}\_\{i\}\)\\\|p\(\\mathbf\{u\}\_\{i\}\)\)\\right\)\}\_\{\\mathcal\{L\}\_\{i\}\}\-\\mathrm\{KL\}\(q\(\\mathbf\{u\}\_\{g\}\)\\\|p\(\\mathbf\{u\}\_\{g\}\)\)\.
##### Expected log\-likelihood with variational covariance\.

For clientii, the conditional likelihoodp⁡\(𝐲i∣𝐮g,𝜹i,𝐮i\)=𝒩⁡\(𝐀g,i​\(𝐮g\+𝜹i\)\+𝐀i​𝐮i,σ2​𝐈\)p\(\\mathbf\{y\}\_\{i\}\\mid\\mathbf\{u\}\_\{g\},\\bm\{\\delta\}\_\{i\},\\mathbf\{u\}\_\{i\}\)=\\mathcal\{N\}\(\\mathbf\{A\}\_\{g,i\}\(\\mathbf\{u\}\_\{g\}\+\\bm\{\\delta\}\_\{i\}\)\+\\mathbf\{A\}\_\{i\}\\mathbf\{u\}\_\{i\},\\sigma^\{2\}\\mathbf\{I\}\)is linear in the inducing variables\. Withq⁡\(𝐮g\)=𝒩⁡\(𝐦g,𝐒g\)q\(\\mathbf\{u\}\_\{g\}\)=\\mathcal\{N\}\(\\mathbf\{m\}\_\{g\},\\mathbf\{S\}\_\{g\}\),q⁡\(𝜹i\)=𝒩⁡\(𝐦δi,𝐒δi\)q\(\\bm\{\\delta\}\_\{i\}\)=\\mathcal\{N\}\(\\mathbf\{m\}\_\{\\delta\_\{i\}\},\\mathbf\{S\}\_\{\\delta\_\{i\}\}\), andq⁡\(𝐮i\)=𝒩⁡\(𝐦ui,𝐒ui\)q\(\\mathbf\{u\}\_\{i\}\)=\\mathcal\{N\}\(\\mathbf\{m\}\_\{u\_\{i\}\},\\mathbf\{S\}\_\{u\_\{i\}\}\), the predictive mean is𝝁𝐟i\\bm\{\\mu\}\_\{\\mathbf\{f\}\_\{i\}\}in \([17](https://arxiv.org/html/2609.19337#S3.E17)\) and the extra expected\-log\-likelihood terms are

−12​σ2​\(Tr⁡\(𝐀g,i​𝐒g​𝐀g,i⊤\)\+Tr⁡\(𝐀δ,i​𝐒δi​𝐀δ,i⊤\)\+Tr⁡\(𝐀i​𝐒ui​𝐀i⊤\)\),\-\\frac\{1\}\{2\\sigma^\{2\}\}\\Bigl\(\\operatorname\{Tr\}\(\\mathbf\{A\}\_\{g,i\}\\mathbf\{S\}\_\{g\}\\mathbf\{A\}\_\{g,i\}^\{\\top\}\)\+\\operatorname\{Tr\}\(\\mathbf\{A\}\_\{\\delta,i\}\\mathbf\{S\}\_\{\\delta\_\{i\}\}\\mathbf\{A\}\_\{\\delta,i\}^\{\\top\}\)\+\\operatorname\{Tr\}\(\\mathbf\{A\}\_\{i\}\\mathbf\{S\}\_\{u\_\{i\}\}\\mathbf\{A\}\_\{i\}^\{\\top\}\)\\Bigr\),in addition to the residual traceTr⁡\(𝐑i\)\\operatorname\{Tr\}\(\\mathbf\{R\}\_\{i\}\)from the sparse GP approximation\. This yields \([16](https://arxiv.org/html/2609.19337#S3.E16)\) and the explicit ELBO in Theorem[1](https://arxiv.org/html/2609.19337#Thmtheorem1)\.

## Appendix DProof of Proposition[2](https://arxiv.org/html/2609.19337#Thmproposition2)\(Low\-Rank Hierarchical Representation\)

We prove that the conditional mean of the total latent function admits the low\-rank basis expansion stated in Proposition[2](https://arxiv.org/html/2609.19337#Thmproposition2)\.

##### Setup\.

Under the GP priors, each component process is a Gaussian process:

fg​\(⋅\)∼𝒢​𝒫​\(0,kg​\(⋅,⋅\)\),fδ,i​\(⋅\)∼𝒢​𝒫​\(0,ϕ​kg​\(⋅,⋅\)\),fi​\(⋅\)∼𝒢​𝒫​\(0,ki​\(⋅,⋅\)\)\.f\_\{g\}\(\\cdot\)\\sim\\mathcal\{GP\}\(0,k\_\{g\}\(\\cdot,\\cdot\)\),\\quad f\_\{\\delta,i\}\(\\cdot\)\\sim\\mathcal\{GP\}\(0,\\phi k\_\{g\}\(\\cdot,\\cdot\)\),\\quad f\_\{i\}\(\\cdot\)\\sim\\mathcal\{GP\}\(0,k\_\{i\}\(\\cdot,\\cdot\)\)\.The inducing variables𝐮g=fg​\(𝐙g\)\\mathbf\{u\}\_\{g\}=f\_\{g\}\(\\mathbf\{Z\}\_\{g\}\),𝜹i=fδ,i​\(𝐙g\)\\bm\{\\delta\}\_\{i\}=f\_\{\\delta,i\}\(\\mathbf\{Z\}\_\{g\}\), and𝐮i=fi​\(𝐙i\)\\mathbf\{u\}\_\{i\}=f\_\{i\}\(\\mathbf\{Z\}\_\{i\}\)follow the joint Gaussian distributions specified in the proposition\.

##### Conditional Mean Derivation\.

By properties of Gaussian processes, the conditional distribution offg​\(𝐱\)f\_\{g\}\(\\mathbf\{x\}\)given𝐮g\\mathbf\{u\}\_\{g\}is Gaussian with mean

𝔼⁡\[fg​\(𝐱\)∣𝐮g\]=Cov⁡\(fg​\(𝐱\),𝐮g\)​Var​\(𝐮g\)−1​𝐮g=𝐤g​\(𝐱,𝐙g\)⊤​𝐊g​g−1​𝐮g,\\mathbb\{E\}\[f\_\{g\}\(\\mathbf\{x\}\)\\mid\\mathbf\{u\}\_\{g\}\]=\\mathrm\{Cov\}\(f\_\{g\}\(\\mathbf\{x\}\),\\mathbf\{u\}\_\{g\}\)\\mathrm\{Var\}\(\\mathbf\{u\}\_\{g\}\)^\{\-1\}\\mathbf\{u\}\_\{g\}=\\mathbf\{k\}\_\{g\}\(\\mathbf\{x\},\\mathbf\{Z\}\_\{g\}\)^\{\\top\}\\mathbf\{K\}\_\{gg\}^\{\-1\}\\mathbf\{u\}\_\{g\},where𝐤g​\(𝐱,𝐙g\)=\[kg​\(𝐱,𝐳g,1\),…,kg​\(𝐱,𝐳g,M\)\]⊤\\mathbf\{k\}\_\{g\}\(\\mathbf\{x\},\\mathbf\{Z\}\_\{g\}\)=\[k\_\{g\}\(\\mathbf\{x\},\\mathbf\{z\}\_\{g,1\}\),\\ldots,k\_\{g\}\(\\mathbf\{x\},\\mathbf\{z\}\_\{g,M\}\)\]^\{\\top\}is the cross\-covariance vector\.

Similarly, for the client\-specific global component:𝔼⁡\[fδ,i​\(𝐱\)∣𝜹i\]=𝐤g​\(𝐱,𝐙g\)⊤​𝐊g​g−1​𝜹i\\mathbb\{E\}\[f\_\{\\delta,i\}\(\\mathbf\{x\}\)\\mid\\bm\{\\delta\}\_\{i\}\]=\\mathbf\{k\}\_\{g\}\(\\mathbf\{x\},\\mathbf\{Z\}\_\{g\}\)^\{\\top\}\\mathbf\{K\}\_\{gg\}^\{\-1\}\\bm\{\\delta\}\_\{i\}, noting thatfδ,if\_\{\\delta,i\}uses the same kernel structure asfgf\_\{g\}\.

For the local component:𝔼⁡\[fi​\(𝐱\)∣𝐮i\]=𝐤i​\(𝐱,𝐙i\)⊤​𝐊i​i−1​𝐮i\\mathbb\{E\}\[f\_\{i\}\(\\mathbf\{x\}\)\\mid\\mathbf\{u\}\_\{i\}\]=\\mathbf\{k\}\_\{i\}\(\\mathbf\{x\},\\mathbf\{Z\}\_\{i\}\)^\{\\top\}\\mathbf\{K\}\_\{ii\}^\{\-1\}\\mathbf\{u\}\_\{i\}\.

##### Total Mean\.

The total latent function isfitot​\(𝐱\)=fg​\(𝐱\)\+fδ,i​\(𝐱\)\+fi​\(𝐱\)f\_\{i\}^\{\\mathrm\{tot\}\}\(\\mathbf\{x\}\)=f\_\{g\}\(\\mathbf\{x\}\)\+f\_\{\\delta,i\}\(\\mathbf\{x\}\)\+f\_\{i\}\(\\mathbf\{x\}\)\. Since the inducing variables from different components are independent, the conditional expectation of the sum equals the sum of conditional expectations:

𝔼\[fitot\(𝐱\)∣𝐮g,𝜹i,𝐮i\]\\displaystyle\\mathbb\{E\}\[f\_\{i\}^\{\\mathrm\{tot\}\}\(\\mathbf\{x\}\)\\mid\\mathbf\{u\}\_\{g\},\\bm\{\\delta\}\_\{i\},\\mathbf\{u\}\_\{i\}\]=𝔼⁡\[fg​\(𝐱\)∣𝐮g\]\+𝔼⁡\[fδ,i​\(𝐱\)∣𝜹i\]\+𝔼⁡\[fi​\(𝐱\)∣𝐮i\]\\displaystyle=\\mathbb\{E\}\[f\_\{g\}\(\\mathbf\{x\}\)\\mid\\mathbf\{u\}\_\{g\}\]\+\\mathbb\{E\}\[f\_\{\\delta,i\}\(\\mathbf\{x\}\)\\mid\\bm\{\\delta\}\_\{i\}\]\+\\mathbb\{E\}\[f\_\{i\}\(\\mathbf\{x\}\)\\mid\\mathbf\{u\}\_\{i\}\]=𝐤g​\(𝐱,𝐙g\)⊤​𝐊g​g−1​\(𝐮g\+𝜹i\)\+𝐤i​\(𝐱,𝐙i\)⊤​𝐊i​i−1​𝐮i\.\\displaystyle=\\mathbf\{k\}\_\{g\}\(\\mathbf\{x\},\\mathbf\{Z\}\_\{g\}\)^\{\\top\}\\mathbf\{K\}\_\{gg\}^\{\-1\}\(\\mathbf\{u\}\_\{g\}\+\\bm\{\\delta\}\_\{i\}\)\+\\mathbf\{k\}\_\{i\}\(\\mathbf\{x\},\\mathbf\{Z\}\_\{i\}\)^\{\\top\}\\mathbf\{K\}\_\{ii\}^\{\-1\}\\mathbf\{u\}\_\{i\}\.
For the vectorized form at observations𝐗i\\mathbf\{X\}\_\{i\}, define𝐊g​\(𝐗i,𝐙g\)∈ℝni×M\\mathbf\{K\}\_\{g\}\(\\mathbf\{X\}\_\{i\},\\mathbf\{Z\}\_\{g\}\)\\in\\mathbb\{R\}^\{n\_\{i\}\\times M\}as the matrix with entries\[𝐊g​\(𝐗i,𝐙g\)\]j​m=kg​\(𝐱i,j,𝐳g,m\)\[\\mathbf\{K\}\_\{g\}\(\\mathbf\{X\}\_\{i\},\\mathbf\{Z\}\_\{g\}\)\]\_\{jm\}=k\_\{g\}\(\\mathbf\{x\}\_\{i,j\},\\mathbf\{z\}\_\{g,m\}\)\. Then:

𝔼\[𝐟itot∣𝐮g,𝜹i,𝐮i\]\\displaystyle\\mathbb\{E\}\[\\mathbf\{f\}\_\{i\}^\{\\mathrm\{tot\}\}\\mid\\mathbf\{u\}\_\{g\},\\bm\{\\delta\}\_\{i\},\\mathbf\{u\}\_\{i\}\]=𝐊g​\(𝐗i,𝐙g\)​𝐊g​g−1​\(𝐮g\+𝜹i\)\\displaystyle=\\mathbf\{K\}\_\{g\}\(\\mathbf\{X\}\_\{i\},\\mathbf\{Z\}\_\{g\}\)\\mathbf\{K\}\_\{gg\}^\{\-1\}\(\\mathbf\{u\}\_\{g\}\+\\bm\{\\delta\}\_\{i\}\)\+𝐊i​\(𝐗i,𝐙i\)​𝐊i​i−1​𝐮i,\\displaystyle\+\\mathbf\{K\}\_\{i\}\(\\mathbf\{X\}\_\{i\},\\mathbf\{Z\}\_\{i\}\)\\mathbf\{K\}\_\{ii\}^\{\-1\}\\mathbf\{u\}\_\{i\},which is the low\-rank representation stated in \([8](https://arxiv.org/html/2609.19337#S3.E8)\)\.□\\square

## Appendix EProof of Proposition[3](https://arxiv.org/html/2609.19337#Thmproposition3)\(Hierarchical Conditional Covariance Decomposition\)

We prove the additive decomposition of the conditional covariance stated in Proposition[3](https://arxiv.org/html/2609.19337#Thmproposition3)\.

##### Conditional Covariance Formula\.

For a Gaussian processf∼𝒢​𝒫​\(0,k\)f\\sim\\mathcal\{GP\}\(0,k\)with inducing variables𝐮=f⁡\(𝐙\)\\mathbf\{u\}=f\(\\mathbf\{Z\}\), the conditional covariance at observations𝐗\\mathbf\{X\}is given by:

Cov⁡\(f⁡\(𝐗\)∣𝐮\)=k⁡\(𝐗,𝐗\)−k⁡\(𝐗,𝐙\)​k​\(𝐙,𝐙\)−1​k​\(𝐙,𝐗\)\.\\mathrm\{Cov\}\(f\(\\mathbf\{X\}\)\\mid\\mathbf\{u\}\)=k\(\\mathbf\{X\},\\mathbf\{X\}\)\-k\(\\mathbf\{X\},\\mathbf\{Z\}\)k\(\\mathbf\{Z\},\\mathbf\{Z\}\)^\{\-1\}k\(\\mathbf\{Z\},\\mathbf\{X\}\)\.This is the residual covariance after conditioning on the inducing points\.

##### Component\-wise Decomposition\.

For the fixed global component:Cov⁡\(fg​\(𝐗i\)∣𝐮g\)=𝐊g,i​i−𝐐g,i​i\\mathrm\{Cov\}\(f\_\{g\}\(\\mathbf\{X\}\_\{i\}\)\\mid\\mathbf\{u\}\_\{g\}\)=\\mathbf\{K\}\_\{g,ii\}\-\\mathbf\{Q\}\_\{g,ii\}, where𝐊g,i​i=kg​\(𝐗i,𝐗i\)\\mathbf\{K\}\_\{g,ii\}=k\_\{g\}\(\\mathbf\{X\}\_\{i\},\\mathbf\{X\}\_\{i\}\)and𝐐g,i​i=𝐊g​\(𝐗i,𝐙g\)​𝐊g​g−1​𝐊g​\(𝐙g,𝐗i\)\\mathbf\{Q\}\_\{g,ii\}=\\mathbf\{K\}\_\{g\}\(\\mathbf\{X\}\_\{i\},\\mathbf\{Z\}\_\{g\}\)\\mathbf\{K\}\_\{gg\}^\{\-1\}\\mathbf\{K\}\_\{g\}\(\\mathbf\{Z\}\_\{g\},\\mathbf\{X\}\_\{i\}\)\.

For the client\-specific deviation component withkδ=ϕ​kgk\_\{\\delta\}=\\phi k\_\{g\}:𝐊δ,i​i=ϕ​𝐊g,i​i\\mathbf\{K\}\_\{\\delta,ii\}=\\phi\\mathbf\{K\}\_\{g,ii\}and𝐐δ,i​i=ϕ​𝐐g,i​i\\mathbf\{Q\}\_\{\\delta,ii\}=\\phi\\mathbf\{Q\}\_\{g,ii\}, soCov⁡\(fδ,i​\(𝐗i\)∣𝜹i\)=𝐊δ,i​i−𝐐δ,i​i=ϕ⁡\(𝐊g,i​i−𝐐g,i​i\)\\mathrm\{Cov\}\(f\_\{\\delta,i\}\(\\mathbf\{X\}\_\{i\}\)\\mid\\bm\{\\delta\}\_\{i\}\)=\\mathbf\{K\}\_\{\\delta,ii\}\-\\mathbf\{Q\}\_\{\\delta,ii\}=\\phi\(\\mathbf\{K\}\_\{g,ii\}\-\\mathbf\{Q\}\_\{g,ii\}\)\. The scaleϕ\\phitherefore enters the residual covariance of the deviation layer even though the conditional mean shares the same interpolation matrix asfgf\_\{g\}\.

For the local component:Cov⁡\(fi​\(𝐗i\)∣𝐮i\)=𝐊i,i​i−𝐐i,i​i\\mathrm\{Cov\}\(f\_\{i\}\(\\mathbf\{X\}\_\{i\}\)\\mid\\mathbf\{u\}\_\{i\}\)=\\mathbf\{K\}\_\{i,ii\}\-\\mathbf\{Q\}\_\{i,ii\}, where𝐊i,i​i=ki​\(𝐗i,𝐗i\)\\mathbf\{K\}\_\{i,ii\}=k\_\{i\}\(\\mathbf\{X\}\_\{i\},\\mathbf\{X\}\_\{i\}\)and𝐐i,i​i=𝐊i​\(𝐗i,𝐙i\)​𝐊i​i−1​𝐊i​\(𝐙i,𝐗i\)\\mathbf\{Q\}\_\{i,ii\}=\\mathbf\{K\}\_\{i\}\(\\mathbf\{X\}\_\{i\},\\mathbf\{Z\}\_\{i\}\)\\mathbf\{K\}\_\{ii\}^\{\-1\}\\mathbf\{K\}\_\{i\}\(\\mathbf\{Z\}\_\{i\},\\mathbf\{X\}\_\{i\}\)\.

##### Independence and Additivity\.

Under the mean\-field variational approximation, the three componentsfgf\_\{g\},fδ,if\_\{\\delta,i\}, andfif\_\{i\}are conditionally independent given their respective inducing variables\. Therefore, for the total latent function𝐟itot=fg​\(𝐗i\)\+fδ,i​\(𝐗i\)\+fi​\(𝐗i\)\\mathbf\{f\}\_\{i\}^\{\\mathrm\{tot\}\}=f\_\{g\}\(\\mathbf\{X\}\_\{i\}\)\+f\_\{\\delta,i\}\(\\mathbf\{X\}\_\{i\}\)\+f\_\{i\}\(\\mathbf\{X\}\_\{i\}\):

Cov⁡\(𝐟itot∣𝐮g,𝜹i,𝐮i\)\\displaystyle\\mathrm\{Cov\}\(\\mathbf\{f\}\_\{i\}^\{\\mathrm\{tot\}\}\\mid\\mathbf\{u\}\_\{g\},\\bm\{\\delta\}\_\{i\},\\mathbf\{u\}\_\{i\}\)=Cov⁡\(fg​\(𝐗i\)∣𝐮g\)\+Cov⁡\(fδ,i​\(𝐗i\)∣𝜹i\)\\displaystyle=\\mathrm\{Cov\}\(f\_\{g\}\(\\mathbf\{X\}\_\{i\}\)\\mid\\mathbf\{u\}\_\{g\}\)\+\\mathrm\{Cov\}\(f\_\{\\delta,i\}\(\\mathbf\{X\}\_\{i\}\)\\mid\\bm\{\\delta\}\_\{i\}\)\+Cov⁡\(fi​\(𝐗i\)∣𝐮i\)\\displaystyle\+\\mathrm\{Cov\}\(f\_\{i\}\(\\mathbf\{X\}\_\{i\}\)\\mid\\mathbf\{u\}\_\{i\}\)=\(𝐊g,i​i−𝐐g,i​i\)\+\(𝐊δ,i​i−𝐐δ,i​i\)\+\(𝐊i,i​i−𝐐i,i​i\)\\displaystyle=\(\\mathbf\{K\}\_\{g,ii\}\-\\mathbf\{Q\}\_\{g,ii\}\)\+\(\\mathbf\{K\}\_\{\\delta,ii\}\-\\mathbf\{Q\}\_\{\\delta,ii\}\)\+\(\\mathbf\{K\}\_\{i,ii\}\-\\mathbf\{Q\}\_\{i,ii\}\)=𝐑i\.\\displaystyle=\\mathbf\{R\}\_\{i\}\.

##### Diagonal Residual Structure via Trace Term\.

The VFE objective function \(ELBO\) contains the expected log\-likelihood term𝔼q​\[log⁡p​\(𝐲\|𝐟\)\]\\mathbb\{E\}\_\{q\}\[\\log p\(\\mathbf\{y\}\|\\mathbf\{f\}\)\], which analytically reduces to a form involving−12​σ2​Tr⁡\(𝐊n​n−𝐐n​n\)\-\\frac\{1\}\{2\\sigma^\{2\}\}\\operatorname\{Tr\}\(\\mathbf\{K\}\_\{nn\}\-\\mathbf\{Q\}\_\{nn\}\)\. Since the trace operator sums only the diagonal elements, the optimization objective depends only on the marginal variances of the residual process\([Titsias, 2009](https://arxiv.org/html/2609.19337#bib.bib29)\)\. This allows us to compute the residual term efficiently using only the diagonal elements of\(𝐊⋅,i​i−𝐐⋅,i​i\)\(\\mathbf\{K\}\_\{\\cdot,ii\}\-\\mathbf\{Q\}\_\{\\cdot,ii\}\), without assuming the off\-diagonal elements are zero\. This yields the computationally efficient form of𝐑i\\mathbf\{R\}\_\{i\}mentioned in the proposition\.□\\square

## Appendix FProof of Proposition[7](https://arxiv.org/html/2609.19337#Thmproposition7)\(Global ELBO given fixed local blocks\)

Starting from the explicit ELBO decomposition in Theorem[1](https://arxiv.org/html/2609.19337#Thmtheorem1),

ℒ\\displaystyle\\mathcal\{L\}=∑i=1T\[−ni2log\(2πσ2\)−12​σ2∥𝐲i−𝝁𝐟i∥2\\displaystyle=\\sum\_\{i=1\}^\{T\}\\Bigl\[\-\\frac\{n\_\{i\}\}\{2\}\\log\(2\\pi\\sigma^\{2\}\)\-\\frac\{1\}\{2\\sigma^\{2\}\}\\\|\\mathbf\{y\}\_\{i\}\-\\bm\{\\mu\}\_\{\\mathbf\{f\}\_\{i\}\}\\\|^\{2\}−12​σ2​\(Tr⁡\(𝐀g,i​𝐒g​𝐀g,i⊤\)\+Tr⁡\(𝐀δ,i​𝐒δi​𝐀δ,i⊤\)\+Tr⁡\(𝐀i​𝐒ui​𝐀i⊤\)CLOSE\\displaystyle\\qquad\\qquad\-\\frac\{1\}\{2\\sigma^\{2\}\}\\bigl\(\\operatorname\{Tr\}\(\\mathbf\{A\}\_\{g,i\}\\mathbf\{S\}\_\{g\}\\mathbf\{A\}\_\{g,i\}^\{\\top\}\)\+\\operatorname\{Tr\}\(\\mathbf\{A\}\_\{\\delta,i\}\\mathbf\{S\}\_\{\\delta\_\{i\}\}\\mathbf\{A\}\_\{\\delta,i\}^\{\\top\}\)\+\\operatorname\{Tr\}\(\\mathbf\{A\}\_\{i\}\\mathbf\{S\}\_\{u\_\{i\}\}\\mathbf\{A\}\_\{i\}^\{\\top\}\)\+Tr\(𝐊g,i​i−𝐐g,i​i\)\+Tr\(𝐊δ,i​i−𝐐δ,i​i\)\+Tr\(𝐊i,i​i−𝐐i,i​i\)\)\]\\displaystyle\\qquad\\qquad\\quad\+\\operatorname\{Tr\}\(\\mathbf\{K\}\_\{g,ii\}\-\\mathbf\{Q\}\_\{g,ii\}\)\+\\operatorname\{Tr\}\(\\mathbf\{K\}\_\{\\delta,ii\}\-\\mathbf\{Q\}\_\{\\delta,ii\}\)\+\\operatorname\{Tr\}\(\\mathbf\{K\}\_\{i,ii\}\-\\mathbf\{Q\}\_\{i,ii\}\)\\bigr\)\\Bigr\]−KL\(q\(𝐮g\)∥p\(𝐮g\)\)−∑i=1TKL\(q\(𝜹i\)∥p\(𝜹i\)\)−∑i=1TKL\(q\(𝐮i\)∥p\(𝐮i\)\),\\displaystyle\\quad\-\\mathrm\{KL\}\\bigl\(q\(\\mathbf\{u\}\_\{g\}\)\\\|p\(\\mathbf\{u\}\_\{g\}\)\\bigr\)\-\\sum\_\{i=1\}^\{T\}\\mathrm\{KL\}\\bigl\(q\(\\bm\{\\delta\}\_\{i\}\)\\\|p\(\\bm\{\\delta\}\_\{i\}\)\\bigr\)\-\\sum\_\{i=1\}^\{T\}\\mathrm\{KL\}\\bigl\(q\(\\mathbf\{u\}\_\{i\}\)\\\|p\(\\mathbf\{u\}\_\{i\}\)\\bigr\),
we treat\{Θi,local\(r\+1\)\}\\\{\\Theta\_\{i,\\mathrm\{local\}\}^\{\(r\+1\)\}\\\}as fixed and collect all terms that do not depend onΘglobal\\Theta\_\{\\mathrm\{global\}\}into a constant\. The terms depending onΘglobal\\Theta\_\{\\mathrm\{global\}\}are the least squares term, the variational covariance tracesTr⁡\(𝐀g,i​𝐒g​𝐀g,i⊤\)\\operatorname\{Tr\}\(\\mathbf\{A\}\_\{g,i\}\\mathbf\{S\}\_\{g\}\\mathbf\{A\}\_\{g,i\}^\{\\top\}\)andTr⁡\(𝐀δ,i​𝐒δi​𝐀δ,i⊤\)\\operatorname\{Tr\}\(\\mathbf\{A\}\_\{\\delta,i\}\\mathbf\{S\}\_\{\\delta\_\{i\}\}\\mathbf\{A\}\_\{\\delta,i\}^\{\\top\}\)\(with𝐀δ,i\\mathbf\{A\}\_\{\\delta,i\}depending on𝐙g\\mathbf\{Z\}\_\{g\}andkgk\_\{g\}\), the global and deviation residual traces, the global KL, and the deviation KLs∑iKL\(q\(𝜹i\)∥p\(𝜹i;ϕ,𝐊g​g\)\)\\sum\_\{i\}\\mathrm\{KL\}\(q\(\\bm\{\\delta\}\_\{i\}\)\\\|p\(\\bm\{\\delta\}\_\{i\};\\phi,\\mathbf\{K\}\_\{gg\}\)\)becausep⁡\(𝜹i\)p\(\\bm\{\\delta\}\_\{i\}\)depends on\(ϕ,𝐊g​g\)\(\\phi,\\mathbf\{K\}\_\{gg\}\)\. This yields \([21](https://arxiv.org/html/2609.19337#S3.E21)\)–\([22](https://arxiv.org/html/2609.19337#S3.E22)\) up to an additive constant independent ofΘglobal\\Theta\_\{\\mathrm\{global\}\}\.□\\square

## References

- Achituve et al\. \(2021\)Achituve I, Shamsian A, Navon A, Chechik G, Fetaya E \(2021\) Personalized federated learning with gaussian processes\.*Advances in Neural Information Processing Systems*34:8392–8406\.
- Bruinsma et al\. \(2020\)Bruinsma W, Perim E, Tebbutt W, Hosking S, Solin A, Turner R \(2020\) Scalable exact inference in multi\-output gaussian processes\.*International Conference on Machine Learning*, 1190–1201 \(PMLR\)\.
- Chakraborty et al\. \(2020\)Chakraborty S, Adhikari S, Ganguli R \(2020\) The role of surrogate models in the development of digital twins of dynamic systems\.*arXiv preprint arXiv:2001\.09292*\.
- Chellapandi et al\. \(2023\)Chellapandi VP, Yuan L, Brinton CG, Żak SH, Wang Z \(2023\) Federated learning for connected and automated vehicles: A survey of existing approaches and challenges\.*IEEE Transactions on Intelligent Vehicles*9\(1\):119–137\.
- Chung and Al Kontar \(2024\)Chung S, Al Kontar R \(2024\) Federated multi\-output gaussian processes\.*Technometrics*66\(1\):90–103\.
- Deisenroth and Ng \(2015\)Deisenroth M, Ng JW \(2015\) Distributed gaussian processes\.*International conference on machine learning*, 1481–1490 \(PMLR\)\.
- Deisenroth et al\. \(2013\)Deisenroth MP, Fox D, Rasmussen CE \(2013\) Gaussian processes for data\-efficient learning in robotics and control\.*IEEE transactions on pattern analysis and machine intelligence*37\(2\):408–423\.
- Deng et al\. \(2020\)Deng Y, Kamani MM, Mahdavi M \(2020\) Adaptive personalized federated learning\.*arXiv preprint arXiv:2003\.13461*\.
- Gal et al\. \(2014\)Gal Y, Van Der Wilk M, Rasmussen CE \(2014\) Distributed variational inference in sparse gaussian process regression and latent variable models\.*Advances in neural information processing systems*27\.
- Gao and Chung \(2024\)Gao J, Chung S \(2024\) Federated automatic latent variable selection in multi\-output gaussian processes\.*arXiv preprint arXiv:2407\.16935*\.
- Guo et al\. \(2022\)Guo Y, Hu T, Zhou Y, Zhao K, Zhang Z \(2022\) Multi\-channel data fusion and intelligent fault diagnosis based on deep learning\.*Measurement Science and Technology*34\(1\):015115\.
- Hu et al\. \(2025\)Hu J, Shi N, Al Kontar R, Yan H \(2025\) Personalized tucker decomposition: Modeling commonality and peculiarity on tensor data\.*Technometrics*1–17\.
- Jeong et al\. \(2025\)Jeong C, Yue X, Chung S \(2025\) Fed\-joint: Joint modeling of nonlinear degradation signals and failure events for remaining useful life prediction using federated learning\.*arXiv preprint arXiv:2503\.13404*\.
- Kennedy and O’Hagan \(2001\)Kennedy MC, O’Hagan A \(2001\) Bayesian calibration of computer models\.*Journal of the Royal Statistical Society: Series B \(Statistical Methodology\)*63\(3\):425–464\.
- Kessels et al\. \(2022\)Kessels B, Korver J, Fey R, van de Wouw N \(2022\) Model updating for digital twins using gaussian process inverse mapping models\.*ENOC 2020\+2*\.
- Kontar et al\. \(2021\)Kontar R, Shi N, Yue X, Chung S, Byon E, Chowdhury M, Jin J, Kontar W, Masoud N, Nouiehed M, et al\. \(2021\) The internet of federated things \(ioft\)\.*IEEE Access*9:156071–156113\.
- Kontoudis and Stilwell \(2022\)Kontoudis GP, Stilwell DJ \(2022\) Fully decentralized, scalable gaussian processes for multi\-agent federated learning\.*CoRR*\.
- LEI et al\. \(2010\)LEI Y, ZHANG Z, JIN J \(2010\) Automatic tonnage monitoring for missing part detection in multi\-operation forging processes\.*Journal of manufacturing science and engineering*132\(5\)\.
- Liu et al\. \(2018\)Liu H, Cai J, Wang Y, Ong YS \(2018\) Generalized robust bayesian committee machine for large\-scale gaussian process regression\.*International Conference on Machine Learning*, 3131–3140 \(PMLR\)\.
- Llorente et al\. \(2025\)Llorente F, Waxman D, Jantre S, Urban NM, Minkoff SE \(2025\) Robust, online, and adaptive decentralized gaussian processes\.*arXiv preprint arXiv:2509\.18011*\.
- Matthews et al\. \(2017\)Matthews AGdG, van der Wilk M, Nickson T, Fujii K, Boukouvalas A, León\-Villagrá P, Ghahramani Z, Hensman J \(2017\) Scalable gaussian process inference using variational methods\.*Advances in Neural Information Processing Systems*, volume 30\.
- Ng and Deisenroth \(2014\)Ng JW, Deisenroth MP \(2014\) Hierarchical mixture\-of\-experts model for large\-scale gaussian process regression\.*arXiv preprint arXiv:1412\.3078*\.
- Nguyen and Bonilla \(2014\)Nguyen T, Bonilla E \(2014\) Fast allocation of gaussian process experts\.*International conference on machine learning*, 145–153 \(PMLR\)\.
- Peng et al\. \(2017\)Peng H, Zhe S, Zhang X, Qi Y \(2017\) Asynchronous distributed variational gaussian process for regression\.*International Conference on Machine Learning*, 2788–2797 \(PMLR\)\.
- Pruckovskaja et al\. \(2023\)Pruckovskaja V, Weissenfeld A, Heistracher C, Graser A, Kafka J, Leputsch P, Schall D, Kemnitz J \(2023\) Federated learning for predictive maintenance and quality inspection in industrial applications\.*2023 Prognostics and Health Management Conference \(PHM\)*, 312–317 \(IEEE\)\.
- Seeger \(2004\)Seeger M \(2004\) Gaussian processes for machine learning\.*International journal of neural systems*14\(02\):69–106\.
- Shi and Al Kontar \(2024\)Shi N, Al Kontar R \(2024\) Personalized pca: Decoupling shared and unique features\.*Journal of machine learning research*25\(41\):1–82\.
- Thorgeirsson and Gauterin \(2020\)Thorgeirsson AT, Gauterin F \(2020\) Probabilistic predictions with federated learning\.*Entropy*23\(1\):41\.
- Titsias \(2009\)Titsias MK \(2009\) Variational learning of inducing variables in sparse gaussian processes\.*International Conference on Artificial Intelligence and Statistics*, 567–574 \(PMLR\)\.
- Tresp \(2000\)Tresp V \(2000\) A bayesian committee machine\.*Neural computation*12\(11\):2719–2741\.
- Wilson and Adams \(2013\)Wilson A, Adams R \(2013\) Gaussian process kernels for pattern discovery and extrapolation\.*International conference on machine learning*, 1067–1075 \(PMLR\)\.
- Xu et al\. \(2021\)Xu J, Glicksberg BS, Su C, Walker P, Bian J, Wang F \(2021\) Federated learning for healthcare informatics\.*Journal of healthcare informatics research*5:1–19\.
- Xu et al\. \(2019\)Xu Y, Yin F, Xu W, Lin J, Cui S \(2019\) Distributed gaussian process: New paradigm and application to wireless traffic prediction\.*ICC 2019\-2019 IEEE International Conference on Communications \(ICC\)*, 1–6 \(IEEE\)\.
- Yang et al\. \(2024\)Yang H, Li J, Hao M, Zhang W, He H, Sangaiah AK \(2024\) An efficient personalized federated learning approach in heterogeneous environments: a reinforcement learning perspective\.*Scientific Reports*14\(1\):28877\.
- Ye et al\. \(2019\)Ye F, Zhang Z, Xia Z, Zhou Y, Zhang H \(2019\) Monitoring and diagnosis of multi\-channel profile data based on uncorrelated multilinear discriminant analysis\.*The International Journal of Advanced Manufacturing Technology*103\(9\):4659–4669\.
- Ye et al\. \(2023\)Ye M, Fang X, Du B, Yuen PC, Tao D \(2023\) Heterogeneous federated learning: State\-of\-the\-art and research challenges\.*ACM Computing Surveys*56\(3\):1–44\.
- Yin and Gunnarsson \(2017\)Yin F, Gunnarsson F \(2017\) Distributed recursive gaussian processes for rss map applied to target tracking\.*IEEE Journal of Selected Topics in Signal Processing*11\(3\):492–503\.
- Yin et al\. \(2020\)Yin F, Lin Z, Kong Q, Xu Y, Li D, Theodoridis S, Cui SR \(2020\) Fedloc: Federated learning framework for data\-driven cooperative localization and location data processing\.*IEEE Open Journal of Signal Processing*1:187–215\.
- Yu et al\. \(2022\)Yu H, Guo K, Karami M, Chen X, Zhang G, Poupart P \(2022\) Federated bayesian neural regression: A scalable global federated gaussian process\.*arXiv preprint arXiv:2206\.06357*\.
- Yue et al\. \(2024\)Yue X, Kontar RA, Gomez AME \(2024\) Federated data analytics: A study on linear models\.*IISE Transactions*56\(1\):16–28\.
- Yurochkin et al\. \(2019\)Yurochkin M, Agarwal M, Ghosh S, Greenewald K, Hoang N, Khazaeni Y \(2019\) Bayesian nonparametric federated learning of neural networks\.*International conference on machine learning*, 7252–7261 \(PMLR\)\.
- Zhou et al\. \(2015\)Zhou C, Liu K, Zhang X, Zhang W, Shi J \(2015\) An automatic process monitoring method using recurrence plot in progressive stamping processes\.*IEEE Transactions on Automation Science and Engineering*13\(2\):1102–1111\.