# Mortality Prediction using Survival Energy Models with Functional Data
Analysis

Daiki<EMAIL_ADDRESS>
Graduate School of Fundamental Science and Engineering, Waseda University
Yasutaka<EMAIL_ADDRESS>
Department of Applied Mathematics, Waseda University

###### Abstract

The Survival Energy Model (SEM), as originally introduced by Shimizu et al.
[11], is designed to characterize human bioenergetics by employing diffusion
processes or inverse Gaussian processes. While parametric models have been
employed to articulate the SEM, they exhibit inherent sensitivity in their
parameters and hyperparameters, which in turn introduces issues of instability
in estimation and prediction. In this paper, we demonstrate that the
utilization of functional data analysis techniques for nonparametric
estimation and prediction of critical functions within the SEM leads to a
substantial enhancement in prediction performance.

Keywords: Survival energy models, ID-SEM, IG-SEM, functional principal
component analysis

MSC2020: 62M20; 62P05; 91D20.

##  1 Introduction

Survival Energy Model (SEM) is a stochastic model proposed by Shimizu et al.
[11] for cohort-wise prediction of mortality. It assumes that human death is
caused by the disappearance of hypothetical ‘survival energy’ and expresses
the dynamics as a stochastic process.

On a stochastic basis $(\Omega,\mathcal{F},{\mathbb{P}};\mathbb{F})$ with the
right-continuous filtration $\mathbb{F}:=(\mathcal{F}_{t})_{t\geq 0}$, a
survival energy process of a certain cohort, say $c$, is defined as a càdlàg
process $X^{c}=(X^{c}_{t})_{t\geq 0}$. Then the time of death, say $\tau$, is
defined as the first-hitting time of $X$ to zero:

$\tau_{c}:=\inf\\{t>0\,|\,X^{c}_{t}<0\\},$

which is a $\mathcal{F}_{t}$-stopping time. Thus, the probability of death up
to time $t>0$ of the corresponding cohort is given by

$\displaystyle q_{c}(t):={\mathbb{P}}(\tau_{c}\leq t),$ (1.1)

which is called the mortality function at cohort $c$.

The first SEM by Shimizu et al. [11] is an inhomogeneous diffusion process,
called ID-SEM:

$\displaystyle
X^{c}_{t}=x_{c}+\int_{0}^{t}U_{c}(s)\,\mathrm{d}s+\int_{0}^{t}V_{c}(s)\,\mathrm{d}W_{s},$
(1.2)

where $x_{c}>0$ is the ‘initial energy’,
$U_{c},\,V_{c}:[0,\infty)\to\mathbb{R}$ are deterministic functions, and $W$
is a Wiener process. The model features intuitive modeling by interpreting
$U_{c}(t)\mathrm{d}t$ as representing the average drift at time $t$, and
$V_{c}^{2}(t)\mathrm{d}t$ as representing its dispersion, which can cause
sudden death if it values larger. Note that this SEM has a continuous path,
and $X^{c}_{\tau_{c}}=0$ almost surely.

On the other hand, the second SEM by Shimizu et al. [12] is based on an
inverse-Gaussian process, say $Y^{c}=(Y^{c}_{t})_{t\geq 0}$:

$\displaystyle X^{c}_{t}=x_{c}-Y^{c}_{t},$ (1.3)

where the probability density of $Y^{c}_{t}$ is given by

$\displaystyle
f_{Y_{t}^{c}}(y)={\color[rgb]{0,0,0}\sqrt{\frac{\sigma_{c}\Lambda_{c}^{2}(t)}{2\pi
y^{3}}}\exp\left(-\frac{\sigma_{c}(y-\Lambda_{c}(t))^{2}}{2y}\right),}$

where $\sigma_{c}>0$ is a constant and $\Lambda_{c}$ is a nonnegative
increasing function with $\Lambda_{c}(0)=0$. This SEM is called IG-SEM, and
the path may have jumps such that $X^{c}_{\tau_{c}}<0$. In particular, if
$\Lambda_{c}(t)=t$, then $X^{c}$ is a spectrally negative Lévy process with
infinite activity jumps. Although the interpretation of $\Lambda_{c}$ is not
as intuitively easy as in the case of ID-SEM, $\Lambda_{c}(t)$ approximately
represents the degree of energy change at time $t$, and the larger
$\Lambda_{c}$, the higher the mortality rate in the model.

The utilization of the SEM offers numerous computational advantages in
actuarial science, particularly in the context of cohort-specific premium
calculations, sensitivity analyses for premiums with respect to the change of
mortality, and longevity assessments. This is well-documented in the works of
Shimizu et al. [13] and Shirai and Shimizu [14].

In employing these models, a suitable choice of $U_{c},V_{c}$, and
$\Lambda_{c}$ is important. Occasionally, the predictive performance can be
influenced by the hyperparameters embedded within these functions. While a
well-constructed parametric model has demonstrated the potential to outperform
not only the classical Lee-Carter model [6] but also other established cohort
models such as the CBD model by Cairns et al. [1] and the RH model by Renshaw
and Haberman [8, 9] — as detailed in Shimizu et al. [11, 12] numerical
comparisons — the converse outcome can easily transpire if executed
improperly. The selection of such parametric models requires a universal
measure, like an information criterion. However, constructing such a criterion
poses challenges since $X^{c}$ represents a fictitious quantity that was
originally unobservable. To overcome this difficulty, embracing the
contemporary paradigm of Functional Data Analysis (FDA) holds great promise.

In this manuscript, we aim to estimate these functions and generate forecasts
for prospective cohorts through the application of the FDA. First, we employ
nonparametric estimation of the functions $U_{c},V_{c}$, and $\Lambda_{c}$,
using the least-squares method at specific time points $t=t_{1},\dots,t_{d}$.
Plotting these estimates with respect to $t$, achieved through appropriate
smoothing techniques, furnishes us with the functional data for $U_{c},V_{c}$
and $\Lambda_{c}$. This process is undertaken for each cohort, facilitating
our comprehension of how these functional data evolve with cohort progression.
Subsequently, by applying the FDA’s guidelines to these developments, we can
predict the functions $U_{c},V_{c}$, and $\Lambda_{c}$ for forthcoming
cohorts. These predictions are then complemented by our original refinements,
culminating in the computation of a predictive mortality function within each
SEM.

Our research shows that using SEM predictions with the FDA is better than
using parametric frameworks. Importantly, the FDA technique we use is a widely
accepted and easy-to-use method available in the R programming package.

In Section 2, we explain the specific steps for estimating and predicting
using the FDA for some critical functions, say ‘key’ functions, that represent
each SEM. It is important to note that the predictions from the ‘key’ function
determine the projected mortality function. We also provide an overview of
Functional Principal Component Analysis (FPCA) in Subsection 2.2, as it
constitutes a primary tool in our forecasting procedure. Section 3 is
dedicated to the analysis of real data, sourced from the Human Mortality
Database (HMD) [5], which is a prominent open-source repository. Our analysis
conclusively demonstrates the superior performance of our proposed methodology
in contrast to parametric alternatives. The paper is closed with concluding
remarks, presented in Section 4, in which we synthesize our findings and offer
insights into the broader implications of our research.

##  2 Prediction procedure with FDA

###  2.1 Functional datatization of ‘key’ functions in SEM

Our initial undertaking involves the process of ”datatization” of the ‘key’
functions, which constitute unobservable entities within the SEM. The
fundamental concept is to estimate these functions at discrete time points for
each cohort and subsequently apply a basis expansion to transform them into
(continuous) functional data.

Recall the explicit expressions for mortality functions given in (1.1) in ID-
and IG-SEM.

ID-SEM: The mortality function of (1.2) for cohort $c$ is given as follows
(Shimizu et al. [11]).

$\displaystyle
q_{c}^{ID}(t;M_{c})=1-\Phi\left(\frac{x_{c}+M_{c}(t)}{\sqrt{2S_{c}(t)}}\right)+e^{-\kappa_{c}x_{c}}\Phi\left(\frac{-x_{c}+M_{c}(t)}{\sqrt{2S_{c}(t)}}\right),$
(2.1)

where $\Phi$ is the cumulative distribution function of the standard normal
distribution, $M_{c}(t)=\int_{0}^{t}U_{c}(s)\,\mathrm{d}s$ and
$S_{c}(t)=\frac{1}{2}\int_{0}^{t}V_{c}^{2}(s)\,\mathrm{d}s$ under the
assumption that

$\displaystyle\frac{M_{c}(t)}{S_{c}(t)}\equiv\kappa_{c}<0\ (\mbox{constant}).$
(2.2)

Due to the last assumption, it suffices to know the negative constant
$\kappa_{c}$ and, e.g., the function $M_{c}$ (or $S_{c}$) to have
$q_{c}^{ID}$. Therefore, we consider $M_{c}$ a ‘key’ function in ID-SEM.

IG-SEM: The mortality function of (1.3) for cohort $c$ is given as follows
(Shimizu et al. [12]).

$\displaystyle
q_{c}^{IG}(t;\Lambda_{c})=\Phi\left(\sqrt{\frac{\sigma_{c}}{x_{c}}}(\Lambda_{c}(t)-x_{c})\right)-e^{2\sigma_{c}\Lambda_{c}(t)}\Phi\left(-\sqrt{\frac{\sigma_{c}}{x_{c}}}(\Lambda_{c}(t)+x_{c})\right).$
(2.3)

In this model, we need to know the constant $\sigma_{c}$ and the function
$\Lambda_{c}$, which is the ‘key’ function in IG-SEM.

######  Remark 2.1.

In this paper, we adopt a fixed approach with respect to the parameters
$\kappa_{c}$ in the ID-SEM and $\sigma_{c}$ in the IG-SEM, setting them to
specific values that are determined through empirical selection. These values
are referenced in (3.1) within the data analysis in Section 3. The empirical
selection of these values is informed by the outcomes of parametric inference,
with the primary objective being the attainment of stability in the estimation
of $M_{c}$ or $\Lambda_{c}$. While the determination of $\kappa_{c}$ and
$\sigma_{c}$ remains an important consideration, we do not discuss it further
in this paper, but it is worth noting that even if these values are chosen in
an ad hoc manner, the predictive efficacy of the mortality function remains
good, provided that the ‘key’ functions can be accurately estimated.

Our primary objective is to acquire the ”functional data” pertaining to the
‘key’ functions. We denote an empirical mortality function by
$q_{c}^{DATA}(t)$ $(t=S,S+1,\dots,w;\ c=c_{1},c_{2},\dots,c_{m})$. These
empirical functions are derived from data in the Human Mortality Database, as
detailed in Shimizu et al. [11]. Here, $w$ denotes the final age recorded in
the Human Mortality Database, and the value of $S$ is determined on a case-by-
case basis by the observer.

In the sequel, we will explain the procedure to obtain functional data for
$M_{c}$ in the case of ID-SEM. The one for $\Lambda_{c}$ in the case of IG-SEM
is the same. Moreover, in practice, we use the conditional version of the
mortality function following the manner of Shimizu et al. [11]: e.g., in the
case of ID-SEM,

$\displaystyle q^{ID}_{c}(t,{\color[rgb]{0,0,0}M_{c}}|S)$
$\displaystyle:=\frac{q_{c}^{ID}(t;{\color[rgb]{0,0,0}M_{c}})-q_{c}^{ID}(S;{\color[rgb]{0,0,0}M_{c}})}{1-q_{c}^{ID}(S;{\color[rgb]{0,0,0}M_{c}})}={\mathbb{P}}(\tau_{c}\leq
t|\tau_{c}>S),\quad t>S,$

and its empirical version

$\displaystyle q^{DATA}_{c}(t|S)$
$\displaystyle:=\frac{q_{c}^{DATA}(t)-q_{c}^{DATA}(S)}{1-q_{c}^{DATA}(S)},\quad
t>S.$

######  Remark 2.2.

We assume we have mortality data for individuals born in cohort $c_{m}$ until
they reach the age of $w$ years. In other words, we have access to mortality
data extending up to the calendar year $c_{m+w}$. Consequently, for cohort
$c_{m+L}$, we have mortality data available up to the age of $w-L$ years, and
our objective for this specific cohort is the prediction of the mortality
function $q_{c}(t|S)$ for $t>(w-L)\vee S$

Smoothing for discretely estimated ‘key’ functions:

  * •

Estimate $M_{c}(t)$ for each $t=S,S+1,\dots,w;\ c=c_{1},c_{2},\dots,c_{m}$ by
the least squares method:

$\widehat{M}_{c}(t)=\arg\min_{M\in\mathbb{R}}\left|q^{ID}_{c}(t,M|S)-q^{DATA}_{c}(t|S)\right|^{2}.$

Note that the error on the right-hand side can be zero. In this stage, our
focus is on determining the value of $M$ at the specific time point $t$ where
$q_{c}^{ID}(t;M|S)\approx q_{c}^{DATA}(t|S)$. In essence, we are seeking
discrete observations of the latent function $M$ at the moment when the
identified quantile function closely approximates the data quantile function.

  * •

For each cohort $c=c_{1},c_{2},\ldots,c_{m}$, we treat the values
$\\{M_{c}(t)\,|\,t=S,S+1,\dots,w\\}$ as discrete data sampled from the
trajectory of $(M_{c}(t))_{S\leq t\leq w}$, and proceed to interpolate these
values using continuous functions $h_{l}:\mathbb{R}+\to\mathbb{R}$, which form
an orthogonal basis in the space $L^{2}([S,w])$. In practical terms, we
truncate the infinite series, limiting it to a suitably “large” value of $L$:
for $c=c_{1},c_{2},\dots,c_{m}$ and
${\color[rgb]{0,0,0}\alpha_{c}=(\alpha_{c,1},\dots,\alpha_{c,L})}$,

$\displaystyle
M_{c}(t;{\color[rgb]{0,0,0}\alpha_{c}}):=\sum_{l=1}^{L}{\color[rgb]{0,0,0}\alpha_{c,l}}h_{l}(t),\quad
S\leq t\leq w.$

######  Remark 2.3.

Regarding the orthonormal basis $\\{h_{l}\\}_{l=1,2,\dots}$, the choice of the
Fourier series expansion is typically suitable for handling periodic data,
whereas $B$-splines are preferred for non-periodic data. In our research, we
use the $B$-spline expansion method.

The coefficient ${\color[rgb]{0,0,0}\alpha_{c,l}}$ in each summand is
estimated by the least squares method:

${\color[rgb]{0,0,0}\widehat{\alpha}_{c}}:=({\color[rgb]{0,0,0}\widehat{\alpha}_{c,1},\dots,\widehat{\alpha}_{c,L}})=\arg\min_{\alpha\in\mathbb{R}^{L}}\sum_{t=S}^{w}\left|M_{c}(t;\alpha)-\widehat{M}_{c}(t)\right|^{2}.$

In this way, we have $m$-discrete obervations as functional data for the ‘key’
function $M_{c}$ in ID-SEM, and we obtain a basis expnsion of centerd data
$\widehat{M}^{*}_{c}=(\widehat{M}^{*}_{c}(t))_{S\leq t\leq w}$ such that

$\displaystyle\widehat{M}^{*}_{c}(t)=M_{c}(t;{\color[rgb]{0,0,0}\widehat{\alpha}_{c}})-\overline{M}(t)=:\sum_{l=1}^{L}{\color[rgb]{0,0,0}\widehat{\alpha}^{*}_{c,l}}h_{l}(t),$
(2.4)

where
$\overline{M}(t)=\frac{1}{m}\sum_{c=c_{1}}^{c_{m}}M_{c}(t;{\color[rgb]{0,0,0}\widehat{\alpha}_{c}})$,
and
${\color[rgb]{0,0,0}\widehat{\alpha}^{*}_{c,l}=\widehat{\alpha}_{c,l}-\frac{1}{m}\sum_{c=c_{1}}^{c_{m}}\widehat{\alpha}_{c,l}}$.

Similarly, we also have centred functional data of
$\widehat{\Lambda}^{*}_{c}=(\widehat{\Lambda}^{*}_{c}(t))_{S\leq t\leq w}$ in
IG-SEM, where

$\displaystyle\widehat{\Lambda}^{*}_{c}(t)=\Lambda_{c}(t;{\color[rgb]{0,0,0}\widehat{\beta}_{c}})-\overline{\Lambda}(t):=\sum_{l=1}^{L}{\color[rgb]{0,0,0}\widehat{\beta}^{*}_{c,l}}h_{l}(t),$
(2.5)

with
${\color[rgb]{0,0,0}\widehat{\beta}_{c}=(\widehat{\beta}_{c,1},\dots,\widehat{\beta}_{c,L})}$
and
${\color[rgb]{0,0,0}\widehat{\beta}^{*}_{c,l}=\widehat{\beta}_{c,l}-\frac{1}{m}\sum_{c=c_{1}}^{c_{m}}\widehat{\beta}_{c,l}}$.

###  2.2 Functional principle component analysis

Functional Principal Component Analysis (FPCA) is an extension of classical
Principal Component Analysis (PCA) to handle random elements that take on
functional values.

Consider a stochastic process $Y_{c}=(Y_{c}(t))_{t\in[0,T]}$ with continuous
paths, where $T>0$. These processes are indexed by $c=1,2,\dots,m$, and their
values reside in a Hilbert space $L^{2}([0,T])$. We assume that the random
sequence $Y_{1},Y_{2},\dots,Y_{m}$ comes from weakly stationary function-
valued time series, with a mean function $\mu(t)=\mathbb{E}[Y_{1}(t)]$. In
this context, the covariance kernel, denoted as
$K(t,s):=\mathrm{Cov}(Y^{c}(t),Y^{c}(s))$, is invariant with respect to $c$.
It has a spectral-type decomposition, which can be characterized using
Mercer’s Theorem (see, e.g., Riesz and Sz.-Nagy [10]):

$K(t,s)=\sum_{j=1}^{\infty}\lambda_{j}e_{j}(t)e_{j}(s),\quad\lambda_{1}\geq\lambda_{2}\geq\dots\geq\lambda_{n}\geq\dots,$

where $\lambda_{j}$ is the $j$th eigenvalue and $e_{j}$ is the corresponding
eigenfunction in $L^{2}([0,T])$ space. Then the Karhunen-Lòeve expansion
allows us to decompose the random function $Y_{c}(\cdot)$ as

$\displaystyle
Y^{*}_{c}(t)=\sum_{j=1}^{\infty}{\color[rgb]{0,0,0}Z_{cj}}e_{j}(t),\quad{\color[rgb]{0,0,0}Z_{cj}}=\langle
Y_{c}-\mu,e_{j}\rangle_{L^{2}([0,T])},$

where $Y^{*}_{c}(t):=Y_{c}(t)-\mu(t)$ and ${\color[rgb]{0,0,0}Z_{cj}}$ is the
inner product in $L^{2}([0,T])$ of $Y^{*}_{c}$ and $e_{j}$:

${\color[rgb]{0,0,0}Z_{cj}}:=\int_{0}^{T}Y^{*}_{c}(t)e_{j}(t)\,\mathrm{d}t,$

which corresponds to the principal component score in the classical (finite-
dimensional) PCA. According to this manner, we can assume that

$\displaystyle
Y_{c}(t)=\mu(t)+\sum_{j=1}^{L}{\color[rgb]{0,0,0}Z_{cj}}e_{j}(t)+{\epsilon}_{c}(t),$
(2.6)

where $L$ is an integer suitably chosen and ${\epsilon}_{c}$ is a noise
process depending on the index $c$.

Our aim is to find an approximation of the expansion (2.6). First, we shall
estimate the covariance kernel $K$ by a sample version:

$\widehat{K}(s,t)=\frac{1}{m-1}\sum_{c=1}^{m}Y^{*}_{c}(s)Y^{*}_{c}(t),$

and we would like to find the eigenvalues $\lambda_{j}$’s and eigenfunctions
$e_{j}$’s of $K$:

$\int_{0}^{T}K(s,t)e_{j}(s)\,\mathrm{d}s=\lambda_{j}e_{j}(t),\quad t\in[0,T].$

This practical procedure is seen in Ramsay and Silverman [7], Section 8.4.2,
or a more general argument is found in Hyndman and Ullah [4]. With the
expansion of our data using an orthonormal system, as detailed in (2.4) and
(2.5), we can now consider our data in a centered form, expressed as follows:

$\widehat{Y}^{*}_{c}(\cdot):=Y_{c}(\cdot)-\overline{Y}(\cdot)=\sum_{l=1}^{L}{\color[rgb]{0,0,0}\widehat{\alpha}^{*}_{c,l}}\,\xi_{j}(\cdot),\quad
c=1,2,\dots,m,$

where $\overline{Y}(t)=\frac{1}{m}\sum_{c=1}^{m}Y_{c}(t)$, and
$\\{\xi_{j}\\}_{j=1\dots}$ is an orthonormal basis in $L^{2}([0,T])$.

Let
$\mbox{\boldmath$A$}=({\color[rgb]{0,0,0}\widehat{\alpha}^{*}_{c,l}})_{1\leq
c\leq m;1\leq l\leq L}$ is the $m\times L$ matrix of coefficients. Following
the manner of Ramsay and Silverman [7], the eigenfunctions $\widehat{e}_{j}$’s
of the kernel function $\widehat{K}$ has the following basis expansion:

$\widehat{e}_{j}(t)=\sum_{k=1}^{L}b^{j}_{k}\xi_{k}(t),\quad j=1,2,\dots,L,$

where $\mbox{\boldmath$b$}^{(j)}=(b^{j}_{1},\dots,b^{j}_{L})$ is the $j$th
eigenvector of $(m-1)^{-1}\mbox{\boldmath$A$}^{\top}\mbox{\boldmath$A$}$. The
corresponding eigenvalue becomes $\widehat{\lambda}_{j}$:

$\frac{1}{m-1}\mbox{\boldmath$A$}^{\top}\mbox{\boldmath$A$}\mbox{\boldmath$b$}^{(j)}=\widehat{\lambda}_{j}\mbox{\boldmath$b$}^{(j)}.$

Then, the estimated principal component score is given by

$\displaystyle{\color[rgb]{0,0,0}\widehat{Z}_{cj}}=\int_{0}^{T}\widehat{Y}^{*}_{c}(t)\widehat{e}_{j}(t)\,\mathrm{d}t,\quad
c=1,2,\dots,m.$

As a consequence, we have an approximation of (2.6) as follows: for
$L\in\mathbb{N}$ large enough,

$\displaystyle\widehat{Y}^{*}_{c}(t)=\sum_{j=1}^{L}{\color[rgb]{0,0,0}\widehat{Z}_{cj}}\,\widehat{e}_{j}(t).$

For the theoretical details for FPCA and the inference, see also Horváth and
Kokoszka [2], or Hsing and Eubank [3].

###  2.3 Prediction of the key function for future cohorts

We implement the procedure outlined in the preceding section, applying it to
the functional data $Y_{c}=M_{c}$ and $Y_{c}=\Lambda_{c}$. These datasets take
the form:

$\\{\widehat{M}^{*}_{c_{1}},\dots,\widehat{M}^{*}_{c_{m}}\\}\quad\mbox{and}\quad\\{\widehat{\Lambda}^{*}_{c_{1}},\dots,\widehat{\Lambda}^{*}_{c_{m}}\\}$

as in (2.4) and (2.5), respectively. We explain mainly the procedures of
estimation and prediction for $M_{c}$. It is important to note that the
procedures for $\Lambda_{c}$ follows precisely the same way.

Let us start by recalling that we have the smoothly estimated functional data:

$\\{\widehat{M}^{*}_{c}(t)\,|\,t=S,S+1,\dots w;\
c=c_{1},c_{2},\dots,c_{m}\\}.$

Step 1: Approximating the Karhunen-Lòeve expansion.

Following the steps in the previous subsection, we have the eigenfunctions
$\widehat{e}_{1j}$’s and the corresponding eigenvalues $\lambda_{1j}$’s of

$\widehat{K}_{1}(s,t)=\frac{1}{m-1}\sum_{c=c_{1}}^{c_{m}}\widehat{M}^{*}_{c}(s)\widehat{M}^{*}_{c}(t).$

Without loss of generality, we may assume that
$\widehat{\lambda}_{1}\geq\widehat{\lambda}_{2}\geq\dots\geq\widehat{\lambda}_{L}$.
Then, we have

$\widehat{M}^{*}_{c}(t)=\sum_{j=1}^{L}{\color[rgb]{0,0,0}\widehat{Z}^{(1)}_{cj}}\widehat{e}_{1j}(t),$

where
${\color[rgb]{0,0,0}\widehat{Z}^{(1)}_{cj}}=\int_{0}^{T}\widehat{M}^{*}_{c}(t)\widehat{e}_{1j}(t)\,\mathrm{d}t$
is the $j$th FPC-score. The number $L$ is numerically chosen as
$L=\mathrm{rank}\widehat{K}_{1}$, where
$\widehat{K}_{1}=(\widehat{K}_{1}(t,s))_{t,s=S,S+1,\dots,w}$, and we can
reduce it so that the contribution rate of the score exceeds a certain
threshold $\vartheta\in(0,1)$ (usually close to 1):

$\displaystyle
K_{1\vartheta}=\inf\left\\{k\in\mathbb{N}\,\Big{|}\,\frac{\sum_{j=1}^{k}\widehat{\lambda}_{1j}}{\sum_{j=1}^{L}\widehat{\lambda}_{1j}}\geq\vartheta\right\\}.$
(2.7)

Then the estimated (approximated) $M_{c}$ is given by

$\displaystyle\widehat{M}_{c}(t;{\color[rgb]{0,0,0}\widehat{\alpha}_{c}})\approx\overline{M}(t)+\sum_{j=1}^{K_{1\vartheta}}{\color[rgb]{0,0,0}\widehat{Z}^{(1)}_{cj}}\widehat{e}_{1j},\quad
c=c_{1},\dots,c_{m}.$

Following the identical procedure used for
$\widehat{M}_{c}(t;{\color[rgb]{0,0,0}\widehat{\alpha}_{c}})$, we arrive at an
approximation for $\Lambda_{c}$ as follows:

$\displaystyle\widehat{\Lambda}_{c}(t;{\color[rgb]{0,0,0}\widehat{\beta}_{c}})\approx\overline{\Lambda}(t)+\sum_{j=1}^{K_{2\vartheta}}{\color[rgb]{0,0,0}\widehat{Z}^{(2)}_{cj}}\widehat{e}_{2j}(t),\quad
c=c_{1},\dots,c_{m}.$

where
${\color[rgb]{0,0,0}\widehat{Z}^{(2)}_{cj}}=\int_{0}^{T}\widehat{\Lambda}^{*}_{c}(t)\widehat{e}_{2j}(t)\,\mathrm{d}t$,
$\widehat{e}_{2j}$’s are the eigenfunctions of
$\widehat{K}_{2}(s,t)=\frac{1}{m-1}\sum_{c=c_{1}}^{c_{m}}\widehat{\Lambda}^{*}_{c}(s)\widehat{\Lambda}^{*}_{c}(t).$
The selection of $K_{2\vartheta}$ is determined in a manner similar to that
outlined in (2.7).

Step 2: Prediction of mortality function.

We now possess a set of time series data for the principal component scores,
represented as
$\\{{\color[rgb]{0,0,0}\widehat{Z}^{(1)}_{c_{1}j},\widehat{Z}^{(1)}_{c_{2}j},\dots,\widehat{Z}^{(1)}_{c_{m}j}}\\}$.
Considering this dataset as historical data, our objective is to make
predictions for $\\{{\color[rgb]{0,0,0}\widehat{Z}^{(1)}_{cj}}\,|\,c>c_{m}\\}$
for each $j=1,\dots,K_{1\vartheta}$. To accomplish this, we will employ time
series models and forecast future values based on the distribution of these
scores. This approach aligns with the methodology described in Hyndman and
Ullah [4].

######  Remark 2.4.

Subsequently, we will employ the ARIMA model to fit the data and utilize the
95% prediction intervals to forecast the SEM mortality function for future
cohorts. The process of implementing ARIMA predictions is straightforward,
particularly when using the R package forecast.

Once we obtain the predictor
$\\{{\color[rgb]{0,0,0}\widehat{Z}^{(1)}_{cj}}\,|\,c>c_{m}\\}$ for each
$j=1,\dots,K_{1\vartheta}$, we have the (conditional) mortality functions of
future cohorts $c>c_{m}$:

$\displaystyle
q_{c}^{ID}(t,\widehat{M}_{c}(\cdot;{\color[rgb]{0,0,0}\widehat{\alpha}_{c}})|S),\quad
c=c_{m+1},c_{m+2},\dots.$

Regarding IG-SEM, we take the same procedures as above. We employ ARIMA models
to forecast the future values of
$\\{{\color[rgb]{0,0,0}\widehat{Z}^{(2)}_{cj}}\,|\,c>c_{m}\\}$ for each
$j=1,\dots,K_{2\vartheta}$ and obtain

$\displaystyle
q_{c}^{IG}(t,\widehat{\Lambda}_{c}(\cdot;{\color[rgb]{0,0,0}\widehat{\beta}_{c}})|S),\quad
c=c_{m+1},c_{m+2},\dots.$

Step 3: Modification of predicted mortality function.

This modification step is a specialized technique unique to the SEM approach
and deviates from conventional forecasting procedures. In this regard, it
represents a novel contribution to the field.

When performing our estimations, we utilize (empirical) mortality data up to
the age of $w$ for the $c_{m}$-cohort. Therefore, when we seek to estimate
mortality for a cohort $c\,(>c_{m})$, individuals born in that cohort have
already reached the age of $c_{m}+w-c$, and we have access to mortality data
up to this age. As a result, we can draw upon information from actual
mortality data up to the age of $c_{m}+w-c$ for the cohort $c$ under
consideration to enhance our predictions.

In the previous step, where we predicted
$\\{{\color[rgb]{0,0,0}\widehat{Z}^{(k)}_{cj}}\,|\,c>c_{m}\\}$ using the ARIMA
model, we also obtained the corresponding prediction intervals. These
intervals, denoted as $I_{\delta}^{c,k,j}$, represent the $\delta$-prediction
intervals:

${\mathbb{P}}\left({\color[rgb]{0,0,0}\widehat{Z}^{(k)}_{cj}}\in
I_{\delta}^{c,k,j}\right)=\delta,\quad(c>c_{m};\ k=1,2;\
j=1,\dots,K_{k\vartheta})$

For instance, with $\delta=0.95$, the true score $Z^{c}_{kj}$ falls into the
$I_{\delta}^{c,k,j}$ with a probability of $\delta$. This assumption allows us
to fine-tune the fitting of the mortality function, as exemplified below for
the case where $k=1$:

$({\color[rgb]{0,0,0}\widetilde{Z}^{(1)}_{c1},\dots,\widetilde{Z}^{(1)}_{cK_{1\vartheta}}})=\arg\min_{{\color[rgb]{0,0,0}\widehat{Z}_{c}}}\sum_{t=S}^{c_{m+w-c}}\left|q_{c}^{ID}(t,\widehat{M}_{c}(\cdot;{\color[rgb]{0,0,0}\widehat{\alpha}_{c}})|S)-q_{c}^{DATA}(t|S)\right|^{2},$

where $\arg\min_{{\color[rgb]{0,0,0}\widehat{Z}_{c}}}$ is taken over all
${\color[rgb]{0,0,0}\widehat{Z}^{(1)}_{c}}$ satisfying that

${\color[rgb]{0,0,0}\widehat{Z}^{(1)}_{c}:=(\widehat{Z}^{(1)}_{c1},\dots,\widehat{Z}^{(1)}_{cK_{1\vartheta}})\in
I_{\delta}^{c,1,1}\times\dots\times I_{\delta}^{c,1,K_{1\vartheta}}}.$

Thus, we have the modified version of the conditional mortality function:

$\displaystyle q_{c}^{ID}(t,\widetilde{M}_{c}^{\delta}|S),\quad
c=c_{m+1},c_{m+2},\dots,c_{m+w-S}.$ (2.8)

with the modified ‘key’ function

$\widetilde{M}_{c}^{\delta}(t)=\overline{M}(t)+\sum_{j=1}^{K_{1\vartheta}}{\color[rgb]{0,0,0}\widetilde{Z}^{(1)}_{cj}}\widehat{e}_{1j}(t).$

In the same way, we also have

$\displaystyle q_{c}^{IG}(t,\widetilde{\Lambda}_{c}^{\delta}|S),\quad
c=c_{m+1},c_{m+2},\dots,c_{m+w-S}.$ (2.9)

with the modified ‘key’ function

$\widetilde{\Lambda}_{c}^{\delta}(t)=\overline{\Lambda}(t)+\sum_{j=1}^{K_{2\vartheta}}{\color[rgb]{0,0,0}\widetilde{Z}^{(2)}_{cj}}\widehat{e}_{2j}(t).$

######  Remark 2.5.

The above modification step is applicable to cohorts
$c=c_{m+1},c_{m+2},\ldots,c_{m+w-S}$ because no mortality data are available
beyond the age of $S$ years for cohorts beyond $c_{m+w-S}$.

##  3 Data analysis

In this section, we leverage the modified mortality functions described in
(2.8) and (2.9) to make predictions for selected future cohorts. In our
comparative analysis, we evaluate not only these two modified approaches but
also consider parametric models proposed in previous works by Shimizu et al.
[11, 12].

###  3.1 Swedish mortality

The Swedish mortality data available in the Human Mortality Database are among
the oldest records, rendering them well-suited for long-term projections. We
make use of empirical mortality functions $q_{c}^{DATA}(t)$ for cohorts
$c_{1}=1781,c_{2}=1782,\ldots,c_{m}=1830$ (with $m=50$) at ages
$t=20,21,\ldots,110$ (with $S=20$ and $w=110$). Our objective is to forecast
future mortality functions $q_{c}^{ID}(t|S)$ and $q_{c}^{IG}(t|S)$ for cohorts
$c=1870,1890$, and $1910$ (which correspond to 40, 60, and 80 years following
the last cohort, $c_{m}=1830$).

It is important to note that, as discussed in Remark 2.2, we already have
mortality data available up to the year $c_{m}+w=1940$ (for the purpose of
this analysis, we assume the current year to be 1940). Consequently,
individuals born in cohort $c$ are already $(1940-c)$ years old, and our
forecasting takes into account the age $(1940-c)+1$, assuming the individual
has lived up to age $S$. For example, when predicting the cohort $c=1910$, our
analysis pertains to individuals over the age of 31, and we compare the curve
$q_{c}(t|30)$ (for $31\leq t\leq 110)$.

The parametric models for key functions used for comparison below are the
following setup, similar to Shimizu et al. [12]:

  * •

ID-SEM (1.2):
$U_{c}(t,\vartheta_{c})=\alpha_{c}+\beta_{c}\exp\left(\gamma_{c}(t-T_{c})\right)\mathbf{1}_{\\{t\geq
T_{c}\\}}$, where $\vartheta_{c}=(\alpha_{c},\beta_{c},\gamma_{c})$, and
$T_{c}$ is a change point where the mortality changes drastically, fixed at
$T_{c}={\color[rgb]{0,0,0}50}$. Note that it is too hard to estimate and
predict $T_{c}$ in practice; see Shimizu et al. [11].

  * •

IG-SEM (1.3): $\Lambda_{\vartheta_{c}}(t)=e^{a_{c}t}+b_{c}t-1$, where
$\vartheta_{c}=(a_{c},b_{c},\sigma_{c})$.

All parameter estimation and forecasting methods were performed according to
the methodologies by Shimizu et al. [12].

Subsequently, we maintain the values of $x_{c}$ representing the initial
energy, $\kappa_{c}$ in (2.2), and $\sigma_{c}$ in (2.3) at a fixed set of
values as follows:

$\displaystyle x_{c}\equiv
1000;\quad\kappa_{c}\equiv-0.25;\quad\sigma_{c}\equiv 0.001,$ (3.1)

These specific values have been empirically selected based on insights
obtained from data analysis, as detailed in Shimizu et al. [12, 13]. They are
chosen to ensure that the estimation and optimization processes remain
numerically stable across all cohorts.

###  3.2 Forecasting of $M_{c}$ and $\Lambda_{c}$

We employed the R package fda to estimate the eigenfunctions
$\widehat{e}_{kj}$ and eigenvalues $\widehat{\lambda}_{kj}$ $(k=1,2)$. In
determining the number of terms, $K_{1\vartheta}$ in equation (2.7), we set
$\vartheta=0.995$, which yielded $K_{1\vartheta}=3$. The output includes the
mean function $\mu(t)=\overline{M}(t)$ and the eigenfunctions
$\widehat{e}_{1j}$ for $j=1,2,3$, as displayed in Figure 1 (a). The same
procedure for $\Lambda_{c}$ obtains Figure 1 (b).

The predictions for the PCA scores
${\color[rgb]{0,0,0}\widehat{Z}^{(1)}_{cj}}$ for cohorts $c>c_{m}$ are
conducted using forecast with ARIMA models. The results of these predictions
are given in Figure 2, (a)-(c), and the forecasted trajectories of
$\widetilde{M}^{\delta}_{c}(t)$ modified based on the prediction intervals are
illustrated in (d).

To obtain a modified key function, $\widetilde{M}^{\delta}_{c}(t)$, we employ
the $\delta=95$%-prediction interval (represented by the grey range in the
Figure 2 (a)-(c)). The graphs reveal that the function for $M_{c}$ exhibits a
tendency to increase year by year. Given that $M_{c}$ originally represents a
drift of survival energy, this suggests that the trajectory of survival energy
is on an upward trajectory year by year. Consequently, future cohorts are less
likely to experience mortality.

Similarly, we employ analogous procedures to forecast $\Lambda_{c}$, and the
corresponding results are presented in Figure 4.

###  3.3 Parametric vs. Nonparametric

Using the key function $M_{c}$ derived from the preceding step, we compute a
predicted mortality function and compare it with the results of the parametric
model proposed by Shimizu et al. [12] (using the same model setup as in their
paper). Our predictions are cohorts born in $c=1870,1890$, and $1910$. In our
data configuration, we used data up to the cohort born in 1830 (covering ages
20-110), placing us in the year 1940. Consequently, our predictions target
individuals from the $c=1870,1890$, and $1910$ cohorts at ages 70+, 50+, and
30+, respectively. Figures 3 and 5 display $q^{ID}_{c}(t|30)$ and
$q^{IG}_{c}(t|30)$ for the male cohort born in $c=1910$, which represents the
longest projection. Additionally, Tables 1 and 2 provide a comparison of the
Mean Squared Error (MSE) of males and females predicted mortality functions
for these cohorts: for $c=1870,1890,1910$ with $\widetilde{c}=1940-c$, and for
$w=110$,

$MSE=\frac{1}{w-\widetilde{c}+1}\sum_{t=\widetilde{c}}^{w}\left|q_{c}^{**}(t|30)-q_{c}^{DATA}(t|30)\right|^{2},\quad**=ID\
\ \mbox{or}\ \ IG$

The outcomes clearly indicate that our FDA approach outperforms the parametric
models (except for a few rare cases) and offers superior accuracy in mortality
prediction. As illustrated in Figures 3 and 5, our methodology demonstrates a
propensity for heightened accuracy in mortality prediction, particularly when
applied to the elderly demographic. Furthermore, this nonparametric approach
eliminates the need to select change points $T_{c}$ required in ID-SEM. This
constitutes a significant practical contribution, as accurate forecasting for
the elderly demographic carries substantial overarching significance.

Occasionally, ID-SEM outperforms the FDA approach in terms of predictions
(e.g., 1890-cohort for females). Nevertheless, the selection of the
hyperparameter $T_{c}$ in ID-SEM bears significant importance. An erroneous
choice in this regard can lead to a sudden deterioration in predictions,
rendering the parametric approach notably unstable.

######  Remark 3.1.

As discussed by Shimizu et al. [12], the introduced ”Modification” in Section
2.3 does not invariably lead to improved predictions. This is attributed to
the risk of entering an ”over-learning” state, induced by excessive LSE-
fitting to existing data, which subsequently diminishes the model’s predictive
capability for future data. One potential remedy involves the reduction of the
prediction threshold, e.g., $\delta=0.8$, etc. Nevertheless, based on our
extensive experience in numerous data analyses, it is evident that, in most
cases, this ”Modification” indeed results in a reduction of MSE in the
predictions.

Males | ID-SEM |  | IG-SEM |
---|---|---|---|---
Cohort | Parametric | FDA | Parametric | FDA
1870 | 6.6747E-4 | 9.6269E-5 | 5.1827E-4 | 7.3466E-6
1890 | 3.4059E-3 | 3.2618E-3 | 2.6057E-2 | 4.7426E-4
1910 | 1.8461E-3 | 7.6211E-4 | 4.8067E-2 | 1.2011E-3
Table 1: Parametric vs. FDA (Males): MSEs for predictions of Swedish mortality functions of cohort $c=1870,1890$ and $1910$, in ID-SEM and IG-SEM.  Females | ID-SEM |  | IG-SEM |
---|---|---|---|---
Cohort | Parametric | FDA | Parametric | FDA
1870 | 9.7386E-4 | 2.8153E-4 | 1.1994E-3 | 2.1915E-4
1890 | 1.0991E-4 | 1.2383E-4 | 1.8389E-2 | 1.1757E-3
1910 | 5.8036E-3 | 1.0184E-3 | 1.8928E-2 | 1.1845E-2
Table 2: Parametric vs. FDA (Females): MSEs for predictions of Swedish
mortality functions of cohort $c=1870,1890$ and $1910$, in ID-SEM and IG-SEM.

##  4 Concluding remarks

We used a nonparametric approach with functional data analysis for cohort-
specific mortality rate prediction based on Survival Energy Modeling (SEM). In
contrast to existing research with parametric models, which lack objective
model selection criteria and make it difficult to identify data-driven models,
our FDA approach allows for easy estimation and prediction of key functions.
In fact, this approach was able to significantly outperform predictions made
by existing parametric models in both the ID-SEM and IG-SEM. Additionally, the
convenient environment of R package for the FDA makes these approaches readily
applicable in practice without the need for extensive theoretical expertise,
making it a valuable tool. While finding more powerful candidates than ID-SEM
or IG-SEM for future predictions is a potential challenge, our data analysis
suggests that ID-SEM already possesses sufficient predictive power at this
stage.

However, amending forecasts through the utilization of the prediction interval
associated with the PCA score can occasionally yield counterproductive
outcomes. This issue bears a resemblance to the concept of ”overlearning”,
where excessive adaptation to historical data results in inadequate
performance when applied to future data. One potential remedy involves the
reduction of the prediction threshold denoted as $\delta$. Nonetheless, based
on our empirical findings, even when errors worsen, they tend to be of minor
magnitude, and, in the majority of instances, such adjustments lead to an
enhancement in forecasting accuracy.

As discussed by Shimizu et al. [13], cohort-specific mortality rate prediction
with SEM can be a powerful tool for calculating insurance premiums for each
cohort, predicting changes in pension amounts, and financial forecasting.
Particularly in countries with increasing life expectancy, such as Japan, the
analysis of longevity risk is a crucial task, and SEM is likely to become a
valuable tool for addressing this challenge.

---
(a)
(b)
Figure 1: Figure (a) is the estimated mean functions $\overline{M}(t)$ and
eigenfunctions for $M^{*}_{c}(t;{\color[rgb]{0,0,0}\widehat{\alpha}_{c}})$.
Figure (b) is the one for $\overline{\Lambda}_{c}(t)$ and
$\Lambda^{*}_{c}(t;{\color[rgb]{0,0,0}\widehat{\beta}_{c}})$.

(a) (b)
---

(c) (d)
---

Figure 2: Figures (a)-(c) illustrate the predictions for the 1st to 3rd PCA
scores. In each figure, the blue line represents the mean, while the purple
and grey ranges denote the 80% and 95% prediction intervals, respectively.
Figure (d) provides the forecast for $\widetilde{M}^{\delta}_{c}$ for cohorts
$1840\leq c\leq 1910$ with $\delta=0.95$. The grey curves correspond to the
estimated past curves $\widehat{M}_{c}$ for cohorts $1781\leq c\leq 1830$.
Figure 3: Comparison between the Swedish Male ID-SEM $q_{1910}^{ID}(t|30)$
predictions using parametric (Shimizu et al. [12]) and nonparametric (FDA)
methods. It’s evident that the FDA approach outperforms the parametric one,
particularly in predicting the elderly population.

(a) (b)
---

(c) (d)
---

Figure 4: Figures (a)-(c) illustrate the predictions for the 1st to 3rd PCA
scores. In each figure, the blue line represents the mean, while the purple
and grey ranges denote the 80% and 95% prediction intervals, respectively.
Figure (d) provides the forecast for $\widetilde{\Lambda}^{\delta}_{c}$ for
cohorts $1840\leq c\leq 1910$ with $\delta=0.95$. The grey curves correspond
to the estimated past curves $\widehat{\Lambda}_{c}$ for cohorts $1781\leq
c\leq 1830$.  Figure 5: Comparison between the Swedish Male IG-SEM
$q_{1910}^{IG}(t|30)$ predictions using parametric (Shimizu et al. [12]) and
nonparametric (FDA) methods. In this case, the FDA method outperforms the
parametric one overall.

Acknowledgement. This work is partially supported by JSPS KAKENHI Grant-in-Aid
for Scientific Research (C) #21K03358.

Conflict of interest statement. On behalf of all authors, the corresponding
author states there is no conflict of interest.

## References

  * [1] Cairns, A. J. G.; Blake, D. and Dowd, K. (2006a). A two-factor model for stochastic mortality with parameter uncertainty: theory and calibration. Journal of Risk and Insurance, 73, 687–718.
  * [2] Horváth, L. and Kokoszka, P. (2012). Inference for functional data with applications, Springer New York, NY.
  * [3] Hsing, T. and Eubank, R. (2015). Theoretical Foundations of Functional Data Analysis, with an Introduction to Linear Operators, John Wiley & Sons, UK.
  * [4] Hyndman, R. J. and Ullah, M S. (2007). Robust forecasting of mortality and fertility rates: A functional data approach. Computational Statistics and Data Analysis, 51, 4942–4956.
  * [5] Human Mortality Database: https://www.mortality.org/
  * [6] Lee, R.D. and Carter, L. (1992). Modeling and forecasting the time series of US mortality, Journal of the American Statistical Association, 87, pp.659–671.
  * [7] Ramsay, J. O. and Silverman, B. W. (2005). Functional data analysis, Springer New York, NY.
  * [8] Renshaw, A.E. and Haberman, S. (2003). Lee-Carter mortality forecasting with age-specific enhancement. Insurance: Mathematics and Economics, 33, 255–272.
  * [9] Renshaw, A. and Haberman, S. (2006) A cohort-based extension to the Lee-Carter model for mortality reduction factors. Insurance: Mathematics and Economics, 38, 556–570.
  * [10] Riesz, F. and Sz.-Nagy, B. (1990). Functional analysis. Dover.
  * [11] Shimizu, Y.; Minami, Y. and Ito, R. (2020). Why does a human die? A structural approach to cohort-wise mortality prediction under the survival energy hypothesis. ASTIN bulletin, 51, (1) 191–219.
  * [12] Shimizu, Y.; Shirai, K.; Kojima, Y.; Mitsuda, D. and Inoue, M. (2023). Survival energy models for mortality prediction and the future prospects, ASTIN Bulletin, 53, (2), 377–391.
  * [13] Shimizu, Y; Shirai, K.; Kojima, Y.; Mitsuda, D. and Inoue, M. (2023). A new survival energy model and SEM project. SSRN:4127900: http://dx.doi.org/10.2139/ssrn.4127900.
  * [14] Shirai, K. and Shimizu, Y. (2022). On the prediction of full life expectancy via Survival Energy Models, Master thesis, Graduate School of Fundamental Science and Engineering, Waseda University (in Japanese); to appear in JARIP Journal.


ifaamas [AAMAS ’24]Proc. of the 23rd International Conference on Autonomous
Agents and Multiagent Systems (AAMAS 2024)May 6 – 10, 2024 Auckland, New
ZealandN. Alechina, V. Dignum, M. Dastani, J.S. Sichman (eds.) 2024 2024 933
Purdue University West Lafayette United States Purdue University West
Lafayette United States Purdue University West Lafayette United States

# Holonic Learning: A Flexible Agent-based Distributed Machine Learning
Framework

Ahmad Esmaeili<EMAIL_ADDRESS>,  Zahra Ghorrati<EMAIL_ADDRESS>and
Eric T. Matson<EMAIL_ADDRESS>
###### Abstract.

Ever-increasing ubiquity of data and computational resources in the last
decade have propelled a notable transition in the machine learning paradigm
towards more distributed approaches. Such a transition seeks to not only
tackle the scalability and resource distribution challenges but also to
address pressing privacy and security concerns. To contribute to the ongoing
discourse, this paper introduces Holonic Learning (HoL), a collaborative and
privacy-focused learning framework designed for training deep learning models.
By leveraging holonic concepts, the HoL framework establishes a structured
self-similar hierarchy in the learning process, enabling more nuanced control
over collaborations through the individual model aggregation approach of each
holon, along with their intra-holon commitment and communication patterns.
HoL, in its general form, provides extensive design and flexibility
potentials. For empirical analysis and to demonstrate its effectiveness, this
paper implements HoloAvg, a special variant of HoL that employs weighted
averaging for model aggregation across all holons. The convergence of the
proposed method is validated through experiments on both IID and Non-IID
settings of the standard MNISt dataset. Furthermore, the performance behaviors
of HoL are investigated under various holarchical designs and data
distribution scenarios. The presented results affirm HoL’s prowess in
delivering competitive performance particularly, in the context of the Non-IID
data distribution.

###### Key words and phrases:

Distributed learning, Holonic Learning, Collaborative Learning, Edge Computing

##  1\. Introduction

Today’s interconnected world, marked by the proliferation of information and
computational capabilities across vast digital landscapes, has led to critical
challenges concerning the scalability, adaptability, speed, and support for
diversity in the traditional centralized learning models. As a solution to
such limitations and harnessing the collective power of the networked
devices—whether servers or edge devices—recent decade has witnessed a
substantial growth in the rise of distributed learning methodologies that not
only expedite training process but also enhances robustness and privacy.

Federated Learning (FL) McMahan et al. (2017) stands out as one of the
prominent distributed learning approaches attracting attentions from both
academia and industry Yang et al. (2019). By empowering end users to
participate in the collaborative learning process while retaining ownership of
their data, FL ensures both data security and privacy on one hand, and
leveraging the potential of computationally constrained devices on the other.
Despite its apparent simplicity in its original format, FL introduces a series
of intricate open problems and challenges tackled by a plethora of research
work. Communication efficiency, convergence rates, privacy preservation,
addressing potential adversarial attacks, robustness across heterogeneous
devices and data distributions, non-identically and independently distributed
(non-IID) data, mitigating potential aggregation bias, and ensuring fairness
in model updates are among the nuanced and critical open problems that arise
within the realm of Federated Learning Kairouz et al. (2021).

Emerging real-world scenarios and applications have given rise to new trends
in methodologies that focus on relaxing some of FL’s core assumptions. Some of
the noteworthy examples include peer-to-peer Learning Jiang et al. (2017);
Vanhaesebrouck et al. (2017), edge-based learning Wang et al. (2019); Tran et
al. (2019), and split learning Gupta and Raskar (2018); Vepakomma et al.
(2018). In fully decentralized Peer-to-Peer (P2P) learning the objective
revolves around replacing the central orchestration server with a peer-to-peer
communications topology between the clients to mitigate the risk of single-
point of failure and address the potential communication bottlenecks in large-
scaled systems. In edge-based FL system, on the other hand, the emphasis lies
in relocating the parameter server to the vicinity of edge devices to enhance
latency and strike a balance between computation and communication. Unlike the
previous approaches, split learning introduces another dimension of
parallelism and shifts the focus from data partitioning to splitting the
execution of the learning model into layers across multiple clients and the
server.

Hierarchical approaches have consistently been instrumental in shaping the
development and structure of intricate systems, and FL is no exception to this
pattern. In Hierarchical FL approaches, the clients are clustered into groups
to reduce communication overheads and add an extra privacy layer through
periodic model updates within and across groups Rana et al. (2023). The
studies reported in Abad et al. (2020); Liu et al. (2020); Xu et al. (2021);
Hosseinalipour et al. (2022) exemplify noteworthy contributions that center on
enhancing energy consumption, latency, and communication efficiency. In Abad
et al. (2020), a two-tier hierarchical FL architecture is employed for
wireless heterogeneous cellular networks. This structure clusters mobile users
into localized cells, leveraging local FL within micro-cell clusters and
aggregating model updates at the macro-level through base stations. The
HierFAVG model introduced in Liu et al. (2020) presents a three-layer
hierarchy—client-edge-cloud—extending FL with a two-level average-based model
aggregation approach. In Xu et al. (2021), a dual-level framework is proposed,
comprising mobile devices and edge/cloud servers. This framework optimizes
resource allocation and aggregation intervals adaptively. Finally, the work
outlined in Hosseinalipour et al. (2022), unlike the proceeding models,
employs an arbitrary number of levels in its design. Their proposed tree-like
structure has a fixed hierarchical topology and is constructed based on the
physical proximity and connectivity of devices. Notably, their model utilizes
a device-to-device (D2D) consensus method to update the local models at the
lowest level, followed by a sequence of recursive aggregations through the
upstream path, culminating in a central server.

This paper presents a holonic solution for the distributed learning problem.
The multi-level architecture of holonic systems embodies an arrangement of
self-contained entities known as holons, forming a self-similar structure
rooted in the whole/part relationship. Benefitting from their balance between
local autonomy and coordinated decision-making, holonic systems adeptly tackle
challenges such as adaptability, fairness, and scalability. Additionally,
their decentralized authority and decision-making approach mitigate the risk
of a single point of failure, further enhancing their resilience. Holonic
solutions have already been employed in the design and management of machine
learning systems. As an example, the Hierarchical Agent-based Machine LEarning
PlaTform (HAMLET) presented in Esmaeili et al. (2022) leverages holons to
organize machine learning resources distributed across an array of networked
devices and orchestrates automatic parallel training and testing tasks. Our
paper, in contrast, primarily focuses on distributed learning and introduces a
novel holonic collaborative learning methodology.

Figure 1 portrays an illustrative multi-level holonic model built over a set
of connected devices through grouping them together at different levels of
granularity. The construction of such structure indeed constitutes a highly
challenges phase, potentially contingent upon various factors ranging from
physical constraints to trust and security considerations. To maintain our
main focus and streamline the analysis, we assume that the holonic
construction is preexisting and accessible. We will return to this figure in
the subsequent section as we succinctly review the fundamentals of holonic
concepts.

Figure 1. An example HoL model built over an IoT network

The presented Holonic Learning (HoL) model generalizes the conventional
distributed approaches by viewing the holons of each level, as a network of
self-contained agents (workers or clients, depending on the application) that
collaborate to optimize the objective model parameters regardless of paying
attention to the details of learning process that each agent pursues. In other
words, each holon member of the level has the autonomy to chose to learn while
being committed to the goal of the superordinate it is a member of. This model
differs from fully decentralized peer-to-peer learning models with direct
inter-agent coordination as it facilitates enforcing diverse set of policies
at different levels. On the other hand, thanks to its self-similar structure,
it is distinct from hierarchical approaches that require clear lines of
authority on one hand and relying on a central server to dictate the next
steps. The main contributions of this paper are highlighted as follows:

  * •

We present a collaborative learning framework that consolidates existing
distributed machine learning solutions. By harnessing both vertical and
horizontal interactions, and drawing inspiration from the self-similar
concepts of holonic structures, our framework enables the integration of more
sophisticated collaboration patterns while maintaining an intuitive design.

  * •

Through your empirical analysis, you provide concrete evidence of HoL’s
effectiveness in achieving consistent convergence behavior across various
holarchies and data distribution settings. The comparison of performance
metrics, such as training loss and test accuracy, demonstrates the advantages
of the holarchical structure in improving learning outcomes compared to plain
collaboration networks.

The paper is organized as follows: Section 2 reviews the fundamental
principles of holonic systems and the holonic multi-agent concepts employed in
this study. Section 3 delves into the formulation of our holonic learning
model and presents the details of a more specific use-case analyzed in this
paper. Section 4 elaborates on the experiments conducted and examines the
model’s characteristics from multiple perspectives. Finally, Section 5
concludes the paper and offers suggestions for future research directions.

##  2\. The Preamble

###  2.1. Holonic Structures

Formally introduced by Arthur Koestler, the term “holon” refers to a
philosophical concept describing entities that are simultaneously whole in
themselves and parts of larger systems Koestler (1968). In other words, a
self-contained holon embodies both individuality and integration, which helps
to model complex systems and understand how they evolve in response to
changing circumstances. Although it was originally used to describe the self-
similar hierarchical structure of biological and social system, it has found
its way into many disciplines such as computer science Esmaeili et al. (2022),
engineering Benabdellah et al. (2021), physics Capra (1984), and manufacturing
Valckenaers et al. (1997), to name a few. This section highlights the key
definitions related to the holonic systems that are heavily used in the
subsequent context.

As stated above, the key idea of holons is their dual nature, or as described
in Koestler (1978), how they behave like a two-faced Janus: the face
representing the whole looks inward and the other looks up, or outward. That
is at the same time that each holon is independent with its own
characteristics and functions, it is also a relatively more concrete component
of larger systems. Applying such a whole-part relationship to different
entities of a system create a multi-level interconnected structure of holons,
where the behavioral variations of each holon at any level is characterized by
operational variation of its subordinate holons in a lower level. In the
holonic terminology, such a hierarchical architecture is called “holarchy”,
and the higher-level and lower-level holons are respectively referred to as
“super-holons” and “sub-holons”. In this paper, we also use the terms
“terminal holons” and “non-terminal” holons to refer to holons at the lowest
level of the holarchy and the ones in the intermediate layers, respectively.

When incorporating holonic principles into the realm of computer science and
the design of distributed autonomous systems, the “Holonic Multiagent Systems
(HMAS)” emerges as an exceptionally promising and broadly applicable model to
embrace. As originally proposed in Gerber et al. (1999), a holon (holonic-
agent) is formed by a group of agents giving up part of their autonomy to the
super-holon and through the designation of one agent as the representative
called “head”. With administrative and/or directive responsibilities, the head
plans for the holon on basis of its subordinates. With that said, a holon is
formally characterized by the tuple $\left<head,subholons,C\right>$, where $C$
denote commitments, i.e., the itera-holon relationships agreed by all members
Fischer et al. (2003). Based on this formalization, each agent $a$ of a MAS
constitutes an atomic (terminal) holon
$\left<\\{a\\},\\{a\\},\emptyset\right>$ at the lowest level of the holarchy,
and there are exponential ways to form each upper-level holon based on its
subordinates.

There are multitude of research efforts, spanning both domain-specific and
generic contexts, that have delved into the dynamic formation of holarchies,
such the ones reported in Esmaeili et al. (2017); Hilaire et al. (2008);
Barbosa (2015) and Ulieru (2002) to name a few. Due to our focus on learning,
this paper opts for a predetermined, static structure for the utilized
holarchy, while acknowledging that relaxing this assumption could potentially
pave the way for novel research trajectories within the realm of distributed
learning. To emphasize both inter- and intra-holon interactions, we define a
“holarchical network” as depicted in Fig. 2. In constructing a holarchy,
super-holons are virtual concepts overseen by head holons. These head holons
can be specialized terminal holons designated solely for this role, or they
can be existing holons selected (temporarily or permenantly) to serve as heads
while concurrently fulfilling their duties as subordinates at higher levels.
To simplify the discussion in this paper, we assume the former cases in our
HoL framework.

Figure 2. An example holarchical network. The levels are numbered from top to
bottom, and different color shades represent the heterogeneity of the holons.

###  2.2. Notations

We use the following notations throughout this paper. Please note that we use
the term “holon” to refer to the nodes of the HoL model, not only to comply
with the holonic terminologies, but also due to the fact that it facilitates
referring to the both physical and abstract nodes residing at different level
of the proposed HoL model. At the physical level, one can safely assumes each
holon represent an edge device.

  * •

$\mathscr{h}^{(l)}_{i}$: the $i$-th terminal holon at level $l$.

  * •

$\bm{\mathscr{h}}^{(l)}_{i}$: the $i$-th super-holon at level $l$. This
notations is also used to denote the head of the holon.

  * •

$\mathsf{h}^{(l)}_{i}$: the $i$-th holon at level $l$ regardless of specifying
if it is terminal or non-terminal.

  * •

$n^{(l)}_{i}$: the number of subordinates in the holon
$\bm{\mathscr{h}}^{(l)}_{i}$.

  * •

$\mathcal{M}^{(l)}_{i}$: the learning model of holon $\mathscr{h}^{(l)}_{i}$.

  * •

$\bm{\theta}^{(l)}_{i}(t)$: the model parameters of holon
$\mathscr{h}^{(l)}_{i}$ at the end of iteration $t$.

  * •

$\mathcal{F}(\bm{\theta}^{(l)}_{i};\bm{x}^{(j)},y^{(j)})$: the loss function
of holon $\mathscr{h}^{(l)}_{i}$ for data sample $(\bm{x}^{(j)},y^{(j)})$. For
the sake of brevity we write it as
$\mathcal{F}^{(l)}_{i}(\bm{\theta};\bm{\xi}^{(j)})$ or
$\mathcal{F}^{(l)}_{i}(\bm{w})$ depending on the detail level.

  * •

$\mathbb{X}^{(l)}_{i}$: the private dataset used for training holon
$\mathscr{h}^{(l)}_{i}$’s model.

  * •

${C}^{(l)}_{i}$: the local update stopping criteria of holon
$\mathscr{h}^{(l)}_{i}$.

  * •

$\mathcal{G}^{(l)}_{i}$: the communication graph of super-holon
$\bm{\mathscr{h}}^{(l)}_{i}$.

  * •

$\text{Sup}(\mathsf{h}^{(l)}_{i})$: the superior(s) of holon
$\mathsf{h}^{(l)}_{i}\in\\{\mathscr{h}^{(l)}_{i},\bm{\mathscr{h}}^{(l)}_{i}\\}$
in the holarchy.

  * •

$\text{Sub}(\mathsf{h}^{(l)}_{i})$: the subordinate(s) of holon
$\mathsf{h}^{(l)}_{i}\in\\{\mathscr{h}^{(l)}_{i},\bm{\mathscr{h}}^{(l)}_{i}\\}$
in the holarchy.

  * •

$\text{Nei}(\mathsf{h}^{(l)}_{i})$: the neighbor(s) of holon
$\mathsf{h}^{(l)}_{i}\in\\{\mathscr{h}^{(l)}_{i},\bm{\mathscr{h}}^{(l)}_{i}\\}$
according to the graph $\mathcal{G}^{(l-1)}_{j}$ s.t.
$\bm{\mathscr{h}}^{(l-1)}_{j}\in\text{Sup}(\mathsf{h}^{(l)}_{i})$.

  * •

$\mathcal{F}(\bm{\theta})$: the loss function of the system-level model
$\bm{\theta}$.

  * •

$\mathcal{A}^{(l)}_{i}$: the aggregation function used by holon
$\bm{\mathsf{h}}^{(l)}_{i}$.

  * •

$[M_{1}|M_{2}]$: the column-wise concatenation of two matrices $M_{1}$ and
$M_{2}$.

Note that based on the holarchical structure of the system and the presumption
that the root (cloud) of the system is at level 0, we have:
$\bm{\theta}(t)=\bm{\theta}^{(0)}_{1}(t)$ and
$\mathcal{F}(\bm{\theta})=\mathcal{F}^{(0)}_{1}(\bm{\theta})$.

##  3\. Holonic Learning Model

This section delves into the details of the proposed distributed holonic
learning model. As stated in section 2.1 the holons of a holarchy can be seen
as reference points, representing different development stages within the
system. In other words, each level of the holarchy is self-contained and
consists of autonomous entities with specific interactions and problems,
depending on the problem’s abstraction.

###  3.1. Problem Formulation

The presented holonic learning employs a level-based learning model. In the
example depicted in Fig. 1, for instance, the learning problem on the cloud
holon is viewed as the process of optimizing model parameters based on the
information that is shared by its two subordinate collaborators without
worrying about the details of the other levels. With that said, formally, each
individual holon locally seeks the solution to the minimization of its
empirical risk function $\mathsf{F}^{(l)}_{i}$. For the terminal holon
$\mathscr{h}^{(l)}_{i}$, the problem is defined as:

$\min_{\bm{\theta}}\left[\mathsf{F}^{(l)}_{i}(\bm{\theta}):=\frac{1}{|\mathbb{X}^{(l)}_{i}|}\sum_{j=1}^{|\mathbb{X}^{(l)}_{i}|}\mathcal{F}^{(l)}_{i}(\bm{\theta};\bm{\xi}^{(j)})\right]$
(1)

and for the higher level holon $\bm{\mathscr{h}}^{(l)}_{i}$, we have:

$\min_{\bm{\theta}}\left[\mathsf{F}^{(l)}_{i}(\bm{\theta}):=\sum_{\mathscr{h}^{(l+1)}_{k}\in\text{Sub}(\bm{\mathscr{h}}^{(l)}_{i})}\phi^{(l+1)}_{k}\mathsf{F}^{(l+1)}_{k}(\bm{\theta})\right]$
(2)

where $\phi^{(l+1)}_{k}$ denote the weight the holon applies to the empirical
risk value of the subordinate holon $\mathscr{h}^{(l+1)}_{k}$.

The proposed HoL model tries to find the near-optimal solution for the above-
mentioned optimization problem in a sequence of local model updates and
aggregation of the feedback it receives from its neighbors, subordinates,
superiors. Let
$\bm{\theta}_{i}^{(l)}=\left[\theta_{i,1}^{(l)},\theta_{i,2}^{(l)},\dots,\theta_{i,d}^{(l)}\right]^{\top}\in\mathbb{R}^{d}$
specify the model parameter vector of holon $\mathsf{h}^{(l)}_{i}(t)$. For
$n$, $b$, and $u$ updates that the holon has received from its neighbors,
internal/subordinates, and superiors, respectively, we define the
corresponding matrices based the received contributions as follows:

$\bm{\Theta}_{\mathcal{N}_{i}}^{(l)}(t)=\left[\bm{\theta}_{j}^{(l)}|\ldots|\bm{\theta}_{k}^{(l)}\right]\in\mathbb{R}^{d\times
n}$ (3)

and

$\bm{\Theta}_{\mathcal{B}_{i}}^{(l)}(t)=\begin{cases}\left[\bm{\theta}_{j}^{(l+1)}|\ldots|\bm{\theta}_{m}^{(l+1)}\right]\in\mathbb{R}^{d\times
b}&\text{if }\mathsf{h}^{(l)}_{i}(t)\text{ is non-terminal}\\\ \\\
\left[\bm{\theta}_{i}^{(l)}\right]\in\mathbb{R}^{d}&\text{if
}\mathsf{h}^{(l)}_{i}(t)\text{ is terminal}\end{cases}$ (4)
$\bm{\Theta}_{\mathcal{U}_{i}}^{(l)}(t)=\begin{cases}\left[\bm{\theta}_{j}^{(l-1)}|\ldots|\bm{\theta}_{o}^{(l-1)}\right]\in\mathbb{R}^{d\times
u}&\text{if }l>0\\\ \\\ \text{none}&\text{if }l=0\end{cases}$ (5)

Additionally, each holon uses the column vectors
$\bm{w}_{\mathcal{N}_{i}}^{(l)}(t)\in\mathbb{R}^{n}$,
$\bm{w}_{\mathcal{B}_{i}}^{(l)}(t)\in\mathbb{R}^{b}$, and
$\bm{w}_{\mathcal{U}_{i}}^{(l)}(t)\in\mathbb{R}^{u}$ to weigh the
contributions of its neighbors and subordinates, respectively. With that said,
each holon uses its own aggregation function to update the value of its model
parameters for the next step as follows:

$\bm{\theta}^{(l)}_{i}(t+1)=\mathcal{A}^{(l)}_{i}\left(\bm{\Theta}_{i}^{(l)}(t),\bm{W}_{i}^{(l)}(t)\right)$
(6)

where $\bm{\Theta}_{i}^{(l)}(t)$ and $\bm{W}_{i}^{(l)}(t)$ are referred to as
contribution matrix and contribution weights vector, respectively, and are
constructed as follows:

$\bm{\Theta}_{i}^{(l)}(t):=\left[\bm{\Theta}_{\mathcal{B}_{i}}^{(l)}(t)|\bm{\Theta}_{\mathcal{N}_{i}}^{(l)}(t)|\bm{\Theta}_{\mathcal{U}_{i}}^{(l)}(t)\right]\in\mathbb{R}^{d\times(n+b+u)}$
(7)
$\bm{W}_{i}^{(l)}(t):=\left[\bm{w}_{\mathcal{B}_{i}}^{(l)}(t)^{\top}|\bm{w}_{\mathcal{N}_{i}}^{(l)}(t)^{\top}|\bm{w}_{\mathcal{U}_{i}}^{(l)}(t)^{\top}\right]^{\top}\in\mathbb{R}^{(n+b+u)}$
(8)

The proposed holonic learning employs a criteria-based iterative learning
process that dedicates two separate procedures for horizontal and vertical
contributions as presented in Algorithm 1. During its lifetime, holon
$\mathsf{h}_{i}^{l}$ uses lines 1 to 1 to continuously listen to new messages
from the peers, including its subordinates (if non-terminal), neighbors, and
superiors, in parallel and properly updates its knowledge matrices and weight
vectors accordingly. In case the holon is training in sync with its peers, it
also forwards any sent model from the super-holon to the neighbors.
Concurrently, the holon follows a set of aggregation, updating, sharing steps
until its termination criteria is satisfied. Function ReceivedEnoughPeerModels
checks whether it has received expected number of models from its peers and if
it has, it runs aggregation on the models and updates its knowledge base and
initiates its training process if it is terminal. Next, through the steps
between lines 1 and 1, the holon determines proper peers to share its train
results. Choosing the peers to send the model depends on the training stage of
the holon, its position in the holarchy, and the internal messaging policy of
the super-holon it is a member of. Fig. 3 illustrates the steps of the
algorithm in a very simple holarchy. It is worth noting that the
synchronization of the communications in this example are merely for the sake
of clarity and is not the requirement of the proposed model.

1 Procedure  _Learn(__)_

2 while  _$C_{i}^{(l)}=\text{false}$ non-blocking in parallel_ do

3
$\bm{\theta}_{j},w_{j},sender\leftarrow\textsc{Receive}(\bm{\tilde{\theta}})$

4
$\bm{\Theta}_{\mathcal{U}_{i}}^{(l)},\bm{w}_{\mathcal{U}_{i}}^{(l)},\bm{\Theta}_{\mathcal{B}_{i}}^{(l)},\bm{w}_{\mathcal{B}_{i}}^{(l)},\bm{\Theta}_{\mathcal{N}_{i}}^{(l)},\bm{w}_{\mathcal{N}_{i}}^{(l)}\leftarrow\textsc{Process}(\bm{\theta}_{j},w_{j})$

5 if  _$sender\in\text{Sup}(\mathsf{h}_{i}^{l})$ and Synchronized_ then

6 for  _$n\in\text{Nei}(\mathsf{h}_{i}^{l})$ in parallel_ do

7 $\textsc{SendModel}(n,\bm{\theta}_{j})$

8

9 end for

10

11 end if

12

13 end while

14 while  _$C_{i}^{(l)}=\text{false}$_ do

15 if  _ReceivedEnoughPeerModels $()$_ then

16
$\bm{\theta}^{(l)}_{i}(t)\leftarrow\mathcal{A}^{(l)}_{i}\left(\bm{\Theta}_{i}^{(l)},\bm{W}_{i}^{(l)}\right)$

17
$\textsc{UpdateKB}\left(\bm{\Theta}_{\mathcal{U}_{i}}^{(l)},\bm{w}_{\mathcal{U}_{i}}^{(l)},\bm{\Theta}_{\mathcal{B}_{i}}^{(l)},\bm{w}_{\mathcal{B}_{i}}^{(l)},\bm{\Theta}_{\mathcal{N}_{i}}^{(l)},\bm{w}_{\mathcal{N}_{i}}^{(l)}\right)$

18

19 else

20 Wait()

21

22 end if

23 if  _$\mathsf{h}^{(l)}_{i}$ is terminal_ then

24
$\bm{\theta}_{i}^{(l)}(t)\leftarrow\textsc{Train}(\bm{\theta}^{(l)}_{i}(t),\bm{\lambda})$

25
$\bm{\Theta}_{\mathcal{B}_{i}}^{(l)}\leftarrow\left[\bm{\theta}_{i}^{(l)}(t)\right]$

26 end if

27 $\mathcal{P}_{i}^{(l)}\leftarrow\textsc{WhoToSendModel}()$

28 for  _$p\in\mathcal{P}_{i}^{(l)}$ in parallel_ do

29 $\textsc{SendModel}(p,\theta_{i}^{(l)}(t))$

30

31 end for

32

33 end while

34

35 end

Algorithm 1 The generic learning procedure run by each holon
${\mathsf{h}}^{(l)}_{i}$ in parallel. Figure 3. An example demonstrating the
flow of communications in Algorithm 1.

The above-mentioned formulation is defined generically in such a way that it
supports the autonomy of holons and not enforcing any specific aggregation or
communication schemes, it supports the dynamic aspects of the system,
including communication graphs, memberships, peer weighting, and other
evolving factors. To demonstrate the utilization of the proposed model, we
present Holonic Averaging as a homogeneous special embodiment.

###  3.2. Holonic Averaging Learning

There is a wide spectrum of learning models that can be designed given the
proposed generic HoL framework. This section delves into the details of a
specific instance that exemplifies the principles of the presented model and
scrutinizes its corresponding algorithms. The presented Holonic Averaging
Learning (HoAL) mimics the aggregation mechanism used in the FAVG approach
McMahan et al. (2017) and provides a distributed mechanism for a supervised
learning problem. Below, we provide the breakdown of the details of HoAL’s
design components.

#### Holonic Structure

The holons in HoAL are heterogeneous with respect to their private data and
homogeneous in terms of their learning models and aggregation functions.
Additionally, we assume that each holon belongs to a unique super-holon. To
put it formally, $\forall i\neq j,\forall l,l^{\prime}\in\mathbb{N}$ we have:

$\mathcal{B}_{i}^{(l)}\cap\mathcal{B}_{j}^{(l)}=\emptyset,\quad\mathcal{A}_{i}^{(l)}=\mathcal{A}_{j}^{(l^{\prime})},\quad\mathcal{M}_{i}^{(l)}=\mathcal{M}_{j}^{(l^{\prime})}$
(9)

#### Communication Weights

As seen before, the weights that each holon assigns to the communication with
its neighbors, subordinates, and superior allow for more nuanced control over
the influence of other holons and can be dynamically defined based on metrics
such as reliability, trustworthiness, data quantity, etc. Presuming the
locality of the holons belonging the same connected graph
$\mathcal{G}_{i}^{(l)}$ and that they have similar characteristics in terms of
performance, data quality, etc., HoAL defines the inter- and intra-holon
communication weights across different levels base on the size of the training
data. Formally, let $w_{ij}^{(l)}$ denote the weight that holon
$\mathsf{h}_{i}^{(l)}$ uses for a peer
$\mathsf{h}_{j}^{(l^{\prime})}\in\mathcal{B}_{i}^{(l)}\cup\mathcal{N}_{i}^{(l)}\cup\mathcal{U}_{i}^{(l)}$.
We define $w_{ij}^{(l)}=d_{j}^{(l)}$, where:

$d_{j}=\begin{cases}|\mathbb{X}_{j}^{(l^{\prime})}|&\text{if
$\mathsf{h}_{j}^{(l^{\prime})}$ is terminal}\\\
\sum\limits_{\mathsf{h}_{k}^{(l^{\prime}+1)}\in\mathcal{B}_{j}^{(l^{\prime})}}{d_{k}^{(l^{\prime}+1)}}&\text{if
$\mathsf{h}_{j}^{(l^{\prime})}$ is non-terminal}\end{cases}$ (10)

Please note that despite its recursive definition, the weight does not need to
be calculated recursively each time it is needed. Each non-terminal holon
keeps track of its training size whenever it receives updates from the
subordinates.

#### Updates Stopping Criteria

Each holon $\mathsf{h}_{i}^{(l)}$ in HoAL employs a numerical threshold
$\mathsf{k}_{i}^{(l)}$ to control the number of its local updates. To put it
formally, we have:

$C_{i}^{(l)}:=t_{i}^{(l)}>\mathsf{k}_{i}^{(l)}$ (11)

#### Aggregation

As stated before, HoAL employs weighted averaging in each holon to aggregate
the values of model parameters collected in its contribution matrix. Strictly
speaking, we have:

$\mathcal{A}^{(l)}_{i}\left(\bm{\Theta}_{i}^{(l)},\bm{W}_{i}^{(l)}\right):=\frac{1}{{\bm{W}_{i}^{(l)}}^{\top}\bm{1}_{z}}\left(\bm{\Theta}_{i}^{(l)}\bm{W}_{i}^{(l)}\right)$
(12)

where $z$ denotes the size of vector $\bm{W}_{i}^{(l)}$ at the time of
aggregation, and $\bm{1}_{z}$ represents a $z$-dimensional column vector of
ones. In this paper, we utilize mini-batch SGD as the optimizer and assume
each training iteration equal to a specific number of epochs at the terminal
level.

###  3.3. Existing models as special cases

The strength of holonic learning, in its original definition as presented in
section 3.1, is its flexibility to encapsulate an expansive spectrum of
distributed learning scenarios. An example that showcases its versatility is
when we further investigate the structure of the HAL. While being a specific
instantiation of the HoL with some limitations, HoAL still extends the
flexibility of several existing distributed approaches, such as federated
averaging, hierarchical federated averaging, and peer-to-peer learning to name
a few. Figure 4 visualizes a high-level example HoAL model accompanied by its
particular embodiments to support the aforementioned standard models. Assuming
SGD as the optimization method employed by the terminal holons, the depicted
instantations can further be discussed as follows:

(a) HoAL model

(b) HFL model

(c) P2P model

(d) FL model

Figure 4. The holonic averaging learning model together with different special
cases to implement distributed models of federated averaging, hierarchical
federated averaging, and peer-to-peer models. Different shades in the lower
half of the terminal holons represent the difference in their local data.

#### Federated Averaging

: An FL system can be modeled by a one-level HoAL, where the server is modeled
by the root holon $\bm{\mathscr{h}}_{1}^{(0)}$ and each of the clients with
terminal holon ${\mathscr{h}}_{i}^{(1)},i=1,\dots,n_{1}^{(0)}$. Additionally,
$\mathcal{G}_{1}^{(0)}=\overline{K}_{n}$, where $n$ and $\overline{K}_{n}$
denote the number of clients and a null graph of size $n$. As depicted in
Figure 4(d), all clients use their own private data to learn the same ML
model. At each update carried out by holons $\bm{\mathscr{h}}_{1}^{(0)}$ and
${\mathscr{h}}_{i}^{(1)}$ the contribution matrix and weights vectors would
respectively be:

$\bm{\Theta}_{1}^{(0)}=\bm{\Theta}_{\mathcal{B}_{1}}^{(0)}=\left[\bm{\theta}_{1}^{(1)}|\dots|\bm{\theta}_{n_{1}^{(0)}}^{(1)}\right]$
(13)
$\bm{W}_{1}^{(0)}=\bm{w}_{\mathcal{B}_{1}}^{(0)}=\left[|\mathbb{X}_{1}^{(1)}|,\dots,|\mathbb{X}_{n_{1}^{(0)}}^{(1)}|\right]^{\top}$
(14)
$\bm{\Theta}_{i}^{(1)}=\bm{\theta}_{i}^{(1)}=\bm{\theta}_{i}^{(1)}-\eta_{i}^{(1)}\left[\frac{1}{|\mathbb{X}^{(1)}_{i}|}\sum_{j=1}^{|\mathbb{X}^{(1)}_{i}|}\nabla\mathcal{F}^{(1)}_{i}(\bm{\theta};\bm{\xi}^{(j)})\right]$
(15) $\bm{W}_{i}^{(1)}=\left[|\mathbb{X}_{i}^{(1)}|\right]^{\top}$ (16)

#### Peer-to-Peer Learning

: Similar to Fl, a P2P learning model can also be modeled by a one-level HoAL,
as depicted in Figure 4(c). A key distinction her is that the communication
graph connecting all the terminal holons, i.e. $\mathcal{G}_{1}^{(0)}$ is
connected. Additionally, as there is no need for server-side aggregation of
the learned models, the root $\bm{\mathscr{h}}_{1}^{(0)}$, being connected to
at least one client, only serves as the initiator of the learning process.
Once initiated, the model updates of each holon ${\mathscr{h}}_{i}^{(1)}$ will
be calculated using the following contribution matrix and weight vector
configurations:

$\begin{split}\bm{\Theta}_{i}^{(1)}=\left[\bm{\theta}_{i}^{(1)}|\bm{\Theta}_{\mathcal{N}_{i}}^{(1)}\right]=\Biggl{[}\bm{\theta}_{i}^{(1)}&-\frac{\eta_{i}^{(1)}}{|\mathbb{X}^{(1)}_{i}|}\sum_{j=1}^{|\mathbb{X}^{(1)}_{i}|}\nabla\mathcal{F}^{(1)}_{i}(\bm{\theta};\bm{\xi}^{(j)})|\\\
\quad\bm{\theta}_{e}^{(1)}&-\frac{\eta_{e}^{(1)}}{|\mathbb{X}^{(1)}_{e}|}\sum_{j=1}^{|\mathbb{X}^{(1)}_{e}|}\nabla\mathcal{F}^{(1)}_{e}(\bm{\theta};\bm{\xi}^{(j)})|\dots|\\\
\bm{\theta}_{f}^{(1)}&-\frac{\eta_{f}^{(1)}}{|\mathbb{X}^{(1)}_{f}|}\sum_{j=1}^{|\mathbb{X}^{(1)}_{f}|}\nabla\mathcal{F}^{(1)}_{f}(\bm{\theta};\bm{\xi}^{(j)})\Biggr{]}\end{split}$
(17)
$\bm{W}_{i}^{(1)}=\left[|\mathbb{X}_{i}^{(1)}|,\bm{w}_{\mathcal{N}_{i}}^{(1)}\right]^{\top}=\left[|\mathbb{X}_{i}^{(1)}|,|\mathbb{X}_{e}^{(1)}|,\dots,|\mathbb{X}_{f}^{(1)}|\right]^{\top}$
(18)

#### Hierarchical Federated Learning

: The client-edge-server hierarchical federated learning can be modeled by a
two level balanced holarchy, in which the root holon
$\bm{\mathscr{h}}_{1}^{(0)}$ serves as the cloud, the non-terminal holons
$\bm{\mathscr{h}}_{i}^{(1)},i=1,\dots,n_{1}^{(0)}$ serves as the edge model-
aggregators, and finally the terminal holon
${\mathscr{h}}_{j}^{(2)},j=1,\dots,n_{i}^{(1)}$ represent the clients employed
to learn optimize the models using their own data. Additionally, all the
communication graphs of the holons are null, i.e.
$\mathcal{G}_{1}^{(0)}=\overline{K}_{n_{1}^{(0)}}$ and
$\mathcal{G}_{i}^{(1)}=\overline{K}_{n_{i}^{(1)}},i=1,\dots,n_{i}^{(1)}$.
Similar to the above-mentioned federated averaging case, the contribution
matrices and the weight vectors would be as follows:

$\begin{split}\bm{\Theta}_{1}^{(0)}=\bm{\Theta}_{\mathcal{B}_{1}}^{(0)}=\left[\bm{\theta}_{1}^{(1)}|\dots|\bm{\theta}_{n_{1}^{(0)}}^{(1)}\right]=&\left[\mathcal{A}_{1}^{(1)}\left(\bm{\Theta}_{1}^{(1)},\bm{W}_{1}^{(1)}\right),\dots,\right.\\\
&\left.\quad\mathcal{A}_{n_{1}^{(0)}}^{(1)}\left(\bm{\Theta}_{n_{1}^{(0)}}^{(1)},\bm{W}_{n_{1}^{(0)}}^{(1)}\right)\right],\quad\end{split}$
(19)
$\bm{W}_{1}^{(0)}=\bm{w}_{\mathcal{B}_{1}}^{(0)}=\left[\sum_{\mathscr{h}_{j}^{(2)}\in\mathcal{B}_{1}^{(1)}}|\mathbb{X}_{j}^{(2)}|,\dots,\sum_{\mathscr{h}_{j}^{(2)}\in\mathcal{B}_{n_{1}^{(0)}}^{(1)}}|\mathbb{X}_{j}^{(2)}|\right]^{\top}$
(20)
$\bm{\Theta}_{i}^{(1)}=\bm{\Theta}_{\mathcal{B}_{i}}^{(1)}=\left[\bm{\theta}_{1}^{(2)}|\dots|\bm{\theta}_{n_{i}^{(1)}}^{(2)}\right]$
(21)
$\bm{W}_{i}^{(1)}=\bm{w}_{\mathcal{B}_{i}}^{(1)}=\left[|\mathbb{X}_{1}^{(2)}|,\dots,|\mathbb{X}_{n_{i}^{(1)}}^{(2)}|\right]^{\top}$
(22)
$\bm{\Theta}_{i}^{(2)}=\bm{\theta}_{i}^{(2)}=\bm{\theta}_{i}^{(2)}-\eta_{i}^{(2)}\left[\frac{1}{|\mathbb{X}^{(2)}_{i}|}\sum_{j=1}^{|\mathbb{X}^{(2)}_{i}|}\nabla\mathcal{F}^{(2)}_{i}(\bm{\theta};\bm{\xi}^{(j)})\right]$
(23) $\bm{W}_{i}^{(2)}=\left[|\mathbb{X}_{i}^{(2)}|\right]^{\top}$ (24)

##  4\. Experiments

This section details our experiments investigating the convergence and
performance of HoAL under various design scenarios and data distributions. The
HoAL framework was implemented using `pytorch` along with its `multiprocess`
and `distributed` libraries. We have utilized training loss and test accuracy
as the measure of learning convergence and performance.

###  4.1. Settings

The experimental setup comprised 10 terminal holons (agents) organized into
four distinct holarchical structures (see Figure 6). These structures were
chosen based on their shared characteristic of utilizing a collaboration
network among terminal holons, representing various facets of collaborations.
For instance, HoAL1P resembles a P2P learning setting, primarily relying on
inter-holon communications. In contrast, the remaining three structures differ
in the number of holarchical levels while confining collaborations to peers
within the holons at each level.

Throughout the experiments, personal data for terminal holons and the
collaboration network at the terminal level remained constant. We utilized the
standard `MNIST` handwritten digits classification dataset, with all terminal
holons equipped with a Convolutional Neural Netwrok Model as in McMahan et al.
(2017) initialized with the same values for its 21,840 trainable parameters.
Terminal holons used mini-batch Stochastic Gradient Descent (SGD) with a batch
size of 32 and an initial learning rate of 0.01. Both terminal and non-
terminal holons use separate CPU processes to train/aggregate their models in
parallel and employ threads to handle inter-holon communications.

(a) HoAL1P

(b) HoAL2L

(c) HoAL3L

(d) HoAL4L

Figure 5. The different holarchical HoAL structures used for experiments.

Data distribution, a pivotal factor in distributed machine learning, was
evaluated in three settings: `IID`, `EqNIID`, and `UEqNIID`. In `EqNIID`, each
terminal holon possessed 500 samples from 2 class labels, while in `UEqNIID`,
terminals had varying numbers of samples from sorted groups.

Collaboration schemes varied across holarchy levels: each terminal holon
underwent 250 rounds of a 5-epoch budget, collaborating with peers and
superordinates at the end of each round. Non-terminal holons, on the other
hand, had a 2-round local budget for aggregations and guiding subordinates
before communicating with their own superordinates. Non-terminal holons budget
is reset after each update from their superordinates. To make comparisons
fair, we made the duration of experiments be controlled by terminal holons.
That is, each terminal holon quits as soon as its budget of 200 rounds is
exhausted, and the non-terminal holons quit once all their subordinates
finished their job.

###  4.2. Results

In conducting an empirical analysis of HoL’s convergence, we observed
consistent behavior across various holarchies and data distribution settings.
Figures 6(a) and 6(b) present a comparative view of the training loss and test
accuracy, respectively, for the experimented holarchies on `IID` data
distribution. These metrics are averaged over all terminal holons. As seen,
the holarchical structure has yielded relatively improved performance, in both
training loss and test accuracy, compared to the plain usage of collaboration
network among terminal holons as represented by HoAL1P model.

This trend persists under the `UEqNIID` data distribution, as demonstrated in
Figures 7(a) and 7(b). Here, the terminal holons exhibit a similar pattern of
enhanced performance compared to the collaboration network model, though with
comparatively smaller improvements than the `IID` case. The nuances observed
in these results indicate the insensitivity of the models to unbalanced data
distribution. Such resilience to uneven data distributions further underscores
the potential robustness of the HoAL learnings framework.

(a) Training Loss

(b) Test Accuracy

Figure 6. The average performance of all terminal holons of different
holarchical structures trained on MNIST IID distribution..

(a) Training Loss

(b) Test Accuracy

Figure 7. The average performance of all terminal holons of different
holarchical structures trained on MNIST UENIID distribution..

In the `EqNIID` setting, illustrated in Figures 8(a) and 8(b), we persistently
observe similar behavior, albeit with more noticeable distinctions in the
performance of the examined structures. Unlike the `IID` setting, the
introduction of more holarchical levels, as in HoAL4L and HoAL3L cases, has
led to a degradation in the generalization capability of terminal holons on
average. This observed phenomenon can be primarily attributed to the
collaboration delays enforced by extra levels of the holarchy.

(a) Training Loss

(b) Test Accuracy

Figure 8. The average performance of all terminal holons of different
holarchical structures trained on MNIST EqNIID distribution..

While we acknowledge that the presented experimental settings are on a small
scale, limiting our ability to fully characterize the extensive potential of
the introduced HoL framework, the observed trends and insights serve as a
foundational exploration. These initial findings provide a compelling basis
for further, larger-scale investigations on the capabilities and robustness of
the HoL framework in diverse and complex real-world scenarios.

##  5\. Conclusion

Inspired by the self-similar holonic paradigm, this paper introduced Holonic
Learnings (HoL) as a flexible agent-based distributed machine learning
framework. Within HoL, holons are self-contained entities, responsible for
managing their individual datasets and model aggregation processes, and its
unique characteristic lies in the collaborative synergy across the holarchy,
fostering both horizontal and vertical collaborations among holons. This
holistic approach not only allows each holon to maintain autonomy over its
local data and aggregation strategies but also facilitates the integration of
diverse set of policies and learning strategies at different segments and
levels of the system. In our exploration of HoL, we presented HoAL, as a
homogenous instantiation, featuring a weighted parameter averaging aggregation
strategy and showed how it extends the flexibility of existing distributed
approaches such as federated learning. The reported investigation of HoAL’s
capability in various design and data distribution settings revealed its
convergence and effectiveness, providing insights into the applicability of
HoL.

As we delve into the future of distributed machine learning, Holonic Learnings
stands out as a promising model that harmonizes individual agency with
collaborative intelligence and, further leveraging the characteristics of
holonic multiagent systems, offers vast potentials for the development of
scalable, adaptable, privacy-focused, and heterogenous learning environments.
This paper used empirical results to scrutinize HoAL’s behavior and
demonstrated the pivotal role of holarchical design in convergence rate and
generalization, even in the presented small-scale settings. In contemplating
the next steps, conducting theoretical investigation of HoL’s convergence in
both convex and non-convex optimization settings and integrating the gained
insight into the design of large-scale holarchies are crucial for its real-
world applications.

## References

  * (1)
  * Abad et al. (2020) Mehdi Salehi Heydar Abad, Emre Ozfatura, Deniz Gunduz, and Ozgur Ercetin. 2020.  Hierarchical federated learning across heterogeneous cellular networks. In _ICASSP 2020-2020 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP)_. IEEE, 8866–8870.
  * Barbosa (2015) José Barbosa. 2015\.  _Self-organized and evolvable holonic architecture for manufacturing control_.  Ph.D. Dissertation. Université de Valenciennes et du Hainaut-Cambresis.
  * Benabdellah et al. (2021) Abla Chaouni Benabdellah, Imane Bouhaddou, and Asmaa Benghabrit. 2021\.  SmartDFRelevance: a holonic agent based system for engineering industrial projects in concurrent engineering context. In _Artificial Intelligence and Industrial Applications: Smart Operation Management_. Springer, 103–123.
  * Capra (1984) Fritjof Capra. 1984\.  _The turning point: Science, society, and the rising culture_.  Bantam.
  * Esmaeili et al. (2022) Ahmad Esmaeili, John C. Gallagher, John A. Springer, and Eric T. Matson. 2022.  HAMLET: A Hierarchical Agent-Based Machine Learning Platform.  _ACM Trans. Auton. Adapt. Syst._ 16, 3–4, Article 9 (jul 2022), 46 pages.  https://doi.org/10.1145/3530191
  * Esmaeili et al. (2017) Ahmad Esmaeili, Nasser Mozayani, Mohammad Reza Jahed Motlagh, and Eric T. Matson. 2017.  A Socially-Based Distributed Self-Organizing Algorithm for Holonic Multi-Agent Systems.  _Cogn. Syst. Res._ 43, C (2017), 21–44.  https://doi.org/10.1016/j.cogsys.2016.12.001
  * Fischer et al. (2003) Klaus Fischer, Michael Schillo, and Jörg Siekmann. 2003\.  Holonic multiagent systems: A foundation for the organisation of multiagent systems. In _International conference on industrial applications of holonic and multi-agent systems_. Springer, 71–80.
  * Gerber et al. (1999) Christian Gerber, Jörg Siekmann, and Gero Vierke. 1999\.  _Holonic multi-agent systems_.  Technical Report. German Research Center for Artificial Intelligence (DFKI).
  * Gupta and Raskar (2018) Otkrist Gupta and Ramesh Raskar. 2018.  Distributed learning of deep neural network over multiple agents.  _Journal of Network and Computer Applications_ 116 (2018), 1–8.
  * Hilaire et al. (2008) Vincent Hilaire, Abder Koukam, and Sebastian Rodriguez. 2008\.  An Adaptative Agent Architecture for Holonic Multi-Agent Systems.  _ACM Trans. Auton. Adapt. Syst._ 3, 1, Article 2 (mar 2008), 24 pages.  https://doi.org/10.1145/1342171.1342173
  * Hosseinalipour et al. (2022) Seyyedali Hosseinalipour, Sheikh Shams Azam, Christopher G Brinton, Nicolo Michelusi, Vaneet Aggarwal, David J Love, and Huaiyu Dai. 2022\.  Multi-stage hybrid federated learning over large-scale D2D-enabled fog networks.  _IEEE/ACM transactions on networking_ 30, 4 (2022), 1569–1584.
  * Jiang et al. (2017) Zhanhong Jiang, Aditya Balu, Chinmay Hegde, and Soumik Sarkar. 2017.  Collaborative Deep Learning in Fixed Topology Networks. In _Proceedings of the 31st International Conference on Neural Information Processing Systems_ (Long Beach, California, USA) _(NIPS’17)_. Curran Associates Inc., Red Hook, NY, USA, 5906–5916.
  * Kairouz et al. (2021) Peter Kairouz, H Brendan McMahan, Brendan Avent, Aurélien Bellet, Mehdi Bennis, Arjun Nitin Bhagoji, Kallista Bonawitz, Zachary Charles, Graham Cormode, Rachel Cummings, et al. 2021\.  Advances and open problems in federated learning.  _Foundations and Trends® in Machine Learning_ 14, 1–2 (2021), 1–210.
  * Koestler (1968) Arthur Koestler. 1968\.  _The ghost in the machine._ Macmillan.
  * Koestler (1978) Arthur Koestler. 1978\.  _Janus: A Summing Up._ Random House.
  * Liu et al. (2020) Lumin Liu, Jun Zhang, SH Song, and Khaled B Letaief. 2020\.  Client-edge-cloud hierarchical federated learning. In _ICC 2020-2020 IEEE International Conference on Communications (ICC)_. IEEE, 1–6.
  * McMahan et al. (2017) Brendan McMahan, Eider Moore, Daniel Ramage, Seth Hampson, and Blaise Aguera y Arcas. 2017\.  Communication-efficient learning of deep networks from decentralized data. In _Artificial intelligence and statistics_. PMLR, 1273–1282.
  * Rana et al. (2023) Omer Rana, Theodoros Spyridopoulos, Nathaniel Hudson, Matt Baughman, Kyle Chard, Ian Foster, and Aftab Khan. 2023.  Hierarchical and Decentralised Federated Learning.  _arXiv preprint arXiv:2304.14982_ (2023).
  * Tran et al. (2019) Nguyen H Tran, Wei Bao, Albert Zomaya, Minh NH Nguyen, and Choong Seon Hong. 2019.  Federated learning over wireless networks: Optimization model design and analysis. In _IEEE INFOCOM 2019-IEEE conference on computer communications_. IEEE, 1387–1395.
  * Ulieru (2002) Mihaela Ulieru. 2002\.  Emergence of holonic enterprises from multi-agent systems: A fuzzy-evolutionary approach.  _Soft Computing Agents: A new perspective on dynamic information systems_ (2002), 187–215.
  * Valckenaers et al. (1997) PAUL Valckenaers, Hendrik Van Brussel, Luc Bongaerts, and Jozef Wyns. 1997.  Holonic manufacturing systems.  _Integrated Computer-Aided Engineering_ 4, 3 (1997), 191–201.
  * Vanhaesebrouck et al. (2017) Paul Vanhaesebrouck, Aurélien Bellet, and Marc Tommasi. 2017\.  Decentralized Collaborative Learning of Personalized Models over Networks. In _Proceedings of the 20th International Conference on Artificial Intelligence and Statistics_ _(Proceedings of Machine Learning Research, Vol. 54)_ , Aarti Singh and Jerry Zhu (Eds.). PMLR, 509–517.  https://proceedings.mlr.press/v54/vanhaesebrouck17a.html
  * Vepakomma et al. (2018) Praneeth Vepakomma, Otkrist Gupta, Tristan Swedish, and Ramesh Raskar. 2018.  Split learning for health: Distributed deep learning without sharing raw patient data.  _arXiv preprint arXiv:1812.00564_ (2018).
  * Wang et al. (2019) Shiqiang Wang, Tiffany Tuor, Theodoros Salonidis, Kin K Leung, Christian Makaya, Ting He, and Kevin Chan. 2019.  Adaptive federated learning in resource constrained edge computing systems.  _IEEE journal on selected areas in communications_ 37, 6 (2019), 1205–1221.
  * Xu et al. (2021) Bo Xu, Wenchao Xia, Wanli Wen, Pei Liu, Haitao Zhao, and Hongbo Zhu. 2021\.  Adaptive hierarchical federated learning over wireless networks.  _IEEE Transactions on Vehicular Technology_ 71, 2 (2021), 2070–2083.
  * Yang et al. (2019) Qiang Yang, Yang Liu, Tianjian Chen, and Yongxin Tong. 2019\.  Federated machine learning: Concept and applications.  _ACM Transactions on Intelligent Systems and Technology (TIST)_ 10, 2 (2019), 1–19.