# Asynchronous Multi-Server Federated Learning for Geo-Distributed Clients

Yuncong Zuo<EMAIL_ADDRESS>Delft University of
TechnologyDelftNetherlands ,  Bart Cox<EMAIL_ADDRESS>0000-0001-5209-6161
Delft University of TechnologyDelftNetherlands ,  Lydia Y. Chen
<EMAIL_ADDRESS>0000-0002-4228-6735 Delft University of
TechnologyDelftNetherlands and  Jérémie Decouchant<EMAIL_ADDRESS>0000-0001-9143-3984 Delft University of TechnologyDelftNetherlands

(2018)

###### Abstract.

Federated learning (FL) systems enable multiple clients to train a machine
learning model iteratively through synchronously exchanging the intermediate
model weights with a single server. The scalability of such FL systems can be
limited by two factors: server idle time due to synchronous communication and
the risk of a single server becoming the bottleneck. In this paper, we propose
a new FL architecture, to our knowledge, the first multi-server FL system that
is entirely asynchronous, and therefore addresses these two limitations
simultaneously. Our solution keeps both servers and clients continuously
active. As in previous multi-server methods, clients interact solely with
their nearest server, ensuring efficient update integration into the model.
Differently, however, servers also periodically update each other
asynchronously, and never postpone interactions with clients. We compare our
solution to three representative baselines – FedAvg, FedAsync and HierFAVG –
on the MNIST and CIFAR-10 image classification datasets and on the WikiText-2
language modeling dataset. Our solution converges to similar or higher
accuracy levels than previous baselines and requires 61% less time to do so in
geo-distributed settings.

Byzantine Learning, Asynchronous Learning, Resource Heterogeneity

††copyright: acmlicensed††journalyear: 2018††doi: XXXXXXX.XXXXXXX††conference:
Make sure to enter the correct conference title from your rights confirmation
emai; June 03–05, 2018; Woodstock, NY††isbn: 978-1-4503-XXXX-X/18/06

##  1\. Introduction

Federated Learning (FL) is the emerging paradigm to train machine learning
models iteratively on extensive data that are dispersed across many users and
cannot be openly shared due to privacy concerns or regulations (McMahan et
al., 2017). The FL paradigm is widely adopted in applications such as mobile
keyboard prediction (Hard et al., 2018), spoken language comprehension (Huang
et al., 2020), and digital healthcare (Li et al., 2020a; Xu et al., 2021;
Brisimi et al., 2018; Rieke et al., 2020). The most prevalent FL architecture
is composed of a single-server and multiple clients, and assumes a synchronous
training process over multiple rounds (McMahan et al., 2017). In each round,
the server selects a set of clients that receive the current global model,
train it with their local dataset, and send a model update back to the server.
Once it has received all model updates, the server aggregates them and
computes the new global model that is used in the subsequent round. The
training process is repeated until the global model converges.

Figure 1. A training round of a synchronous FL system with heterogeneous devices and networks. The times at which the server receives client updates are indicated with colored spheres. The server waits for all client updates to be received to update the model, which results in a mostly idle server. Table 1. Comparison of representative multi-server FL algorithms with Spyker. gray!25 | Mutli-Server | Async Client/Server | Async Server/Server
---|---|---|---
FedAvg (McMahan et al., 2017) | ✗ | ✗ | N/A
FedAsync (Xie et al., 2019) | ✗ | ✓ | N/A
HierFAVG (Abad et al., 2020) | ✓ | ✓ | ✗
Spyker | ✓ | ✓ | ✓
Sync-Spyker | ✓ | ✓ | ✗

The duration of a training round in a synchronous FL framework is determined
by the longest client update time, because a round only completes when all
updates of selected clients are received. Clients exhibiting heterogeneous
computational capacity and network bandwidth can unfortunately prolong the
update time, degrading the efficiency of FL systems. The computational
heterogeneity stems from the difference in their hardware and software stacks;
whereas network latency difference is rooted in the geo-distributed settings
where clients are spread across several countries or continents (Cano et al.,
2016; Hsieh et al., 2017; Zhu et al., 2022). The result of such heterogeneity
is often system inefficiency: the server, along with clients near it or those
with superior computing power, remain underutilized. Fig. 1 illustrates this
problem with a synchronous FL round that involves five clients distributed
over two continents. The server has to collect all five client updates before
updating the global model. The first four clients become idle after they send
an update to the server and the server remains idle until the end of the
round.

Asynchronous FL frameworks (Xie et al., 2019; Wang et al., 2022; Chen et al.,
2020) aim at eliminating performance bottlenecks arising from heterogeneity in
client’s computation and networking. Under such frameworks, the server needs
to immediately process a client’s model update upon receiving it and return
the update global model to the corresponding client. Thus, slow clients in an
asynchronous FL framework no longer hamper the overall performance. However,
this can inadvertently shift the bottleneck to the server. Specifically, if
the server is continuously busy processing client updates, it might cause
waiting periods for clients, in particular as the number of clients increases.

To prevent potential server-based bottlenecks, one can deploy multiple servers
that each interact with a subset of the clients. In particular, hierarchical
FL frameworks (Liu et al., 2020; Abad et al., 2020; Li et al., 2019) rely on
secondary servers that interact with clients, and on a primary server that
merges the models of the secondary servers. This multi-server structure
decreases the average communication latency between a client and a server,
given that clients upload their model updates to the nearest server. However,
in such hierarchical FL frameworks, the primary server still relies on a
synchronous procedure to update its global model, which in geo-distributed
settings implies that secondary servers stop replying to clients.

We introduce Spyker 111The Spyker was the first four-wheel-drive car. The
wheels of such a car are similar to our FL framework’s servers: they make
progress asynchronously and independently from each other, and occasionally
synchronize., a novel multi-server and fully asynchronous FL framework that
relies on a flat architecture of geo-distributed servers. Moreover,
interactions between clients and their nearest server, and between servers
themselves are asynchronous to maintain low response times, and reduce both
client and server idle times. Spyker accelerates the training of a model in
geo-distributed settings, where heterogeneous clients and servers are located
in different geographical regions.

Table 1 compares Spyker to the relevant related work and highlight its main
characteristics

Contributions. This work makes the following contributions to realize the
asynchronous multi-server FL vision:

  * •

To support multi-server asynchronous FL, we define the age of a server model
and the staleness of a client update. Based on these definitions, a server
determines the weight that server models or client updates should be given
when aggregating them into its local model.

  * •

To maintain high accuracy despite resource heterogeneity, Spyker triggers an
exchange and aggregation of server updates when the ages of two server models
differ too much or when the age of a model is sufficiently changed since the
last server model aggregation. To minimize the complexity of the server model
aggregation process, in particular in asynchronous networks, a single server
at a time can trigger the exchange of server models. This decision is helped
by a token that servers circulate among themselves and that collects the age
of all server models.

  * •

To further maintain accuracy despite data heterogeneity, i.e., clients having
non-independent and identically distributed (non-iid) datasets, we detail a
local learning rate decay strategy that prevents server models from becoming
biased towards the data distributions of fast clients.

  * •

We evaluate Spyker’s performance and compare it to three recent and
representative FL frameworks, namely FedAvg (McMahan et al., 2017), HierFAVG
(Abad et al., 2020) and FedAsync (Xie et al., 2019). We also consider a
variant of Spyker, which we call Sync-Spyker, where servers use a synchronous
model exchange protocol. We evaluate all frameworks in emulated geo-
distributed settings using the MNIST and CIFAR-10 image datasets and the
WikiText2 language modeling dataset. We report the accuracy and convergence
speed of all frameworks, and evaluate their ability to scale with the numbers
of clients and servers.

Our results show that Spyker converges to similar or higher accuracy levels
than previous baselines with a 61% shorter running time in geo-distributed
settings. Spyker also scales better than other baselines with the number of
clients thanks to its flat multi-server architecture. For example, increasing
the number of clients from 100 to 200 multiplies the time our three baselines
require to reach 90% accuracy by at least 1.64, while Spyker’s convergence
time is only multiplied by 1.21 (i.e., scalability is improved by 26%).

This paper is organized as follows. Sec. 2 provides some background on FL
systems and introduces the most commonly used synchronous, asynchronous, and
hierarchical FL algorithms. Sec. 3 describes our system model and gives an
overview of Spyker. Sec. 4 discusses Spyker in detail, including the local
update and global model exchange strategies. Sec. 5 presents the results of
our performance evaluation. Sec. 6 discusses the related work, and Sec. 7
concludes this paper.

##  2\. Background

Table 2. Notations and parameter values

Symbol |  |  Description
---|---|---
$D_{k}$ |  |  Local dataset of client $k$
$d_{k}$ |  |  Number of data points in $D_{k}$
$d$ |  |  Total number of data points $\sum_{k}d_{k}$
$n$ |  |  Number of servers
$n_{C}$ |  |  Number of clients
$W^{t}_{k}$ |  |  Weight vector of a model at client $k$
$\eta_{k}$ |  |  Learning rate of client/server $k$
$T_{k}$ |  |  Number of local epochs on client $k$
$A_{k}$ |  |  Model age of client/server $k$
Symbol | Value |  Description
$\eta_{init}$ | 0.5 |  Initial client learning rate
$\eta_{min}$ | $10^{-6}$ |  Minimum client learning rate
$\beta$ | 0.05 |  Decaying rate
$h_{\mathrm{inter}}$ | $n_{C}/5n$ |  Age drift threshold between server models
$h_{\mathrm{intra}}$ | 350 |  Age drift threshold since last global aggregation
$\phi$ | 1.5 |  Activation rate
$\eta_{\alpha}$ | 0.6 |  Aggregation rate (server side)

(a) Centralized

(b) Hierarchical Multi-Server

(c) Flat Multi-Server (used in this work)

Figure 2. Architectures of FL systems

This section provides necessary background on synchronous, asynchronous and
multi-server Federated Learning. Fig. 2 illustrates the architectures these
paradigms rely on, and Tab. 2 summarizes the notations we use throughout the
paper.

###  2.1. Synchronous Federated Learning

The first FL framework was proposed by McMahan et al. in 2017 (McMahan et al.,
2017). It relies on a classical centralized architecture where a server
interacts with all clients (2(a)). The first FL algorithms, namely FedAvg and
FedSgd, are synchronous. In each round, the server selects a group of clients
and sends them the latest global model. Upon receiving a global model $W^{t}$,
a client $k$ trains it on its local dataset $D_{k}$ to obtain $W_{k}^{t+1}$
using Eq. 1, where $F$ is the loss function of the classification, and
$\eta_{k}$ is the learning rate.

(1) $W_{k}^{t+1}=W_{k}^{t}-\eta_{k}\frac{\partial F(W^{t},D_{k})}{\partial
W^{t}},$

Afterwards, client $k$ sends its model update $W_{k}^{t+1}$ to the server.
Upon receiving all client updates it expects in a given round, the server
aggregates them to compute the next global model $W^{t+1}$ using Eq. 2. In
this equation, $d_{k}$ is the number of data points of client $k$ and $d$ is
the total number of data points across all clients.

(2) $W^{t+1}=\sum_{k=1}^{K}{\frac{d_{k}}{d}W_{k}^{t}}.$

The ratio $\frac{d_{k}}{d}$ is used to weight the participation of a client
during aggregation. At the beginning of the following round, the server sends
$W^{t+1}$ to the next selected clients. The training process is repeated over
multiple rounds until the global model converges. In practical settings, where
clients have different computing speeds and are connected to the server with
heterogeneous network links, synchronous FL algorithms may lead to the server
and some clients staying idle. For example, one client in 2(a) would
communicate with the server within only 37 ms. This client and the server
would have to stay idle waiting for the clients with 435 ms and 510 ms
communication latency with the server to send their update.

###  2.2. Asynchronous Federated Learning

Asynchronous FL systems, such as FedAsync (Xie et al., 2019), accelerate
convergence in heterogeneous networks and under client heterogeneity. To do
so, the server aggregates an update into a global model immediately after it
receives it, and returns the new global model to the client. For example, a
server using the asynchronous version of FedAvg in classical centralized
settings (2(a)) would update the global model using Eq. 3. In this equation,
$W^{t}$ is the $t$-th global model, $W^{k}_{t}$ is the local update client $k$
computed on $W^{t}$, $\frac{d_{k}}{d}$ is the data proportion of the update,
and $s(\tau)$ is a staleness parameter that dampens the effect of client $k$’s
update if it relates to a model that is older than the one the server
currently has.

(3) $W^{t+1}=W^{t}-s(\tau)\frac{d_{k}}{d}\left(W^{t}_{k}-W^{t+1}_{k}\right)$

An asynchronous FL scheme aims at keeping the clients busy training a model on
their datasets, and the server updating its model more frequently, both of
which are expected to speed up the model convergence. However, asynchronous
schemes also increase the computational and communication loads on the server,
which might become a performance bottleneck if the number of clients is too
high.

###  2.3. Multi-server Federated Learning

Multi-server FL systems rely on multiple servers that interact with disjoint
subsets of the clients to limit their computational and communication loads.
Hierarchical FL systems (2(b)), such as HierFAVG (Abad et al., 2020), are
multi-server and rely on a principal server to maintain the global model (Abad
et al., 2020; Liu et al., 2020; Wang and Wang, 2022; Li et al., 2019). In
those systems, a model is updated in two times in each round. The first
aggregation is managed by the edge servers that aggregate model updates from a
group of clients and send their own model update to the principal server. Upon
receiving all model updates from the servers, the principal server executes
the second-level aggregation and generates the global model for the new round.
The global model is then distributed to all servers, and eventually to all
clients. Previous hierarchical FL frameworks are synchronous. In geo-
distributed settings, where the servers would be located far from each other,
a round in a hierarchical FL framework involves several successive long
distance communications and model aggregations, which would limit its speed.

##  3\. Overview of Spyker

In this section, we start the presentation of Spyker by providing the key
ideas behind its design.

Flat multi-server architecture. To get the best of both the asynchronous and
hierarchical FL paradigms, Spyker relies on multiple servers that directly
interact with each other and with clients in a fully asynchronous manner. The
servers are organized in a flat multi-server infrastructure (2(c)). Each
server is the center of a star sub-network, and interacts with a group of
clients that are assigned to it based on geographical proximity to minimize
communication latency. In this way, communication latency between clients and
servers are reduced, and servers are always able to process client or server
updates. On the contrary, the related studies mainly consider that clients
interact with a single server.

Asynchronous training and communications. Clients interact with their assigned
server in the same way they would in a classical asynchronous FL system: they
receive a global model from the server, train it over their local dataset and
send their model update back to the server. Spyker’s servers rely on
asynchronous communications to exchange their models and benefit from the
dataset of their respective clients. To minimize the complexity of this
procedure, servers rely on a token-based strategy to trigger the asynchronous
exchange of their models. Model aggregation also happens asynchronously and in
a peer-to-peer fashion among servers. The model exchange can only be initiated
by the server that holds the token, which prevents potential concurrent and
redundant model broadcasts. At a given point in time, servers might maintain
slightly different models.

Handling model staleness and fast clients. Since all interactions in Spyker
are asynchronous, servers take into account model staleness at all levels to
maintain accurate models. To do so, servers maintain the age of their model. A
server increases the age of its model whenever it processes a client model
update, or the model of another server. Spyker uses the age of models to
dampen the contribution of an old model update when processing it. Similarly,
Spyker uses a learning rate decay to limit the influence of clients that
frequently produce updates.

Preventing server models from drifting apart. The frequency at which servers
synchronize their models significantly influences the performance of the
system in presence of data and resource heterogeneity. It is therefore
important to ensure that a server model does not become biased toward the
datasets of its clients, and that it leverages all client datasets. The model
exchange and aggregation procedures allow servers to asynchronously homogenize
their models. To prevent server models from drifting too much from each other
during the training process, servers exchange and aggregate their models
whenever they detect that their ages differ too much, and whenever a server
has updated its model too many times since the last model exchange.

##  4\. Spyker Details

In this section, we present the details of Spyker’s building blocks. We first
detail how servers and clients interact. We then discuss Spyker’s aggregation
algorithm that merges two server models and the token-based algorithm that
triggers the asynchronous exchange of server models. Tab. 2 details our
notations and the values of all parameters.

###  4.1. Local Training and Decaying Learning Rates

Alg. 1 describes the asynchronous interactions between a client and a server
in Spyker, which are divided into two main procedures.

The first procedure, LocalTraining, is executed by a client $C_{k}$ to train a
model over its local dataset (Alg. 1, ll. 4-10), and is triggered upon
receiving a model $W_{i}^{t}$, along with its age $A_{i}$ and a learning rate
$\eta_{k}$ from a server. The client then trains the model over its dataset
using the specified learning rate (Alg. 1, ll. 7-8), and sends the trained
model back to the server along with the age $A_{i}$ (Alg. 1, l. 10). Note that
instead of having the client sends the model’s age with its update, the server
could remember the age of the model it has sent to every client. We adopted
this option to simplify our pseudocode.

The second procedure, Aggregation, is executed by a server upon receiving a
client update (Alg. 1, ll. 12-19). The server first determines the age
difference between its current model and the model it sent to the client (Alg.
1, l. 14), and uses it to possibly decrease the impact of the received update
on its model when updating it (Alg. 1, l. 15). The server then computes the
learning rate $\eta_{k}$ that the client should use for its next local
training (Alg. 1, l. 18), before returning its new model $W_{i}^{t+1}$ along
with $A_{i}$ and $\eta_{k}$ to client $C_{k}$ (Alg. 1, l. 19). Finally, the
server verifies whether it should trigger a synchronization of server models
(Alg. 1, l. 20) if its model age has sufficiently increased since the latest
synchronization by calling function $\textsc{checkSynchronization}()$ (defined
at l. 20 in Alg. 2).

Algorithm 1 Interactions between clients and servers

1:

2:

3:

4:[Client $C_{k}$] Procedure LocalTraining($W^{t}_{i},A_{i},\eta_{k}$)

5:
${\color[rgb]{0.69140625,0.1875,0.375}\definecolor[named]{pgfstrokecolor}{rgb}{0.69140625,0.1875,0.375}\triangleright}$receive
model $W^{t}_{i}$ with age $A_{i}$ from server $S_{i}$

6: $W_{k}^{t}=W^{t}_{i}$

7: for each epoch $\in T_{k}$ do

8: update $W_{k}$ with learning rate $\eta_{k}$

9: $W_{k}^{t+1}=W_{k}^{t}$

10: send $\left(W_{k}^{t+1},A_{i}\right)$ to server $S_{i}$

11:

12:[Server $S_{i}$] Procedure Aggregation($W_{k}^{t^{\prime}},A_{k}$)

13:
${\color[rgb]{0.69140625,0.1875,0.375}\definecolor[named]{pgfstrokecolor}{rgb}{0.69140625,0.1875,0.375}\triangleright}$receive
model $W_{k}^{t^{\prime}}$ with age $A_{k}{=}t^{\prime}$ from client $C_{k}$

14: $w_{k}^{t}=A_{i}-A_{k}$

15: $W^{t+1}_{i}=W^{t}_{i}+\eta_{i}\cdot
w^{t}_{k}\cdot\left(W_{k}^{t^{\prime}}-W_{i}^{t}\right)$

16: $A_{i}=A_{i}+1$

17: $u[k]=u[k]+1$
${\color[rgb]{0.69140625,0.1875,0.375}\definecolor[named]{pgfstrokecolor}{rgb}{0.69140625,0.1875,0.375}\triangleright}$num.
updates received per client

18: $\eta[k]=\mathrm{Decay}(\eta[u[k]],u[k],\overline{u})$

19: send $\left(W^{t+1}_{i},A_{i},\eta[k]\right)$ to client $C_{k}$

20: $\textsc{checkSynchronization}()$

Let us provide more information on the way a server in Spyker updates the
learning rate of a client (Alg. 1, l. 18) to maintain high model accuracy
despite clients and network heterogeneity, i.e., despite the fact that some
fast clients might consistently send updates at short intervals, while slow
clients might contribute infrequently. Spyker adopts an adaptive strategy that
tailors the impact of each client update. More precisely, a server uses
function $\mathrm{Decay}$ (Alg. 1, l. 18) to decrease the impact of the fast
clients that frequently generate updates on its model. To do so, the server
maintains the number of updates it has received from each client in an array
$u$, and therefore increments $u[k]$ (Alg. 1, l. 17) upon receiving a model
from $C_{k}$. The server adjusts client $C_{k}$’s learning rate using the
following $\mathrm{Decay}$ function:

$\mathrm{Decay}(\eta[u[k]],u[k],\overline{u})=$ |
---|---
$\left\\{\begin{aligned} &\eta_{k}\ &\mathrm{if}\ u[k]<\overline{u}\\\ &\max{\left(\eta_{min},\eta\left[u[k]\right]-\beta(u[k]-\overline{u})\right)}&\mathrm{if}\ u[k]\geq\overline{u}\end{aligned}\right.$ |

The $\mathrm{Decay}$ function takes three inputs: $u[k]$ is the number of
updates that the server received from client $C_{k}$; $\eta[u[k]]$ is the
learning rate a client would use without decay, which classically decreases
with $u[k]$; and $\overline{u}$ is the average number of updates that clients
have sent to server $S_{i}$. $\mathrm{Decay}$ also uses two parameters:
$\eta_{min}$ is a lower bound for learning rates; and $\beta$ is the decaying
rate. We empirically determined that $\beta=0.05$ and $\eta_{min}=0.000001$
provide the best results (cf. Tab. 2).

###  4.2. Token-Triggered Server Model Aggregations

Algorithm 2 Token-triggered aggregation of server models.

1:[Server $S_{i}$] ServerInit()

2: $token=(i{\neq}1)?\ NULL:\\{bid{=}1,ages{=}[0,\cdots,0]\\}$

3: $cnt=ages_{i}=\\{\\}$
${\color[rgb]{0.69140625,0.1875,0.375}\definecolor[named]{pgfstrokecolor}{rgb}{0.69140625,0.1875,0.375}\triangleright}$hashmaps

4: $didBroadcast=\emptyset$

5: $ongoingSynchro=\mathrm{False}$

6: $A_{i,prev}=A_{i}=0$

7: set $W_{i}^{0}$ to random model

8:

9:[Server $S_{i}$] Procedure RcvAge($A_{j}$)

10: ${ages}_{i}[j]=\mathrm{max}({ages}_{i}[j],A_{j})$

11: $\textsc{checkSynchronization}()$

12:

13:[Server $S_{i}$] Procedure RcvToken(t)

14: for $j\in[1,N]$ do

15: ${ages}_{i}[j]=\mathrm{max}({ages}_{i}[j],t.ages[j])$

16: $token=t$

17: $t.{bid}=t.{bid}+1$

18: $\textsc{checkSynchronization}()$

19:

20:[Server $S_{i}$] Procedure checkSynchronization()

21: if $\mathrm{max}(ages){-}\mathrm{min}(ages){\geq}h_{\mathrm{inter}}$ or
$A_{i}{-}A_{i,pre}{\geq}h_{\mathrm{intra}}$ then

22: if $hasToken$ and (not $ongoingSynchro$) then

23: $A_{i,pre}=A_{i}$

24: $ongoingSynchro=\mathrm{True}$

25: send $(W_{i}^{t},A_{i},t.bid)$ to all servers

26: $didBroadcast=didBroadcast\cup\\{t.bid\\}$

27: ${cnt}_{i}[token.bid]=1$

28: else

29: send $A_{i}$ to all servers

30:[Server $S_{j}$] Procedure RcvModel($W_{i}^{t},A_{i},{bid}_{i}$)

31: ${ages}_{j}[i]=\mathrm{max}({ages}_{j}[i],A_{i})$

32: if not ${bid}_{i}\in didBroadcast$ then

33: $didBroadcast=didBroadcast\cup\\{{bid}_{i}\\}$

34: $A_{j,pre}=A_{j}$

35: send $\left(W_{j}^{t},A_{j},{bid}_{i}\right)$ to all servers

36: $\textsc{ServerAgg}(W_{i}^{t},A_{i})$

37: if $token\neq\mathrm{NULL}$ and $t.bid={bid}_{i}$ then

38: $cnt[{bid}_{i}]=cnt[{bid}_{i}]+1$

39: if $cnt[{bid}_{i}]=n$ then

40: $t.ages=ages_{j}$

41: send token $t$ to the next server on ring

42: $token=\textrm{NULL}$

43: $ongoingSynchro=\textrm{False}$

44:

45:[Server $S_{i}$] Procedure ServerAgg($W_{j}^{t},A_{j}$)

46: process model $W_{j}^{t}$ of age $A_{j}=t$ from server $S_{j}$

47: $a=\frac{\phi(A_{j}-A_{i})}{A_{i}}$

48: $w_{i,j}={\left(1+e^{-a}\right)}^{-1}$

49: $W_{i}^{t+1}={W_{i}^{t}+\eta_{a}\cdot
w_{i,j}\left(W_{j}^{t}-W_{i}^{t}\right)}$

50: $A_{i}=\left(1-\eta_{a}\cdot w_{i,j}\right)A_{i}+\eta_{a}\cdot
w_{i,j}\cdot A_{j}$

In addition to interacting with their clients, Spyker’s servers synchronize
their models using Alg. 2. Since the age of a server model impacts the way its
client updates are aggregated, Spyker triggers an asynchronous server model
synchronization whenever the ages of server models are drifting too much or
when a server has processes a large enough number of client updates to make
sure that all client datasets are fairly represented in server models.

To avoid the complexity of having to handle multiple concurrent
synchronizations, servers rely on a token so that only one server can trigger
the model broadcasts and aggregations. This token contains a synchronization
ID $bid$ that allows server to broadcast their model only once per
synchronization. Only the server holding the token can trigger a
synchronization, but other servers might indirectly learn about an ongoing
synchronization and broadcast their models to each other. A server broadcasts
its model either because it holds the token and because some conditions are
met (more details below), or by receiving another server’s model with an
unknown synchronization ID. To maintain fairness, the token is circulated
among all servers. To accelerate synchronizations, whenever necessary, servers
can broadcast their age so that the token holder triggers a synchronization.
Note that we assume that links are FIFO, which can easily be enforced if it is
not the case, so that a server receives the models of any other server
according to increasing synchronization IDs and can process them all.

Upon executing the ServerInit initialization procedure (Alg. 2, l. 1), a
server $S_{i}$ initializes several variables to coordinate the exchange of
models between servers. The token initially resides at a randomly chosen
server (here server $S_{1}$), and contains a synchronization ID $bid$ set to
$1$ and a vector of server model ages that are all initially equal to $0$. The
hashmap $cnt$ is used to count the number of models that have been received
for a given synchronization ID. The known age of server models is kept in
${ages}_{i}$. Variable $didBroadcast$ contains the set of synchronization IDs
for which $S_{i}$ has broadcast its model. The Boolean variable
$ongoingSynchro$ is used by the server holding the token to initiate a
synchronization only once before forwarding the token. The model age at which
the latest synchronization happened, and the current model age are
respectively stored in $A_{i,prev}$ and $A_{i}$. Finally, a server randomly
initializes its model.

The token keeps visiting every server following a random ring-based topology.
Upon arrival of the token (Alg. 2, ll. 13– 18), a server $S_{i}$ updates the
ages of the models that it maintains using the token’s information. It then
calls the CheckSynchronization procedure (Alg. 2, l. 20) to determine whether
a model synchronization should be triggered. The frequency at which server
model exchanges happen is controlled by the $h_{\mathrm{inter}}$ and
$h_{\mathrm{intra}}$ parameters (cf. Tab. 2): $h_{\mathrm{inter}}$ is the
threshold for the maximum age difference between different server models; and
$h_{\mathrm{intra}}$ is the threshold for the maximum age difference between a
server model’s current age and the age it had during the last server model
synchronization. If the inter-server age difference exceeds threshold
$h_{\mathrm{inter}}$, or if the intra-server age difference exceeds threshold
$h_{\mathrm{intra}}$, then server $S_{i}$ might trigger a server model
exchange. First, if $S_{i}$ holds the token and has not already triggered a
server model exchange, then it broadcasts its current model $W_{i}^{t}$ of age
$A_{i}$ along with the current broadcast ID $t.bid$ to all the other servers.
Second, if $S_{i}$ does not hold the token, it can broadcast its model age to
all other servers, so that the server that holds the token can execute
procedure RcvAge (Alg. 2, l. 9), which updates the age of a server model and
possibly triggers a model synchronization.

Procedure RcvModel (Alg. 2, ll. 30-41) is executed when a server receives a
server model, and is also in charge of relaying the token. Upon receiving a
server model from a server $S_{i}$, server $S_{j}$ first updates its local age
for server model $i$. It then verifies whether it has already broadcast is
model for the synchronization ID indicated by server $S_{i}$ using its local
$didBroadcast$ set, and if necessary broadcasts its model and updates the
relevant variables. A server that receives a server model then executes the
aggregation procedure ServerAgg (Alg. 2, l. 45), which is described in details
in the next section. Finally, the server that holds the token makes sure that
it receives the model of all other servers (using the $cnt$ hashmap) before
forwarding the token to its successor on the ring of servers.

###  4.3. Aggregation of Server Models

Following their exchange of models, servers aggregate them asynchronously
immediately after they are received. Alg. 2 also details the procedures that a
server $S_{i}$ follows to aggregate a model that it receives from another
server $S_{j}$. This algorithm first calculates the weight that should be used
to update the local server model with the received one. The aggregation weight
$w_{i,j}$ is computed based on the age of the two server models $A_{i}$ and
$A_{j}$ using a sigmoid function that is parameterized using parameter $\phi$
(Alg. 2, ll. 47-48):

$w_{i,j}=\frac{1}{1+e^{-a}},\quad\text{where}\quad
a=\frac{\phi(A_{j}-A_{i})}{A_{i}}$

The age difference $A_{j}-A_{i}$ indicates whether $W_{j}^{t}$ is more mature
than $W_{i}^{t}$ in terms of the number of updates both models have been
trained on. The influence of $W_{j}^{t}$ increases with this age difference.
Denominator $A_{i}$ makes the absolute weight relative to the current age of
server $S_{i}$. As a model age increases, it becomes more stable and the
impact of other server models should be decreased. The sigmoid function
ensures that $w_{i,j}$ remains between $0$ and $1$. The derivative of this
function becomes $0$ and results in a weight of $1$ when the relative model
age difference is too large. Parameter $\phi$ indicates in which range the
sigmoid function is active. A larger $\phi$ leads to a smaller active range,
leaving a larger area for the weight of 0 and 1.

The aggregation rate $\eta_{a}$ also scales the influence of other servers
during model aggregation. When the aggregation with server $S_{i}$ happens
locally at server $S_{j}$, a large $\eta_{a}$ indicates that server $S_{i}$
could greatly influence server $S_{j}$’s model, and $\eta_{a}=1$ means it
could replace server $S_{j}$’s model when their ages are equal. Parameter
$\eta_{a}$ needs to be carefully tuned. If $\eta_{a}$ is too large, the
influence of the server itself is eliminated; if $\eta_{a}$ is too small, the
server learns too little from its peers thus its model could be biased to its
clients’ data distribution. We empirically determined $\phi=1.5$ and
$\eta_{\alpha}=0.6$ to provide the best results (cf. Tab. 2).

Server $S_{i}$ finally updates its age after aggregating another server model.
According to the weight $w_{i,j}$ that has been applied to the model
aggregation, server $S_{i}$ updates its local model age (Alg. 2, l. 50). As
the aggregation of a server model embeds more than one client update, we apply
a weighting strategy to increase a model’s age after aggregating another
server’s model.

##  5\. Performance Evaluation

Figure 3. WikiText2: Perplexity wrt. time. (lower is better)

Figure 4. WikiText2: Perplexity wrt. # updates (lower is better)

In this section, we describe our experimental settings, and compare Spyker to
state-of-the-art representative FL algorithms. We show the superior accuracy
of Spyker over previous baselines, demonstrate its improved scalability with
the number of servers, and report the positive impact of our learning rate
decay function on accuracy. We also report the network consumption of all
algorithms.

###  5.1. Settings

Baselines. We compare Spyker with three recent Federated Learning algorithms:
FedAvg (McMahan et al., 2017), FedAsync (Xie et al., 2019) and HierFedAvg (Liu
et al., 2020). FedAvg is the original synchronous FL framework, while FedAsync
and HierFedAvg are respectively state-of-the-art asynchronous and hierarchical
FL frameworks. We also implemented a synchronous version of Spyker, which we
call Sync-Spyker, that uses a synchronous server-server model aggregation
procedure. Sync-Spyker periodically triggers synchronous model exchanges
between servers (i.e., ignoring the influence of $h_{inter}$ to trigger model
exchanges). Periodically, the servers broadcast their model, along with their
age, and wait for all other server models. Upon receiving all models, a server
aggregates them following a deterministic order (i.e., using the server IDs).
Following this exchange and aggregation of server models, all servers own the
same model. In Sync-Spyker, when the server synchronization algorithm is
started, servers stop processing local updates from clients, and instead store
them. Client updates are then processed when all the server models have
arrived and have been aggregated.

Datasets. We conduct experiments on the MNIST and CIFAR-10 image collections,
and on the WikiText2 language modeling dataset. For a given experiment, a
dataset is split into subsets of equal sizes that are assigned to different
clients. To introduce client data heterogeneity, we assign $l$ labels to each
client. A smaller $l$ indicates a higher degree of non-independent and
identically distributed (IID) data distribution and thus larger heterogeneity
among clients. We set $l$ to 2 in non-IID data experiments.

Computation and Network Delays. Our experiments take into account the
computations that clients execute. Time is maintained at each client, and is
advanced depending on the procedure it executes during the experiments. Tab. 3
details the average computing times of the server aggregation procedure of all
FL frameworks. These values were obtained by benchmarking all algorithms using
the Python $time$ package. To simulate client heterogeneity, we sample each
client’s training delay from a Gaussian distribution $N(\mu,\sigma^{2})$,
where $\mu=150$ and $\sigma=7.5$. We use the same training delays per client
across all experiments. Since servers might be located in different regions of
the world (e.g., Amsterdam and Sydney), we set the communication delays among
servers according to the Amazon Web Service network latency (Services, 2020)
as shown in Tab. 4. We assume that a client and its nearest server are in the
same geographical location and therefore set their communication delay using
the diagonal of Tab. 4. We assume that all network links have a 100 Mbps
bandwidth. For each experiment, we detail the number of servers and clients
that we use along with the parameter values.

Table 3. Computation time required per procedure (ms). Local Training | 200
---|---
Model Aggregation in Sync-Spyker | 2
Model Aggregation in Spyker | 2
Model Aggregation in FedAVg | 15
Model Aggregation in HierFAVG | 15
Model Aggregation in FedAsync | 2
Table 4. Communication delays between geographical locations (ms). | Hongkong | Paris | Sydney | California
---|---|---|---|---
Hongkong | 1.41 | 194.9 | 132.28 | 155.13
Paris | 197.91 | 0.9 | 278.83 | 142.25
Sydney | 132.06 | 280.11 | 2.56 | 138.47
California | 154.96 | 142.79 | 138.57 | 2.14

Figure 5. MNIST: Accuracy wrt. time (higher is better)

Figure 6. MNIST: Accuracy wrt. # updates (higher is better)

Figure 7. CIFAR-10: Accuracy wrt. time (higher is better)

Figure 8. CIFAR-10: Accuracy wrt. # updates (higher is better)

Model. For MNIST, we use a CNN model with two convolutional and two fully
connected layers. For CIFAR-10, we use a CNN model with three convolutional
and two fully connected layers. For WikiText2, we use the next character Long
Short-Term Memory (LSTM) neural network model, which has been designed for
character-level text generation tasks. The language model consists of an
embedding layer that is used to capture the semantic and syntactic properties
of characters, an LSTM layer that captures dependencies between characters and
generates coherent text, and a fully connected layer that transforms the
LSTM’s hidden states into a probability distribution over all possible
characters in the vocabulary. The initial local learning rate of clients
$\eta_{k}$ is $0.05$. In asynchronous settings, we use $\alpha=0.5$ for the
staleness weighting in FedAsync and a global learning rate of $\eta=0.6$ for
the client-server update.

###  5.2. Accuracy

For the following experiments, we use $100$ clients that are equally
distributed across $4$ servers located in the AWS region we consider. Fig. 3,
Fig. 5 and Fig. 7 show that Spyker converges faster than other baselines in
terms of elapsed time for all three datasets. Spyker and Sync-Spyker converge
as quickly as other baselines in terms of number of processed client updates.

Table 5. Multiplicative factor for the amount of time and number of updates necessary for algorithms to reach a $90\%$ accuracy with $200$ and $300$ clients compared with their own execution with $100$ clients. A low value indicates that an algorithm scales well with the number of clients.  | 200 clients | 300 clients
---|---|---
| Time | Updates | Time | Updates
FedAsync | 2.42 | 3.63 | 5.99 | 9.01
FedAvg | 1.98 | 4.05 | 2.44 | 7.26
HierFAVG | 1.64 | 3.18 | 1.75 | 8.19
Spyker | $\mathbf{1.21}$ | $\mathbf{2.43}$ | $\mathbf{1.43}$ | $\mathbf{4.34}$
Sync-Spyker | 1.61 | 2.68 | 1.90 | 4.87

We also evaluated the ability of the $5$ FL frameworks to scale with the
number of clients. We evaluated the time each algorithm needs to reach $90\%$
accuracy with MNIST and $4$ servers using $100$, $200$ and $300$ clients. Tab.
5 reports this time for all algorithms for $200$ and $300$ clients divided by
the time measured with $100$ clients. When the number of clients in the system
increases, the convergence time and required number of updates of the various
methods increase following different trends. Spyker is the method that scales
the best. For example, it requires $21\%$ more time to converge with $200$
clients than with $100$ clients while other baselines (excluding Sync-Spyker)
require at least $64\%$ more time. It also requires $2.43$ times as many
updates to converge while other baselines require at least $3.18$ times as
many. On the other side of the spectrum, FedAsync’s required time and number
of updates to converge increase the most. With $200$ clients, its convergence
time is multiplied by $2.42$, while its required number of updates is
multiplied by $3.63$

Spyker not only exhibits the fastest convergence, but it also scales the best
with the number of clients.

Fig. 4, Fig. 6 and Fig. 8 show the perplexity of all the baselines, Spyker and
Sync-Spyker on WikiText-2, and their accuracy on MNIST and CIFAR-10 depending
on the number of client updates processed. One can see that Spyker and Sync-
Spyker do not always require the least number of updates to reach a given
perplexity or accuracy. This is not a problem in practice as it makes more
sense to evaluate the evolution of perplexity or accuracy based on elapsed
time during a deployment. However, these figures also show that Spyker reaches
similar perplexity or accuracy levels as FedAvg, which generally converges the
fastest based on the number of updates (but not based on elapsed time, where
it is the worst performing baseline).

Figure 9. Number of updates queued at a server with Spyker and FedAsync

Figure 10. Density distribution of the number of updates per client with
Spyker and FedAsync

###  5.3. Impact of Multiple Servers in Asynchronous FL

In asynchronous FL systems, the server aggregates an update immediately after
receiving it. However, a busy server can cause a queueing of updates that
await being processed. To further illustrate the benefits of introducing a
multi-server system to an asynchronous FL framework, we evaluate the update
queuing phenomenon with Spyker ($4$ servers) and FedAsync ($1$ server) with
$200$ clients. The local training delays of clients are generated using a
Gaussian distribution with a mean of $150$ ms and a standard deviation of $60$
ms. Processing delays are set using Tab. 4.

Fig. 9 shows that FedAsync’s queue length increases to nearly $80$ updates
from $300$ ms to $600$ ms and always stays above $20$ during the rest of the
experiment. Client updates are processed more efficiently by the $4$ servers
that Spyker uses, as server queues never contain more than $20$ updates.

We further conducted a group of experiments to evaluate the time needed for
FedAsync and Spyker to reach $90\%$ and $95\%$ accuracy on the MNIST dataset.
For these experiments, the time difference observed is only caused by resource
heterogeneity as we set all network latency’s to the same value. The results
are shown in the bottom three rows of Tab. 6. Spyker reaches a $90\%$ and a
$95\%$ accuracy respectively $38\%$ and $25\%$ faster than FedAsync.

Table 6.  Time to reach $90\%$ or $95\%$ accuracy with FedAsync and Spyker. Lat uses AWS latency (Services, 2020), and No lat assumes uniform latency between servers (equal average in both cases).  Network | Method | Time 90% | Time 95%
---|---|---|---
Lat. | FedAsync | $59s$ | $125s$
Spyker | $22s$ | $51s$
Improvement | -61% | -58%
No lat. | FedAsync | $40s$ | $75s$
Spyker | $25s$ | $56s$
Improvement | -38% | -25%

To demonstrate that Spyker also mitigates the impact of high communication
delays, we conducted a similar group of experiments that, in addition,
introduce network latencies (cf. Tab. 4). The results are shown in the top
three rows of Tab. 6. Compared with the bottom three rows, the top three rows
show a larger difference between Spyker and FedAsync. Spyker appears $61\%$
faster to reach $90\%$ accuracy, and $58\%$ faster to reach $95\%$ accuracy
than FedAsync. This difference can be explained by the fact that Spyker
decreases the distance between servers and clients compared to FedAsync, which
uses a single server.

Figure 11. Accuracy with and without a learning rate decay on MNIST with $4$
servers and $100$ clients ($25$ per server)

Figure 12. Overall network consumption of Spyker compared to the baselines

Table 7. Effect of imbalanced number of clients per server on accuracy.  | Spyker
---|---
| 25 clients | 52 clients | 63 clients | 70 clients
Accuracy | $96.4\%$ | $-14.5\%$ | $-11.9$ | $-7.63\%$
Duration (s) | 120 | $+17.76$ | $+33.04$ | $+53.08$

###  5.4. Number of Clients per Server

We evaluate the effect of imbalanced number of clients allocated to each
server. We consider four different scenarios. The first one consists of 4
servers with each 25 clients. In the three other scenarios, one of the servers
has more clients (i.e., it has 52, 63, or 70 clients) while the remaining
clients are divided evenly over the remaining servers. Table 7 show the
accuracy of the global model and the time necessary for it to converge. These
results show that slight imbalance in the client distribution degrades
accuracy the most (-14.5% with the second most balanced scenario) and that
increasing imbalance always increases the convergence time. In comparison, we
observed that HierFAVG’s performance is less impacted by this kind of
imbalance, but its performance nonetheless is always worse than Spyker’s.

###  5.5. Impact of Learning Rate Decay

To demonstrate the positive impact of our learning rate decay method on
Spyker, we conducted two experiments. We first used Spyker and FedAvg each
with $100$ clients to perform the classification task on the MNIST dataset.
For these two algorithms that represent the single-server and multi-server
baselines, we gathered statistics on the number of updates of every client and
plotted a Kernel Density Estimate (KDE) plot to show its distribution. The
most desirable KDE plot for a FL system should exhibit a single concentrated
peak, which would mean that clients contribute equally. Fig. 10 shows the KDE
plots of Spyker and FedAsync.

FedAsync has a more centralized distribution of updates with a steep peak at
around $200$ updates compared to a gentler peak for Spyker at around $300$
updates. The reason for this difference lies in the fact that the major source
of delay in FedAsync is communication, and the range of possible communication
delays with clients is wider with a single server. On the other hand, by
introducing multiple servers, the range of possible communication delays with
Spyker becomes smaller and the differences between client training times
becomes the major source of heterogeneity. We also observe a peak at around
$1,400$ updates for Spyker and $1,700$ for FedAsync, which comes from the
presence of fast clients near the server, which could cause a biased global
model due to their higher influence on the server’s model if left unaddressed.
This analysis justifies the need for a client learning rate decay mechanism
that customizes the learning rate of clients according to the number of
updates they contribute. The impact of the updates that the most active
clients generate is therefore dampened to balance the overall contribution of
clients. Interestingly however, the decay function would use a smaller decay
ratio if the system as a whole would catch up with the fast clients.

In the second experiment, we evaluated Spyker’s convergence speed with and
without our learning rate decay mechanism. Fig. 11 shows that introducing the
learning rate decay mechanism increases Spyker’s convergence speed, and
therefore that it reduces the impact of client data heterogeneity.

###  5.6. Bandwidth consumption

We evaluated the bandwidth consumption of all FL algorithms over time. We
measured the number of bytes transferred over $110$ s during server-server and
server-client model communications. Fig. 12 details the network consumption of
all algorithms with the MNIST dataset, $4$ servers and $100$ clients equally
distributed over all servers. FedAvg has the lowest bandwidth consumption with
$2.28$ MB consumed, which can be explained by its synchronous communications
and its use of a single server. Spyker has the highest bandwidth requirements
with $63.4$ MB transferred because of its asynchronous server-server and
server-client communications. This network consumption is however very
reasonable for modern networks. For each client update, the aggregated global
model is relayed back to the client, which results in higher bandwidth
requirements. Similarly, Sync-Spyker is the second most network intensive
algorithm and transmits $32.5$ MB over the network since it maintains
asynchronous client-server interactions. Unlike Spyker’s asynchronous multi-
server aggregation, Sync-Spyker’s synchronous multi-server aggregation halts
the client updates, which results in lower communication overhead. FedAsync’s
asynchronous client updates incur $24.5$ MB of data transmissions over the
network, which is higher than with FedAVG. HierFAVG’s synchronous client-
server communication results in lower communication overhead of $16.35$ MB but
the hierarchical server aggregation requires relatively higher bandwidth.

##  6\. Related Work

This section discusses additional related work on asynchronous and multi-
server FL that complements Sec. 2.

###  6.1. Asynchronous FL

ASO-Fed (Chen et al., 2020) is an asynchronous FL framework that enables wait-
free computation and communication. ASO-Fed allows client updates to stream
into the federator in different rounds and each local training is based on the
newly streamed-in data. FedAsync (Xie et al., 2019) is an asynchronous
optimization algorithm for FL that guarantees a near-linear convergence and
address the staleness problem in the asynchronous setting. Each time an update
arrives, the weighted averaging is applied for the current global model and
the update. AsyncFedED is an adaptive asynchronous FL aggregation based on
Euclidean distance (Wang et al., 2022). If the global model is updated before
a client sends back its update, the learning rate of the server should be
adjusted based on the staleness of the current update. AsyncFedED evaluates
the staleness of an update by its Euclidean distance from the server. Larger
staleness leads to a smaller learning rate for its aggregation. The staleness
of weight is the relative weight differences. FedAsync (Xie et al., 2019)
introduces a staleness strategy where the influence of a client’s update
decreases when its staleness grows as the inverse function of $\alpha$, the
adaptive parameter. A larger $\alpha$ indicates the influence of stale updates
should be reduced. FedAsync performs best when the staleness is small, and it
can converge faster than FedAvg.

Overall, asynchronous FL approaches incur a high computational load on the
server and therefore do not scale well with the number of clients. Intelligent
node selection techniques (Zhou et al., 2021; Hao et al., 2020; Chen et al.,
2021) aim at alleviating this burden but are often biased towards specific
clients.

###  6.2. Multi-Server FL

Existing multi-server methods (Azimi-Abarghouyi and Fodor, 2022; Li et al.,
2020b; Wang and Wang, 2022; Abad et al., 2020; Liu et al., 2020) add edge
servers or assign proxy clients to undertake part of the workload.
Hierarchical FL approaches have also been described (Abad et al., 2020; Liu et
al., 2020; Wang and Wang, 2022; Li et al., 2019). Qu et al. demonstrated the
convergence of a multi-server FL algorithm (Qu et al., 2023) where servers
indirectly synchronize their models through overlapping client sets.
Clustering is frequently used in multi-task and unsupervised FL (Adolfsson et
al., 2019; Smith et al., 2017; Dennis et al., 2021). Incorporating clustering
algorithms in our framework for potential refinement is future work.

Xie et al. proposed a multi-center FL algorithm, FeSEM (Xie et al., 2020), to
address the challenge of data heterogeneity. The main idea is to cluster the
updates and assign them to their closest global model (i.e. center). They
showed in experiments that the clusters of local updates characterize clients’
data distribution thus models generated from similar data distributions are
gathered on the same center and aggregated.

HierFAVG (Liu et al., 2020) is a client-edge-cloud hierarchical framework,
which can be considered as a multi-server version of the FedAvg algorithm.
Each edge server employs averaging aggregation for its clients. After every
certain number of rounds, the cloud server applies averaging strategy to the
edge servers. The introduction of the hierarchical structure reduces the
communication burden on servers and improves the capability of the single-
server FedAvg algorithm while maintaining the stability of synchronous
systems.

Despite their ability to support a larger number of clients, current multi-
server and hierarchical-server algorithms do not directly address the
sensitivity to large system heterogeneity of synchronous systems.

###  6.3. Geo-distributed Machine Learning

Many distributed ML systems target large-scale ML applications (Dean et al.,
2012; Chilimbi et al., 2014; Cui et al., 2016), however, they assume that
network communication happens within a datacenter. Cano et al. (Cano et al.,
2016) discussed running a machine learning system in geo-distributed settings.
They show that one can leverage a communication-efficient algorithm for
logistic regression models (Mahajan et al., 2013) to improve performance.
Several works (Jaggi et al., 2014; Zhang and Lin, 2015; Hsieh et al., 2017)
focused on the design of communication-efficient mechanisms that leaves the ML
algorithms unmodified. These works are orthogonal to our approach, since we
design asynchronous multi-server FL algorithms for the geo-distributed
settings.

##  7\. Conclusion

We described Spyker, the first fully asynchronous multi-server FL algorithm.
Spyker focuses on practical scenarios where clients are distributed over the
world, and it addresses the performance bottlenecks that naturally appear in
presence of weak clients or heterogeneous networks. Spyker scales better with
the number of clients than previous works. Spyker relies on an aging
mechanisms that servers use to support asynchrony of all exchanges and on a
token-based algorithm that servers rely on to synchronize their model
exchanges. Our experimental results show that Spyker requires less time to
converge than three representative algorithms – FedAsync, HierFAVG, and
FedAvg, and is also faster than its partially synchronous variant Sync-Spyker
where servers synchronize themselves using a synchronous algorithm. Spyker
also leverages a learning rate decay method to maintain high accuracy despite
client heterogeneity. Future work includes exploring the possibility of
integrating clustering algorithms in Spyker to enable servers to group clients
based on possible similarities in their data distributions.

## References

  * (1)
  * Abad et al. (2020) Mehdi Salehi Heydar Abad, Emre Ozfatura, Deniz Gündüz, and Özgür Erçetin. 2020\.  Hierarchical Federated Learning ACROSS Heterogeneous Cellular Networks. In _ICASSP_.  https://doi.org/10.1109/ICASSP40776.2020.9054634
  * Adolfsson et al. (2019) Andreas Adolfsson, Margareta Ackerman, and Naomi C. Brownstein. 2019.  To cluster, or not to cluster: An analysis of clusterability methods.  _Pattern Recognit._ 88 (2019), 13–26.  https://doi.org/10.1016/j.patcog.2018.10.026
  * Azimi-Abarghouyi and Fodor (2022) Seyed Mohammad Azimi-Abarghouyi and Viktoria Fodor. 2022.  Multi-Server Over-the-Air Federated Learning.  https://doi.org/10.48550/arXiv.2211.16162 arXiv:2211.16162
  * Brisimi et al. (2018) Theodora S. Brisimi, Ruidi Chen, Theofanie Mela, Alex Olshevsky, Ioannis Ch. Paschalidis, and Wei Shi. 2018\.  Federated learning of predictive models from federated Electronic Health Records.  _Int. J. Medical Informatics_ 112 (2018), 59–67.  https://doi.org/10.1016/j.ijmedinf.2018.01.007
  * Cano et al. (2016) Ignacio Cano, Markus Weimer, Dhruv Mahajan, Carlo Curino, and Giovanni Matteo Fumarola. 2016\.  Towards geo-distributed machine learning.  _arXiv:1603.09035_ (2016).
  * Chen et al. (2020) Yujing Chen, Yue Ning, Martin Slawski, and Huzefa Rangwala. 2020\.  Asynchronous Online Federated Learning for Edge Devices with Non-IID Data. In _BigData_.  https://doi.org/10.1109/BigData50022.2020.9378161
  * Chen et al. (2021) Zheyi Chen, Weixian Liao, Kun Hua, Chao Lu, and Wei Yu. 2021.  Towards asynchronous federated learning for heterogeneous edge-powered internet of things.  _Digit. Commun. Networks_ 7, 3 (2021), 317–326.  https://doi.org/10.1016/j.dcan.2021.04.001
  * Chilimbi et al. (2014) Trishul Chilimbi, Yutaka Suzue, Johnson Apacible, and Karthik Kalyanaraman. 2014.  Project adam: Building an efficient and scalable deep learning training system. In _11th USENIX symposium on operating systems design and implementation (OSDI 14)_. 571–582.
  * Cui et al. (2016) Henggang Cui, Hao Zhang, Gregory R Ganger, Phillip B Gibbons, and Eric P Xing. 2016.  Geeps: Scalable deep learning on distributed gpus with a gpu-specialized parameter server. In _EuroSys_.
  * Dean et al. (2012) Jeffrey Dean, Greg Corrado, Rajat Monga, Kai Chen, Matthieu Devin, Mark Mao, Marc’aurelio Ranzato, Andrew Senior, Paul Tucker, Ke Yang, et al. 2012\.  Large scale distributed deep networks.  _NeurIPS_ (2012).
  * Dennis et al. (2021) Don Kurian Dennis, Tian Li, and Virginia Smith. 2021\.  Heterogeneity for the win: One-shot federated clustering. In _ICML_.
  * Hao et al. (2020) Jiangshan Hao, Yanchao Zhao, and Jiale Zhang. 2020.  Time Efficient Federated Learning with Semi-asynchronous Communication. In _ICPADS_.  https://doi.org/10.1109/ICPADS51040.2020.00030
  * Hard et al. (2018) Andrew Hard, Kanishka Rao, Rajiv Mathews, Françoise Beaufays, Sean Augenstein, Hubert Eichner, Chloé Kiddon, and Daniel Ramage. 2018.  Federated Learning for Mobile Keyboard Prediction.  (2018).  arXiv:1811.03604
  * Hsieh et al. (2017) Kevin Hsieh, Aaron Harlap, Nandita Vijaykumar, Dimitris Konomis, Gregory R Ganger, Phillip B Gibbons, and Onur Mutlu. 2017\.  Gaia: Geo-Distributed machine learning approaching LAN speeds. In _NSDI_.
  * Huang et al. (2020) Zhiqi Huang, Fenglin Liu, and Yuexian Zou. 2020.  Federated Learning for Spoken Language Understanding. In _COLING_.  https://doi.org/10.18653/v1/2020.coling-main.310
  * Jaggi et al. (2014) Martin Jaggi, Virginia Smith, Martin Takác, Jonathan Terhorst, Sanjay Krishnan, Thomas Hofmann, and Michael I Jordan. 2014\.  Communication-efficient distributed dual coordinate ascent.  _NeurIPS_ (2014).
  * Li et al. (2019) Li Li, Haoyi Xiong, Zhishan Guo, Jun Wang, and Cheng-Zhong Xu. 2019.  SmartPC: Hierarchical Pace Control in Real-Time Federated Learning System. In _RTSS_.  https://doi.org/10.1109/RTSS46320.2019.00043
  * Li et al. (2020b) Tian Li, Anit Kumar Sahu, Manzil Zaheer, Maziar Sanjabi, Ameet Talwalkar, and Virginia Smith. 2020b.  Federated Optimization in Heterogeneous Networks.  (2020).  https://proceedings.mlsys.org/book/316.pdf
  * Li et al. (2020a) Xiaoxiao Li, Yufeng Gu, Nicha C. Dvornek, Lawrence H. Staib, Pamela Ventola, and James S. Duncan. 2020a.  Multi-site fMRI analysis using privacy-preserving federated learning and domain adaptation: ABIDE results.  _Medical Image Anal._ 65 (2020), 101765.  https://doi.org/10.1016/j.media.2020.101765
  * Liu et al. (2020) Lumin Liu, Jun Zhang, Shenghui Song, and Khaled B. Letaief. 2020\.  Client-Edge-Cloud Hierarchical Federated Learning. In _ICC_.  https://doi.org/10.1109/ICC40277.2020.9148862
  * Mahajan et al. (2013) Dhruv Mahajan, S Sathiya Keerthi, S Sundararajan, and Léon Bottou. 2013.  A functional approximation based distributed learning algorithm.  _arXiv:1310.8418_ (2013).
  * McMahan et al. (2017) Brendan McMahan, Eider Moore, Daniel Ramage, Seth Hampson, and Blaise Agüera y Arcas. 2017\.  Communication-Efficient Learning of Deep Networks from Decentralized Data. In _AISTATS_.
  * Qu et al. (2023) Zhe Qu, Xingyu Li, Jie Xu, Bo Tang, Zhuo Lu, and Yao Liu. 2023.  On the Convergence of Multi-Server Federated Learning With Overlapping Area.  _IEEE Transactions on Mobile Computing_ 22, 11 (2023), 6647–6662.  https://doi.org/10.1109/TMC.2022.3200016
  * Rieke et al. (2020) Nicola Rieke, Jonny Hancox, Wenqi Li, Fausto Milletarì, Holger R. Roth, Shadi Albarqouni, Spyridon Bakas, Mathieu N. Galtier, Bennett A. Landman, Klaus H. Maier-Hein, Sébastien Ourselin, Micah J. Sheller, Ronald M. Summers, Andrew Trask, Daguang Xu, Maximilian Baust, and M. Jorge Cardoso. 2020.  The future of digital health with federated learning.  _npj Digit. Medicine_ 3 (2020).  https://doi.org/10.1038/s41746-020-00323-1
  * Services (2020) Amazon Web Services. 2020\.  AWS Latency Monitoring.  https://www.cloudping.co/grid
  * Smith et al. (2017) Virginia Smith, Chao-Kai Chiang, Maziar Sanjabi, and Ameet Talwalkar. 2017.  Federated Multi-Task Learning. In _NeurIPS_.  https://proceedings.neurips.cc/paper/2017/hash/6211080fa89981f66b1a0c9d55c61d0f-Abstract.html
  * Wang et al. (2022) Qiyuan Wang, Qianqian Yang, Shibo He, Zhiguo Shi, and Jiming Chen. 2022.  AsyncFedED: Asynchronous Federated Learning with Euclidean Distance based Adaptive Weight Aggregation.  https://doi.org/10.48550/arXiv.2205.13797 arXiv:2205.13797
  * Wang and Wang (2022) Xing Wang and Yijun Wang. 2022.  Asynchronous Hierarchical Federated Learning.  https://doi.org/10.48550/arXiv.2206.00054 arXiv:2206.00054
  * Xie et al. (2019) Cong Xie, Sanmi Koyejo, and Indranil Gupta. 2019.  Asynchronous Federated Optimization.  _CoRR_ abs/1903.03934 (2019).  arXiv:1903.03934 http://arxiv.org/abs/1903.03934
  * Xie et al. (2020) Ming Xie, Guodong Long, Tao Shen, Tianyi Zhou, Xianzhi Wang, and Jing Jiang. 2020\.  Multi-Center Federated Learning.  (2020).  arXiv:2005.01026 https://arxiv.org/abs/2005.01026
  * Xu et al. (2021) Jie Xu, Benjamin S. Glicksberg, Chang Su, Peter B. Walker, Jiang Bian, and Fei Wang. 2021\.  Federated Learning for Healthcare Informatics.  _J. Heal. Informatics Res._ 5, 1 (2021), 1–19.  https://doi.org/10.1007/s41666-020-00082-4
  * Zhang and Lin (2015) Yuchen Zhang and Xiao Lin. 2015.  DiSCO: Distributed optimization for self-concordant empirical loss. In _ICML_.
  * Zhou et al. (2021) Chendi Zhou, Hao Tian, Hong Zhang, Jin Zhang, Mianxiong Dong, and Juncheng Jia. 2021\.  TEA-fed: time-efficient asynchronous federated learning for edge computing. In _Computing Frontiers Conference_.  https://doi.org/10.1145/3457388.3458655
  * Zhu et al. (2022) Jie Zhu, Shenggui Li, and Yang You. 2022.  Sky Computing: Accelerating Geo-distributed Computing in Federated Learning.  _arXiv:2202.11836_ (2022).


# An Optical Analysis of Sunspots as Predictors of Geomagnetic Storms

Matthew Shelby, Scott Scharlach, Petar Matejic, RJ Everett, Colton Morgan
Pisgah Astronomical Research Institute
1 PARI Drive
Rosman, NC 28772, USA

###### Abstract

Although a variety of phenomena may create a geomagnetic storm on Earth, the
most severe geomagnetic storms arise from solar activity, and in particular,
coronal mass ejections (CMEs) and solar flares. CMEs and flares originate
primarily from sunspots. The “aa index” is a metric which ranks all of the
strongest geomagnetic storms between 1868 and 2010 based on a variety of
characteristics taken from several sources. This paper examines correlations
between the aa index of the most severe geomagnetic storms and the intrinsic
characteristics of the sunspots from which they originated. We find a
correlation between the total rank of the aa index of the storms and the
“total intensity” of the sunspot, where total intensity is defined as the
sunspot’s mean intensity multiplied by its area. The correlation has an R2 =
0.690 and R2 = 0.855 when a potentially corrupted data point is removed.

Sunspots – Geomagnetic Storms – Coronal Mass Ejections – aa Index

##  1 Introduction

When the Earth’s magnetosphere is sufficiently disturbed, a geomagnetic storm
can potentially arise (Space Weather Prediction Center [Accessed September 23,
2021]). These disturbances mostly emerge from solar flares and Coronal Mass
Ejections (CMEs) which are oriented towards Earth. These CMEs are composed of
magnetized plasma and can form in minutes. When such structures occur, they
are violently discharged from the Sun and sent out into the solar system
(Temmer, M. 2021).

Such an event may give rise to a geomagnetic storm on Earth. The effects of a
geomagnetic storm can be devastating, such as causing entire nations’ power
grids to fail. However, geomagnetic storms seem to be very rare. Studies have
shown that space weather events that originate from the sun occur mostly
during the solar maximum phase (Hayakawa, H. et al., 2020); therefore, space
weather is not as prevalent during the solar minimum phase.

Strong geomagnetic storms may still occur during the solar minimum, such as
the geomagnetic storm of 1903 (Hayakawa, H. et al., 2020). A detailed study
was conducted by (Loewe, C.A. and G.W. Prölss, 1997) which found that the
number of great storms (those with $<$ -350 nT) occur once every 6 years while
3 strong ($<$ -100nT) and 2 severe ($<$ -200nT) storms occur every year. Due
to the large number of geomagnetic storms which occur annually, as well as the
occurrence of massive storms such as the Carrington event of 1859, it is
crucial to study these storms and their effects on Earth so that solar
physicists may more accurately predict when the next geomagnetic storm will
occur.

There are many factors to consider when predicting whether a sunspot will
flare up, such as its background magnetic field strength and spectra (Su, Y.
et al., 2007), (Lites, B.W. et al., 1993). Large and particularly complex
structured sunspots have been known to flare up in the past as well (Hayakawa,
H. et al., 2021). With the Solar Dynamics Observatory (SDO), we now have a
telescope in space that can see the sun in several key wavelengths to measure
and predict the sun’s activity. Yet new technology and theories about
predicting solar activity are still imperfect. Merely one geomagnetic storm as
severe as the Carrington event of 1859 is capable of causing the entire
world’s power grid to fail, which may result in trillions of dollars worth of
damage. The aim of this paper is to discover whether or not the ranks of the
geomagnetic storms of the aa index are correlated with values analyzed from
the sunspots from which these storms originated. If such a correlation exists,
this may give insight into methods that can predict if a sunspot, or group of
sunspots, will flare up and pose a threat to Earth. It is important to analyze
past events to acquire a better grasp on when the next major geomagnetic storm
will arise.

##  2 Procedure

In order to analyze sunspots prior to severe geomagnetic storms, images of the
sun taken 0-4 days before the occurrence of geomagnetic storms were required.
The images of the sun were acquired from two databases: the Astronomical
Photographic Data Archive (APDA) and the Kodaikanal Solar Observatory (KSO).
APDA is the second largest collection of astronomical plates in the world and
is housed at the Pisgah Astronomical Research Institute (PARI) in North
Carolina, United States. The images from both APDA and KSO consisted of
digital scans of emulsion plates; emulsion plates are thin glass squares which
have a negative image imprinted on it made from dried chemicals.

Approximately 0-4 days pass between coronal mass ejections leaving the sun and
arriving at earth (Temmer, M., 2021). Therefore, we analyzed the plates of the
sun corresponding to dates within that time window of the focus event. The
APDA emulsion plates were scanned using a Epson Expression 11000XL Photo
Scanner while the KSO plates were downloaded from KSO’s online database as
.TIFF files. Our project is concerned only with possible correlations between
Coronal Mass Ejections (CMEs) and solar flares with geomagnetic storms;
therefore, geomagnetic storms which were caused by phenomena other than CMEs
and solar flares were not studied.

Some emulsion plates taken prior to geomagnetic storms were in poor condition
or not available in APDA or KSO. Therefore, those geomagnetic storms were not
studied despite their status as some of the strongest geomagnetic storms ever
recorded. Some of the plates that were scanned and analyzed were 118 years old
at the time of digitization, and it cannot be assured that the plates are in
optimal condition. Once a list of relevant emulsion plates with associated
focus dates was compiled, we searched APDA and KSO for all the dates that were
0-4 days before known geomagnetic storms. In an effort to reduce randomized
error, only geomagnetic storms which corresponded to at least three good-
condition emulsion plates were studied; events in which only one or two good-
condition plates were available were not studied.

The aa index rank for a geomagnetic storm accounts for a variety of factors
including the peak aa value, storm duration, and number of peaks in the storm
interval (Vennerstrom et al. 2016). Because of the high number of variables,
definitively ranking the severity of a storm poses a challenge. The aa index
rankings are intended to serve primarily as a general guideline but not a
comprehensive analysis (Vennerstrom et al. 2016). The aa index contains data
beginning in 1868, while the Disturbance Storm Time (dst) records contain data
beginning in 1957. For our study, the aa index was chosen because our earliest
emulsion plates are from 1903, before the earliest dates in the dst (Hayakawa,
H. et al., 2021). In contrast, prior studies provide dst _estimates_ of
pre-1957 strong geomagnetic storms, such as the 1903 storm (Hayakawa, H. et
al., 2020), but there are no estimates for every storm before 1957.

The digitized KSO and APDA emulsion plates were all analyzed using the program
AstroImageJ. For each plate, the following values were recorded: the angle
between the sunspot and Earth relative to the center of the sun (_$\theta$_),
the apparent mean intensity and maximum intensity of the sunspot (Mean
_$I_{\theta}$_ and Max _$I_{\theta}$_ , respectively), the mean intensity and
maximum intensity of the sunspot with limb darkening taken into account (Mean
_$I_{0}$_ and Max _$I_{0}$_ , respectively, or “Mean Head-On Intensity” and
“Max Head-On Intensity,” respectively), the observed area of the sunspot
(_$A_{obs}$_), the area of the sunspot with foreshortening taken into account
(_$A_{0}$_ , or “head-on Area”), and the observed percentage of the solar
surface that the sunspot covered (_% Surf_).

This study also determined the orientation of the sun for each emulsion plate.
For each event, the relative location of the sunspot at the earliest date in
the event was compared to the relative location of the sunspot at the latest
date in the event. The sun rotates West to East as observed from earth;
therefore the proper orientation of the emulsion plates is the orientation in
which the sunspot’s relative location moves from West to East.

A calculation of the angle between the sunspot and Earth relative to the
center of the sun, _$\theta$_ , requires a knowledge of the pixel location
(_x,y_) of the center of the sun on the digitized emulsion plate. To find the
center of the sun, the Ellipse tool in AstroImageJ was used to fit a circle
precisely (to the nearest pixel) over the circumference of the sun. This
indicated two points on the circle which possessed the same x-coordinate as
the center of the sun and two points which possessed the same y-coordinate as
the center of the sun, allowing us to deduce the position ($x_{1}$, $y_{1}$)
of the center (Figure 1).

Figure 1: This is a solar plate from January $18^{th}$ 1938\. On the Western
and Eastern edges of the ellipse are two squares which indicate the _y_
-coordinate of the center of the sun. Similarly, the Northern and Southern
edges indicate the _x_ -coordinate of the center of the sun. The fitted circle
also provided the area and diameter of the sun in pixels. Using the known
radius of the sun $R_{Solar}$ = 1.3927*$10^{6}$ km, we established a rate of
conversion between pixels and kilometers. The position of the sunspot’s center
($x_{2}$, $y_{2}$) was also found with AstroImageJ.

The angle of the sunspot was calculated with the equation $\theta$ =
$sin^{-1}$(d/$R_{Solar}$), where d is the distance between the center of the
sun and the sunspot (with units of pixels) and $R_{Solar}$ is the radius of
the sun (in units of pixels).

The value of _d_ was calculated with the distance formula,

$d=\sqrt{(x_{1}-x_{2})^{2}+(y_{1}-y_{2})^{2}}.$ (1)

The observed maximum pixel intensity of the sunspot, observed mean pixel
intensity of the sunspot, and observed area of the sunspot were recorded using
the Maximum Value, Mean Value, and Area functions in AstroImageJ,
respectively.

The area of the sunspot with foreshortening taken into account was found using
the equation $dA_{obs}$ = $dA_{0}$ $\cos(\theta)$, where $dA_{obs}$ is an
infinitesimal two-dimensional segment of the observed area, $dA_{0}$ is an
infinitesimal segment of the area with foreshortening accounted for, and
$\theta$ is the angle between the sunspot and observer with respect to the
center of the sun (Figure 2). The sunspots were sufficiently small compared to
the surface area of the sun that $dA_{obs}$ $\approx$ $A_{obs}$ and $dA_{0}$
$\approx$ $A_{0}$. The observed percentage of the solar surface was calculated
by dividing the observed area of the sunspot in pixels by the area of the sun
in pixels. 2

Figure 2: The sunspot is positioned at an angle _$\theta$_ relative to the
observer. Therefore, the area covered in the digitized emulsion plate is less
than its true area, or “head-on area.” An infinitesimal area _$dA_{0}$_ on the
surface of the sun can be approximated as a flat surface. Because the area of
the sunspot is small compared to the area of the surface of the sun, we use
this approximation (to an accuracy of 10%). Note that the angle between
_$dA_{obs}$_ and _$dA_{0}$_ is equal to _$\theta$_.

One of the greatest challenges of this study was properly taking into account
limb darkening. Because the sun is a sphere and is of non-uniform density, the
center of the sun appears brighter to an observer than the edges, or limbs, of
the sun. This effect is known as limb darkening (Figure 3). Because the
sunspots closer to the limbs of the sun are surrounded by a darker background
than the sunspots closer to the center, a lack of taking limb darkening into
account would skew our results significantly.

Figure 3: A displacement-to-darkness graph across the diameter of the sun on
September 20th, 1909. The local maximum at approximately -1200 pixels
represents the sunspot. The observed intensity _$I_{\theta}$_ measures the
intensity of the sunspot relative to the x-axis (Pixel Darkness = 0). In
contrast, the head-on intensity _$I_{0}$_ measures the intensity of the
sunspot relative to the intensity of the solar surface that would otherwise be
at that displacement. The sharp peak at Displacement = 0 is due to a black
line drawn onto the emulsion plate to indicate the solar equator; it has no
physical significance.

To account for limb darkening, this study used the equation,

$\frac{I_{\theta}}{I_{0}}=a_{0}+a_{1}cos{\theta}+2a_{2}cos{\theta}^{2}$ (2)

This equation is taken from Introduction to stellar astrophysics, Böhm-Vitense
1989 (p. 233), where $I_{\theta}$ is the apparent brightness of a point on the
sun, $I_{0}$ is the brightness of the sun at the center, _$a_{0}$_ , _$a_{1}$_
, and _$2a_{2}$_ are empirically observed constants provided in (Böhm-Vitense
1989), and ${\theta}$ is the angle between the sunspot and the observer with
respect to the center of the sun. At 501 nm,

$({a_{0}},{a_{1}},{2a_{2}})=(0.2593,0.8724,-0.1136)$ (3)

This paper deals only with the relative brightness of pixels compared to each
other, and thus our equation is unitless. Because this study dealt with the
darkness of sunspots relative to the brightness of the sun, we modified our
equation to be

$\frac{I_{0}}{I_{\theta}}=a_{0}+a_{1}cos{\theta}+2a_{2}cos{\theta}^{2}$ (4)

so that $I_{0}$ and $I_{\theta}$ measure the pixel intensity of darkness, not
brightness. To avoid errors and discrepancies, every plate was analyzed by two
different people to make sure that the method being used was consistent. If
two people were getting completely different numbers, a third party was
brought in to analyze the plate.

##  3 Results and Discussion

Following the collection of all relevant data from the solar plates, we
attempted to find a correlation between sunspots and geomagnetic storms by
comparing certain values from each. For the geomagnetic storms, the primary
focus was on the overall rank of the storm, while a secondary focus was placed
on peak aa value. While the overall rank of the storm does not directly
correlate with any single observable value, we found it to be the best
indication of the strength of the storm as a whole. No singular value could
encompass all of the data. We also examined the peak aa value to determine
whether a correlation existed between the peak intensity of the sunspot and
the peak aa value.

For the sunspots, this paper focused primarily on the “total intensity,” here
defined as the area of the sunspot multiplied by its mean intensity. We
examined both the observed intensity and head-on intensity, as defined in
Section 2. We also examined the peak intensity of the sunspot, as it appeared
to be plausible that the peak intensity of the sunspot might correlate with
the peak aa index or perhaps the overall rank. For the values relating to the
sunspots, we chose the peak value for the entire event. If there were three
dates that were analyzed, only the day with the highest peak intensity or
overall intensity was used as this was most likely when the flare or CME
occured (See the Appendix, Section 6, for daily values).

When looking at the measured peak intensity of the sunspots, we graphed the
values for each event against both the peak aa index and the overall rank of
the geomagnetic storm. There was not a correlation between peak intensity and
either value. The peak intensity and overall rank did not correlate strongly
enough to warrant any confident conclusions about the relation between
sunspots and geomagnetic storms. The comparison of peak intensity to peak aa
index likewise had no statistically significant correlation. We note that
there existed several storms had approximately equal aa indices, but few
geomagnetic storms had approximately equal intensities. Indeed, the peak aa
index of geomagnetic storms correlated with virtually no other characteristic
which we examined.

$0$$2$$4$$6$$8$$10$$12$$14$$0$$2$$4$$6$$\cdot 10^{4}$RankPeak IntensityPeak
Intensity vs. Rank

$500$$550$$600$$650$$700$$750$$0$$2$$4$$6$$\cdot 10^{4}$aa PeakPeak
IntensityPeak Intensity vs. aa Peak

Figure 4: The two graphs above depict the relationship between the peak
intensity of sunspots and both the rank and aa peak of the resulting
geomagnetic storms. There is no statistically significant correlation present
in either graph.

The other values of sunspots which we compared to geomagnetic storms was the
total measured intensity of the sunspot. We modelled the overall intensity of
a sunspot by multiplying the mean intensity of the sunspot by its total area.
Both the overall intensity of the sunspot and the overall rank of the
geomagnetic storm provide an overall (but imperfect) indicator of their
significance, and thus we found it plausible that the two values might
correlate. We utilized three slightly different definitions of total
intensity, all of which we compared with the geomagnetic storm’s overall rank:
the “observed total intensity” ($A_{obs}$ * Mean $I_{obs}$), the “head-on
total intensity” ($A_{0}$ * Mean $I_{0}$), and a hybrid of the previous two
values ($A_{obs}$ * $I_{0}$). The highest total intensities for each sunspot
were graphed in relation to the overall ranks of the resulting geomagnetic
storms.

For each graph, both a linear and exponential fit was tested. An exponential
curve was a better fit for each of the graphs, a conclusion that is not
surprising for two main reasons. Firstly, the gaps in aa indices between the
storms’ ranks are unlikely to be uniform; e.g. the gap between medium-sized
geomagnetic storms should be smaller than the gap between large ones.
Secondly, a linear regression would eventually lead to negative intensity
values, which are not possible. As a result, an exponential fit was used to
test the correlation between total intensities and rank.

$0$$2$$4$$6$$8$$10$$12$$14$$0$$2$$4$$6$$\cdot 10^{4}$Total Rank$A_{0}$ * Mean
$I_{0}$$A_{0}$ * Mean $I_{0}$ vs. Total Rank$R^{2}$=0.388

$0$$2$$4$$6$$8$$10$$12$$14$$0$$2$$4$$6$$\cdot 10^{4}$Total Rank$A_{obs}$ *
Mean $I_{0}$$A_{obs}$ * Mean $I_{0}$ vs. Total Rank$R^{2}$=0.536

$0$$2$$4$$6$$8$$10$$12$$14$$0$$2$$4$$6$$\cdot 10^{4}$Total Rank$A_{obs}$ *
Mean $I_{obs}$$A_{obs}$ * Mean $I_{obs}$ vs. Total Rank$R^{2}$=0.690

Figure 5: The three graphs above depict the relationship between the overall
rank of the geomagnetic storms and several different methods for evaluating
the overall intensity of the preliminary sunspot. An exponential regression,
marked by the black line, is present in each of the graphs. The resulting
correlation between rank and each of the three measures of overall intensity
is depicted by an $R^{2}$ value.

As shown by the graphs, the total intensity measure that had the highest
correlation to the overall rank of the storm was the observed total intensity.
The $R^{2}$ value of the exponential fit was higher than the other two at
0.690. When looking at the graph of observed total intensity vs. rank, it is
clear that the rank six storm, which occurred in 1903, is somewhat of an
outlier. Because 1903 was the oldest date that we analyzed, it is possible
that those plates were in the worst condition and their data was somewhat
corrupted. It is also possible that the 1903 storm is an outlier for some
reason we have not found.

$0$$2$$4$$6$$8$$10$$12$$14$$0$$2$$4$$6$$\cdot 10^{4}$Total Rank$A_{obs}$ *
Mean $I_{obs}$$A_{obs}$ * Mean $I_{obs}$ vs. Total Rank$R^{2}$=0.855 Figure 6:
The graph above depicts the relationship between rank and the total intensity
value that had the best correlation, observed total intensity. In this graph,
however, the sixth ranked storm, an outlier, was removed to show the
correlation without it. It is possible that this graph shows a more accurate
relationship between rank and observed total intensity.

When that data point is removed from the graph, the correlation becomes even
stronger. With an $R^{2}$ value of 0.855, this exponential fit has a strong
correlation, suggesting a legitimate connection between rank and observed
total intensity. This correlation constitutes the main result of this work.

We find it conceptually unsurprising that the aa index rank correlated more
closely to the observed total intensity than the head-on total intensity or
hybrid model. Although the head-on values of intensity and area provide useful
information about the true nature of the sunspot, the geomagnetic storms which
resulted from the sunspots occurred here on Earth. Given that the Earth was
not always directly above the sunspot we analyzed, the strength of geomagnetic
storms would not result from the head-on intensity of the sunspots. Only a
fraction of the emissions would be pointed in the direction of Earth, and thus
only a fraction of the total intensity would be involved in causing a
geomagnetic storm. Therefore, using the observed characteristics of the
sunspot to find a correlation might provide a better representation of how it
affects Earth, while the head-on values of the sunspot are important for
characterizing it.

##  4 Conclusion

In this paper, we aimed to find a correlation between the rank of a
geomagnetic storm from the aa index to values analyzed from the sunspots they
came from. A correlation was found between the rank and the total intensity of
the sunspots ($A_{obs}$ * Mean $I_{obs}$) with an $R^{2}$ value of .690 when
including the outlier from 1903 and an $R^{2}$ 0.855 when excluding it.
Further insight on these storms may be gained from analyzing the physical
structure of the sunspots, spectroscopic analysis, and measuring its
background magnetic field strength. Analyzing more data of geomagnetic storms
that were caused by solar flares and CMEs would prove useful as there are some
dates that were missing. This would entail analyzing events that rank lower on
the overall aa index to see if it fits with the correlation we have found.

We conclude that the aa index rankings are a reasonable measurement of the
severity of the most significant geomagnetic storms, as the rankings correlate
well with our data and findings. One strong geomagnetic storm has the
potential to wipe out an entire country’s power grid, interrupt satellite
communication, and harm astronauts that are in orbit around the Earth. With
the world so heavily reliant on electronics, it may take decades to fully
recover from an event like the Carrington event, particularly if there is no
proper preparation beforehand. Continuing to study and analyze the connection
between solar activity and geomagnetic storms will be highly beneficial to us
on Earth. If humans are sent to other bodies in the solar system, such as
Mars, predicting geomagnetic storms will be especially important, as the
Martian atmosphere is thinner than on Earth and therefore space weather may
cause significant harm to the astronauts and electrical components.

A possible future study may involve converting our unitless data in the
Appendix (Section 6) to metric units or cgs units. A knowledge of the exact
flux of the sunspots relative to the solar surface may yield insightful
results. Future studies may also analyze images of the sun corresponding to
dates that this study could not analyze due to broken or missing plates in
APDA or KSO.

##  5 Acknowledgements

We thank Thurburn Barker for being the advisor for this project and helping
with everything involving APDA. We also would like to thank Sarah Gibson from
the High Altitude Observatory at NCAR who helped answer many important
questions and provided excellent insights. Thank you to the Pisgah
Astronomical Research Institute (PARI) for giving us all the resources
necessary to complete this project. Thank you to the Kodaikanal Solar
Observatory (KSO) for having a free open database that helped supply the data
that was not available in APDA.

##  6 Appendix

Figure 7: The data used in the figures in Section 3 is tabulated here. The
data was acquired from digital scans of emulsion plates of the sun from the
APDA and KSO databases. The asterisk in the data for 5/10/21 indicates that
the sunspot had split into two, and thus presented a difficulty in accurately
representing its area.

##  7 References

Alexander, D., Richardson, I.G., $\&$ Zurbuchen, T.H. 2006, A Brief History of
CME Science. In: Coronal Mass Ejections. Space Sciences Series of ISSI, vol
21. Springer, New York, NY. doi: https://doi.org/10.1007/978-0-387-45088-9_1

Böhm-Vitense, E. (1989) Introduction to Stellar Astrophysics. Cambridge:
Cambridge University Press. doi: https://doi.org/10.1017/CBO9780511623011

Boyd, R.W. (1983). Radiometry and the detection of optical radiation. New
York: John Wiley And Sons.

Hayakawa, H. et al., 2020. The extreme space weather event in 1903
October/November: An outburst from the Quiet Sun. The Astrophysical Journal,
897(1). doi: 10.3847/2041-8213/ab6a18

Hayakawa, H. et al., 2021. The intensity and evolution of the Extreme solar
and geomagnetic storms in 1938 January. The Astrophysical Journal, 909(2),
p.197. doi: 10.3847/1538-4357/abc427

Lites, B.W. et al., 1993. Stokes profile analysis and Vector magnetic Fields.
VI. fine scale structure of a sunspot. The Astrophysical Journal, 418, p.9.
doi: 10.1086/173450

Loewe, C.A. and G.W. Prölss, 1997. Classification and mean behavior of
magnetic storms. Journal of Geophysical Research: Space Physics, 102(A7),
pp.14209–14213. doi: https://doi.org/10.1029/96JA04020

Space Weather Prediction Center. Accessed 23 September 2021. Geomagnetic
storms. National Oceanic and Atmospheric Administration. Available at:
https://www.swpc.noaa.gov/phenomena/geomagnetic-storms.

Su, Y. et al., 2007. What determines the intensity of solar flare/cme events?
The Astrophysical Journal, 665(2), pp.1448–1459. doi: 10.1086/519679

Temmer, M. Space weather: The Solar Perspective. Living Rev Sol Phys 18, 4
(2021). doi: https://doi.org/10.1007/s41116-021-00030-3

Vennerstrom, S. et al., 2016. Extreme geomagnetic storms – 1868 – 2010. Solar
Physics, 291(5), pp.1447–1481. doi: 10.1007/s11207-016-0897-y