weather data. We investigate four scenarios: PV-S (150), PV-S (610), PV-EW
(150), and PV-EW (610). In the PV-S variants, all houses face south, while the
PV-EW variants involve houses oriented either east or west. The PV-EW
scenarios, being bimodal and highly heterogeneous, present greater modeling
challenges. To assess the impact of dataset sizes (denoted as $\mathbf{m}$),
we consider two settings: clients with either $150$ or $610$ training samples,
corresponding to two-week and two-month data windows, respectively. The number
in each dataset’s name indicates the training sample size.

##### Data generation.

To generate simulated datasets for the PV-S and PV-EW experiments, we first
obtained a real dataset of hourly solar radiation and meteorological
measurements in Lausanne from the Photovoltaic Geographical Information System
(PVGIS) online tool [PVG, ]. We then used the pvlib Python library [Holmgren
et al., 2022] to simulate the PV power output based on these measurements. To
simulate the PV panels in different houses, we sampled meteorological data and
installation specifications from normal distributions with specified means and
standard deviations. For example, in the PV-EW experiment, the azimuths of the
PV panels were sampled from a bimodal normal distribution of
$0.5\big{(}\mathcal{N}(90^{\circ},15^{\circ})+\mathcal{N}(270^{\circ},15^{\circ})\big{)}$,
where $\mathcal{N}$ denotes the normal distribution. Fig. 4 illustrates the
power output profiles of 24 houses in the PV-EW experiment over five days.
This figure demonstrates that the curves have noticeable differences; however,
they show similar trends. As expected, the houses facing east or west are more
similar to one another, while the differences among different sub-populations
are higher. For instance, the peak production of the houses facing the east
occurs in the morning, whereas the peak for those facing the west occurs in
the afternoon.

Figure 4:  Power output profile of $24$ houses in the PV-EW experiment over
five days in June $2018$, where each line represents the PV generation of one
house. Green and blue curves correspond to houses facing the east and the west
respectively. Although the curves have noticeable differences, there are
consistent trends present in the data.

To ensure that each client’s dataset is identically distributed, despite the
variations in meteorological patterns throughout the year, we filter the
dataset down to the months of June and July. We employ two distinct training
datasets for our analysis. The first dataset comprises the initial two weeks
of June $2018$, which provides a total of $150$ samples for each client. The
second dataset encompasses the data from both June and July $2018$, resulting
in $610$ training samples per client. For all experiments, the test dataset
consists of the data from June and July $2019$. Furthermore, we exclude
nighttime data points when the PV generation is zero. To normalize the data,
we standardize each house’s features and output, setting the mean to zero and
the standard deviation to one.

##### Features.

The dataset contains various meteorological features, including solar beam and
diffuse irradiances, temperature, wind speed, solar altitude, time, and date.
Solar beam and diffuse irradiances are particularly valuable in predicting PV
generation, but their measurement can be costly. To reflect a realistic
scenario, we assume that all households have access to the irradiance data
recorded at a specific location in the city, such as a weather station.
However, they are unaware of the specific irradiance values at their own
houses, which may differ from the weather station. The temperature and wind
speed may vary slightly among different houses in the dataset. Additionally,
some houses may experience intermittent shadows caused by nearby trees or
buildings, which can appear and disappear at certain times of the day. These
shadows are considered as noise, and no recorded feature in the dataset
provides explicit information about their occurrence or characteristics.

In addition to the meteorological features, we also incorporate autoregressors
in our analysis. Autoregressors are advantageous in time-series prediction
tasks and aim to capture the temporal dependencies in the data. To determine
the autoregressors to include in our analysis, we utilize the partial
autocorrelation function (PACF). The PACF measures the correlation between a
time series and its lagged values while controlling for the correlations with
all shorter lags. We select autoregressors with the highest values of the PACF
among the lagged values up to two weeks ago. By focusing on these highly
correlated autoregressors, we aim to capture the most relevant information
from the past time steps.

####  8.5.2 Polynomial dataset

We examine a bimodal dataset where the data for each client is generated by
sampling a function from one of two GP priors. Each GP prior is characterized
by a polynomial mean function of order $7$ and an SE kernel function. The two
modes have distinct polynomial means and length scales associated with the SE
kernel. Additionally, Gaussian noise is introduced into the generated data to
account for measurement errors and other sources of variability. Fig. 5
illustrates the dataset for a total of $24$ clients, with $12$ clients
belonging to the first mode and another $12$ clients belonging to the second
mode, where each client has $10$ training samples.

Figure 5:  Polynomial dataset for $24$ clients, where $12$ clients belong to
the first and the rest belong to the second mode. The solid lines represent
the unknown true data-generating distribution for each client. Each client has
a training dataset with $10$ samples. The solid lines in the figure represent
the true data-generating distribution, which is unknown to us. The points on
the graph represent the noisy samples in each client’s training dataset.

####  8.5.3 FEMNIST and EMNIST Datasets

The FEMNIST (Federated Extended MNIST) [Caldas et al., 2019] and EMNIST
(Extended MNIST) [Cohen et al., 2017] datasets consist of $28\times 28$ images
of hand-written English characters. Both datasets are used for handwritten
character recognition. Specifically, the FEMNIST dataset is employed in a
$10$-way classification task, focusing on the classification of digits.
Meanwhile, the EMNIST dataset is utilized for a more comprehensive $62$-way
classification task. This broader task encompasses the recognition of both
lower case and upper case English alphabet letters, in addition to digits.

FEMNIST identifies the writer for each character, allowing a natural
partitioning strategy by assigning all images from the same writer to a single
client [Caldas et al., 2019]. This dataset is heterogeneous due to the
inherent diversity in handwriting styles, resulting in skewness within the
feature distribution [Tan et al., 2021]. For our analysis, we subsample $40$
clients and explore two scenarios: FEMNIST ($20$), where each client is
allocated $20$ samples and FEMNIST ($500$), where each client utilizes all
available samples, approximately $500$ per client.

Following [Wang et al., 2020, Marfoq et al., 2021], we subsample $10\%$ of the
EMNIST dataset that amounts to $81425$ total samples. We distribute samples
with the same label across $80$ clients according to a symmetric Dirichlet
distribution with parameter $\alpha=0.4$. This approach introduces a skew in
the label distribution [Tan et al., 2021], thereby synthesizing heterogeneity
in the dataset.

###  8.6 Experiments details

####  8.6.1 Calibration error

Calibration error (CE) is an evaluation metric relevant to probabilistic
models. It is based on the premise that when provided with a test input
$\mathbf{x}_{j}$, the model produces a probability distribution
$\hat{p}(y_{j}|\mathbf{x}_{j})$ over the predicted target $y_{j}$. CE measures
the disparity between the predicted confidence intervals and the actual
proportions of test data falling within those intervals [Kuleshov et al.,
2018]. CE is a nonnegative measure and a smaller CE value is more desirable.

##### Calibration error for regression.

In defining CE for regression, we follow the formula proposed by [Rothfuss et
al., 2021]. Denote a predictor’s cumulative density function (CDF) as
$\hat{F}(y_{j}|\mathbf{x}_{j})=\int_{-\infty}^{y_{j}}\hat{p}(y|\mathbf{x}_{j})d_{y}$.
Given a dataset
$\mathcal{S}=\bigl{\\{}(\mathbf{x}_{j},y_{j})\bigr{\\}}_{j=1}^{m}$ with $m$
samples, we compute the corresponding empirical frequency for confidence
levels $0\leq q_{1}\leq\cdots\leq q_{H}\leq 1$ as follows:

$\hat{q}_{h}\coloneqq\frac{1}{m}\Bigl{|}\bigl{\\{}y_{j}\bigl{|}\hat{F}(y_{j}|\mathbf{x}_{j})\leq
q_{h},\quad j=1,\cdots,m\bigr{\\}}\Bigr{|},$

for $h=1,\cdots,H$.

If the predictions are well-calibrated, we expect that
$\hat{q}_{h}\xrightarrow[]{}q_{h}$ as $m\xrightarrow[]{}\infty$. We adopt the
CE definition proposed by Rothfuss et al. [2021], which formulates CE as a
function of residuals $\hat{q}_{h}-q_{h}$:

$CE\coloneqq\frac{1}{H}\sum_{h=1}^{H}|\hat{q}_{h}-q_{h}|.$ (40)

In our experiments, we evaluate (40) by employing $H=20$ equally spaced
confidence levels ranging from $0$ to $1$.

##### Calibration error for classification.

Given a test input $\mathbf{x}_{j}$, a probabilistic classifier outputs a
categorical probability distribution, $\hat{p}(y_{j}=k|\mathbf{x}_{j})$, where
$k=1,\cdots,C$, with $C$ representing the number of classes. The classifier’s
prediction is the class label with the highest probability, denoted as
$\hat{y}_{j}=\operatorname*{arg\,max}_{k}\hat{p}(y_{j}=k|\mathbf{x}_{j})$.
Correspondingly, we define the classifier’s confidence in the prediction for
the input $\mathbf{x}_{j}$ as
$\hat{p}_{j}\coloneqq\hat{p}(y_{j}=\hat{y}_{j}|\mathbf{x}_{j})$.

For a well-calibrated classifier, the confidence is expected to align with the
probability of correct classification. For example, with $100$ predictions,
each having a confidence of $0.8$, one would anticipate approximately $80$
correct classifications [Guo et al., 2017]. In an empirical comparison of
calibration and accuracy on a dataset
$\mathcal{S}=\bigl{\\{}(\mathbf{x}_{j},y_{j})\bigr{\\}}_{j=1}^{m}$, following
the methodology of Guo et al. [2017], data samples are grouped into $H=20$
intervals, each of length $1/H$, based on their prediction confidence. Let
$B_{h}\coloneqq\bigl{\\{}j|\hat{p}_{j}\in(\frac{h-1}{H},\frac{h}{H}]\bigr{\\}}$
for $h=1,\cdots,H$ be the set of indices of points in $\mathcal{S}$ whose
prediction confidence falls within interval $h$. The accuracy and average
confidence of points in $B_{h}$ are defined as follows:

$acc(B_{h})\coloneqq\frac{1}{|B_{h}|}\sum_{j\in
B_{h}}\mathbf{1}(\hat{y}_{j}=y_{j}),\\\
conf(B_{h})\coloneqq\frac{1}{|B_{h}|}\sum_{j\in B_{h}}\hat{p}_{j}.$

To compare the accuracy and average confidence across all intervals, Guo et
al. [2017] calculates CE as follows:

$CE\coloneqq\sum_{h=1}^{H}\frac{|B_{h}|}{m}|acc(B_{h})-conf(B_{h})|.$ (41)

This formula measures the deviation between accuracy and prediction confidence
level, serving as a metric for assessing the calibration of classification
tasks.

####  8.6.2 Regression experiments details

##### Baselines details.

In the Vanilla approach, the GP hyperparameters for each client are tuned by
maximizing the LML of that specific client, without using FL. For pooled GP,
we use inducing points to handle computational issues due to the large number
of samples involved in this approach. pFedGP optimizes the average LML across
clients to obtain the deep kernel hyper-parameters. We adapted pFedGP,
originally designed specifically for classification, to suit our regression
task. With this adaptation, pFedGP is a special case of our algorithm when
using a single prior in SVGD and a very wide hyper-prior. We try pFedGP with a
zero GP mean, as originally done in Achituve et al. [2021], and a NN mean and
do not use inducing points. We apply pFedBayes [Zhang et al., 2022] only to
the FEMNIST data as the authors only consider classification tasks in their
experiments.

##### Hyper-parameter tuning.

We perform hyper-parameter tuning for all baselines and our method using
$5$-fold cross-validation. For all neural networks, we explore structures with
the same number of neurons per layer. The number of neurons per layer can take
values of $2^{n}$ for $n\in{1,\cdots,6}$, and we consider $2$ or $4$ hidden
layers. For PAC-PFL, we employ $4$ SVGD particles and set $k=4$. The parameter
$\beta$ is set equal to the number of samples for each client, $\beta=m_{i}$.
To determine the value of $\tau$, we conduct a search for values greater than
$1/(1+n\beta)$, ensuring that $\lambda>n_{2}+\upsilon$ and thus, satisfying
the assumption outlined in Theorem 4.1.

The employed hyper-prior is a multivariate Gaussian distribution with a
diagonal covariance matrix. It is a distribution over the weights and biases
of the mean and kernel neural networks, as well as the noise standard
deviation of the likelihood. In all PV experiments, we set the hyper-prior
mean for the neural network weights and biases to $0$ and the hyper-prior mean
for the noise standard deviation to $0.4$. These choices help prevent
overfitting.

##### Results.

In this section, we evaluate our algorithm and compare it to the baseline
methods introduced earlier. We assess the performance based on the RSMSE and
CE error metrics for both existing and new clients. To obtain these metrics
for a given method, we calculate the RSMSE and CE, where applicable, for each
existing and new client using a fixed random seed. Then, we compute $4$ values
for the average RSMSE and CE metrics across the existing and new client
groups. We performed the analysis five times using different random seeds and
calculated the $95\%$ confidence interval for the sample mean of RSMSE and CE
over both the existing and new clients.

The results are presented in Tables 4 to 7. In each column, the number after
the dataset name represents the number of training samples per client. As
shown in Tables 4 and 5, our method outperforms other approaches in the
majority of cases in terms of RSMSE for both existing and new clients.
Furthermore, we observe that our method’s performance improves with an
increase in the number of training samples per client.

Concerning CE, reported in Tables 6 and 7, PAC-PFL consistently exhibits low
CE values across all scenarios, establishing itself as the top-performing or
closely competitive method. This observation aligns with findings in Rothfuss
et al. [2021], where the authors noted that CE improvement is notable when the
meta-learning tasks exhibit greater similarity. In our own experiments, given
the high client heterogeneity present in both PV-EW and PV-S, the improvements
in CE are comparatively modest.

We believe there are three reasons why PAC-PFL outperforms other baselines.
First, by learning the prior with FL and treating posterior inference as
personalization, PAC-PFL exhibits a high level of adaptability to individual
patterns. Second, PAC-PFL’s ability to learn multiple priors enables it to
effectively model clients with heterogeneous data distributions. Finally, the
regularization of the hyper-posterior towards the hyper-prior in PAC-PFL helps
prevent overfitting and allows for learning more complex prior means even with
limited training data. For example, when selecting the GP prior mean structure
through cross-validation, it is observed that the best GP mean for pFedGP is a
linear function, while for PAC-PFL it is a 2-layer neural network with $32$
neurons per layer. The use of a more complex prior mean in PAC-PFL allows for
the potential capture of more intricate patterns in the data.

Table 4:  RSMSE of existing clients for various methods. We report $95\%$ confidence interval for the average test RSMSE across all existing clients. Our proposed algorithm, PAC-PFL, demonstrates superior performance compared to all other baselines across all datasets. Method | PV-S (150) | PV-EW (150) | PV-S (610) | PV-EW (610) | Polynomial (10)
---|---|---|---|---|---
PAC-PFL | 0.429 $\pm$ 0.016 | 0.412 $\pm$ 0.008 | 0.425 $\pm$ 0.014 | 0.403 $\pm$ 0.008 | 0.582 $\pm$ 0.048
pFedGP | 0.479 $\pm$ 0.038 | 0.452 $\pm$ 0.013 | 0.514 $\pm$ 0.034 | 0.457 $\pm$ 0.012 | 0.802 $\pm$ 0.147
MTL | 0.445 $\pm$ 0.000 | 0.426 $\pm$ 0.001 | 0.441 $\pm$ 0.005 | 0.414 $\pm$ 0.000 | 0.851 $\pm$ 0.150
MAML | 0.568 $\pm$ 0.023 | 0.523 $\pm$ 0.027 | 0.580 $\pm$ 0.016 | 0.524 $\pm$ 0.015 | 0.976 $\pm$ 0.129
Vanilla | 0.683 $\pm$ 0.024 | 0.626 $\pm$ 0.033 | 0.488 $\pm$ 0.044 | 0.458 $\pm$ 0.015 | 0.727 $\pm$ 0.067
Pooled | 0.484 $\pm$ 0.034 | 0.467 $\pm$ 0.045 | 0.492 $\pm$ 0.022 | 0.461 $\pm$ 0.030 | 2.328 $\pm$ 0.545
Table 5: RSMSE of new clients for various methods. We report $95\%$ confidence interval for the average test RSMSE across all new clients. Our proposed algorithm, PAC-PFL, outperforms other baselines in the majority of cases. The column name indicates the number of training samples specifically for the existing clients. However, the number of training samples for new clients remains constant at $150$. In the Vanilla case, where each client trains their model individually, the size of the training data for the existing clients does not impact the performance for the new clients. Hence, the RSMSE values for the new clients employing the Vanilla method is the same in both the PV-S (150) and PV-S (610) datasets (similarly for PV-EW (150) and PV-EW (610)). Method | PV-S (150) | PV-EW (150) | PV-S (610) | PV-EW (610) | Polynomial (10)
---|---|---|---|---|---
PAC-PFL | 0.432 $\pm$ 0.018 | 0.417 $\pm$ 0.011 | 0.433 $\pm$ 0.012 | 0.430 $\pm$ 0.010 | 0.654 $\pm$ 0.024
pFedGP | 0.464 $\pm$ 0.042 | 0.447 $\pm$ 0.017 | 0.513 $\pm$ 0.053 | 0.457 $\pm$ 0.012 | 0.874 $\pm$ 0.408
MTL | 0.440 $\pm$ 0.001 | 0.423 $\pm$ 0.001 | 0.449 $\pm$ 0.002 | 0.430 $\pm$ 0.000 | 0.885 $\pm$ 0.123
MAML | 0.546 $\pm$ 0.029 | 0.522 $\pm$ 0.013 | 0.552 $\pm$ 0.014 | 0.530 $\pm$ 0.021 | 1.107 $\pm$ 0.284
Vanilla | 0.687 $\pm$ 0.020 | 0.632 $\pm$ 0.053 | 0.687 $\pm$ 0.020 | 0.632 $\pm$ 0.053 | 0.763 $\pm$ 0.046
Pooled | 0.477 $\pm$ 0.029 | 0.469 $\pm$ 0.046 | 0.481 $\pm$ 0.022 | 0.460 $\pm$ 0.028 | 2.321 $\pm$ 0.371
Table 6: CE of existing clients for various methods. We report $95\%$ confidence interval for the average test CE across all existing clients. Method | PV-S (150) | PV-EW (150) | PV-S (610) | PV-EW (610) | Polynomial (10)
---|---|---|---|---|---
PAC-PFL | 0.072 $\pm$ 0.003 | 0.041 $\pm$ 0.002 | 0.068 $\pm$ 0.002 | 0.044 $\pm$ 0.003 | 0.140 $\pm$ 0.017
pFedGP | 0.067 $\pm$ 0.020 | 0.057 $\pm$ 0.006 | 0.066 $\pm$ 0.002 | 0.052 $\pm$ 0.007 | 0.119 $\pm$ 0.012
Vanilla | 0.122 $\pm$ 0.008 | 0.120 $\pm$ 0.011 | 0.036 $\pm$ 0.002 | 0.028 $\pm$ 0.001 | 0.178 $\pm$ 0.011
Pooled | 0.256 $\pm$ 0.004 | 0.260 $\pm$ 0.010 | 0.263 $\pm$ 0.005 | 0.266 $\pm$ 0.010 | 0.326 $\pm$ 0.036
Table 7: CE of new clients for various methods. We report the $95\%$ confidence interval for the average test CE across all new clients. The CE values for the new clients employing the Vanilla method are not affected by the number of training samples for the existing clients, as indicated in the column names. Method | PV-S (150) | PV-EW (150) | PV-S (610) | PV-EW (610) | Polynomial (10)
---|---|---|---|---|---
PAC-PFL | 0.073 $\pm$ 0.005 | 0.044 $\pm$ 0.003 | 0.071 $\pm$ 0.004 | 0.044 $\pm$ 0.002 | 0.165 $\pm$ 0.008
pFedGP | 0.069 $\pm$ 0.015 | 0.061 $\pm$ 0.007 | 0.064 $\pm$ 0.007 | 0.056 $\pm$ 0.008 | 0.147 $\pm$ 0.023
Vanilla | 0.121 $\pm$ 0.011 | 0.120 $\pm$ 0.010 | 0.121 $\pm$ 0.011 | 0.120 $\pm$ 0.010 | 0.237 $\pm$ 0.033
Pooled | 0.262 $\pm$ 0.002 | 0.263 $\pm$ 0.010 | 0.269 $\pm$ 0.002 | 0.266 $\pm$ 0.014 | 0.356 $\pm$ 0.023

####  8.6.3 Classification experiments details

##### Baselines details.

In our experiments, we utilize the original software by Achituve et al. [2021]
for pFedGP and do not use inducing points. Since the source code for pFedBayes
has not been publicly released, we rely on an unofficial
implementation222https://github.com/AllenBeau/pFedBayes. Our implementations
of FedAvg, MTL, and MAML are based on the code provided in Xie et al. [2023].
For the Vanilla model, we follow the procedure described for the PV dataset.
We utilize a Bayesian Convolutional Neural Network (BCNN) [Gal and Ghahramani,
2016] for all Bayesian approaches and a Convolutional Neural Network (CNN) for
the frequentist methods.

##### Hyper-parameter tuning.

Hyper-parameter tuning for all baselines and our method is conducted using the
cross-validation approach. In the case of PAC-PFL, we select $\beta$ and
$\tau$ following the same procedure as employed for the PV dataset. We use $3$
SVGD particles and a multivariate Gaussian hyper-prior with a diagonal
covariance matrix.

##### FEMNIST results.

We evaluate PAC-PFL in two scenarios: first, with 40 clients, each having 20
training samples, and second, with the same 40 clients, each utilizing their
entire available training dataset, averaging around 500 samples per client. In
Table 8, we present the average accuracy and CE scores, together with $95\%$
confidence intervals, across five different random seeds for the existing
clients. Remarkably, in the low-data scenario with only $20$ samples per
client, PAC-PFL consistently outperforms all other methods.

In the data-rich scenario, pFedGP demonstrates a slight advantage over PAC-PFL
in terms of the mean, albeit with a considerably higher standard deviation,
indicating sensitivity to initialization. Given the closely aligned means for
both pFedGP and PAC-PFL and the fact that the confidence interval of PAC-PFL
is encompassed within that of pFedGP, PAC-PFL is a more reliable method. It
offers nearly identical mean performance. Additionally, pFedGP has a
substantially higher computational cost than PACA-PFL. It is worth noting that
the authors of Achituve et al. [2021] propose alternative variants of pFedGP
that trade off accuracy for reduced computational cost. However, for our
experiments, we employ the original algorithm.

Table 8:  Accuracy and calibration error on the FEMNIST classification dataset. We evaluate the metrics in two scenarios: first, with $40$ clients, each having $20$ training samples, and second, with the same $40$ clients, each utilizing their entire available training dataset, averaging around $500$ samples per client. The number after the metric name signifies the per-client training sample size. For each metric, the sample mean and $95\%$ confidence interval using $5$ random seeds is reported. While it is theoretically possible to calculate CE for pFedBayes, it is not included in the table as the available software only provides prediction means and lacks variances required for CE computation.  Method | Accuracy ($20$) | CE ($20$) | Accuracy ($\approx 500$) | CE ($\approx 500$)
---|---|---|---|---
PAC-PFL | 0.942 $\pm$ 0.026 | 0.080 $\pm$ 0.013 | 0.971 $\pm$ 0.016 | 0.046 $\pm$ 0.007
pFedGP | 0.836 $\pm$ 0.020 | 0.096 $\pm$ 0.039 | 0.977 $\pm$ 0.066 | 0.022 $\pm$ 0.008
pFedBayes | 0.870 $\pm$ 0.020 | - | 0.883 $\pm$ 0.024 | -
FedAvg | 0.881 $\pm$ 0.017 | - | 0.967 $\pm$ 0.005 | -
MTL | 0.758 $\pm$ 0.037 | - | 0.875 $\pm$ 0.001 | -
MAML | 0.820 $\pm$ 0.067 | - | 0.881 $\pm$ 0.026 | -
Vanilla | 0.819 $\pm$ 0.013 | 0.120 $\pm$ 0.014 | 0.925 $\pm$ 0.018 | 0.325 $\pm$ 0.062
Pooled | 0.894 $\pm$ 0.043 | 0.106 $\pm$ 0.048 | 0.943 $\pm$ 0.021 | 0.073 $\pm$ 0.033