# SHED: A Newton-type algorithm for federated learning based on incremental
Hessian eigenvector sharing

Nicolò Dal Fabbro<EMAIL_ADDRESS>
Department of Information Engineering
University of Padova
Padova, PD 35122, Italy Subhrakanti Dey<EMAIL_ADDRESS>
Department of Electrical Engineering
Uppsala University
Uppsala, 751 03, Sweden Michele Rossi<EMAIL_ADDRESS>
Department of Information Engineering
University of Padova
Padova, PD 35122, Italy Luca Schenato<EMAIL_ADDRESS>
Department of Information Engineering
University of Padova
Padova, PD 35122, Italy

###### Abstract

There is a growing interest in the distributed optimization framework that
goes under the name of Federated Learning (FL). In particular, much attention
is being turned to FL scenarios where the network is strongly heterogeneous in
terms of communication resources (e.g., bandwidth) and data distribution. In
these cases, communication between local machines (_agents_) and the central
server (_Master_) is a main consideration. In this work, we present SHED111In
a previous preprint version of this paper, the algorithm was referred to as
IOS, an original communication-constrained Newton-type (NT) algorithm designed
to accelerate FL in such heterogeneous scenarios. SHED is by design robust to
non i.i.d. data distributions, handles heterogeneity of agents’ communication
resources (CRs), only requires sporadic Hessian computations, and achieves
super-linear convergence. This is possible thanks to an _incremental_
strategy, based on eigendecomposition of the local Hessian matrices, which
exploits (possibly) _outdated_ second-order information. The proposed solution
is thoroughly validated on real datasets by assessing (i) the number of
communication rounds required for convergence, (ii) the overall amount of data
transmitted and (iii) the number of local Hessian computations. For all these
metrics, the proposed approach shows superior performance against state-of-the
art techniques like GIANT and FedNL.

Keywords: Federated learning, Newton method, distributed optimization,
heterogeneous networks, edge learning, non i.i.d. data.

##  1 Introduction

With the growing computational power of edge devices and the booming increase
of data produced and collected by users worldwide, solving machine learning
problems without having to collect data at a central server is becoming very
appealing (Shi et al., 2020). One of the main reasons for not transferring
users’ data to cloud central servers is due to privacy concerns. Indeed,
users, such as individuals or companies, may not want to share their private
data with other network entities, while training their machine learning
algorithms. In addition to privacy, distributed processes are by nature more
resilient to node/link failures and can be directly implemented on the servers
at network edge, i.e., within multi-access edge computing (MEC) scenarios
(Pham et al., 2020).

This distributed learning framework goes under the name of federated learning
(FL) (McMahan et al., 2017; Li et al., 2020a), and it has attracted much
research interest in recent years. Direct applications of FL can be found, for
example, in the field of healthcare systems (Rieke et al., 2020; Huang et al.,
2020) or of smartphone utilities (Hard et al., 2018).

One of the most popular FL algorithms is federated averaging (FedAvg) (McMahan
et al., 2017), which achieves good results, but which comes with weak
convergence guarantees and has been shown to diverge under heterogeneous data
distributions (Li et al., 2020b).

Indeed, among the open challenges of FL, a key research question is how to
provide efficient distributed optimization algorithms in scenarios with
heterogeneous communication links (different bandwidth) and non i.i.d. data
distributions (Kairouz et al., 2021; Zhu et al., 2021; Smith et al., 2017;
Zhao et al., 2018). These aspects are found in a variety of applications, such
as the so-called _federated edge learning_ framework, where learning is moved
to the network edge, and which often involves unstable and heterogeneous
wireless connections (Shi et al., 2020; T. Dinh et al., 2022; Nguyen et al.,
2021; Chen et al., 2020). In addition, the MEC and FL paradigms are of high
interest to IoT scenarios such as data retrieval and processing within smart
cities, which naturally entail non i.i.d. data distributions due to inherent
statistical differences in the underlying spatial processes (e.g., vehicular
mobility, user density, etc.) (Liu et al., 2020).

The bottleneck represented by communication overhead is one of the most
critical aspects of FL. In fact, in scenarios with massive number of devices
involved, inter-agent communication can be much slower than the local
computations performed by the FL agents themselves (e.g., the edge devices),
by many orders of magnitude (Li et al., 2020a). The problem of reducing the
communication overhead of FL becomes even more critical when the system is
characterized by non i.i.d. data distributions and heterogeneous communication
resources (CRs) (Li et al., 2020b). In this setting, where the most critical
aspect is inter-node communications, FL agents are usually assumed to have
good computing capabilities, and wisely increasing the computation effort at
the agents is a good strategy to obtain a faster convergence. For this reason,
Newton-type approaches, characterized by strong robustness and fast
convergence rates, even if computationally demanding, have been recently
advocated (Gupta et al., 2021; Safaryan et al., 2021).

In the present work, we present SHED (Sharing Hessian Eigenvectors for
distributed learning), a Newton-type FL framework that is naturally robust to
non i.i.d. data distributions and that intelligently allocates the (per-
iteration) communication resources of the involved FL agents, allowing those
with more CRs to contribute more towards improving the convergence rate.
Differently from prior art, SHED requires FL agents to locally compute the
Hessian matrix sporadically. Its super-linear convergence is here proven by
studying the CRs-dependent convergence rate of the algorithm by analyzing the
dominant Lypaunov exponent of the estimation error. With respect to other
Newton-type methods for distributed learning, our empirical results
demonstrate that the proposed framework is (i) competitive with state-of-the-
art approaches in i.i.d. scenarios and (ii) robust to non i.i.d. data
distributions, for which it outperforms competing solutions.

###  1.1 Related work

Next, to put our contribution into context, we review previous related works
on first and second order methods for distributed learning.

First-order methods. Some methods have been recently proposed to tackle non
i.i.d. data and heterogeneous networks. SCAFFOLD (Karimireddy et al., 2020),
FedProx (Li et al., 2020b) and the work in Yu et al. (2019) propose
modifications to the FedAvg algorithm to face system’s heterogeneity. To deal
with non i.i.d. datasets, strategies like those in Zhao et al. (2018) and
Jeong et al. (2018) have been put forward, although these approaches require
that some data is shared with the central server, which does not fit the
privacy requirements of FL. To reduce the FL communication overhead, the work
in Chen et al. (2018) leverages the use of outdated first-order information by
designing simple rules to detect slowly varying local gradients. With respect
to heterogeneous time-varying CRs problem, Amiri and Gündüz (2020) presented
techniques based on gradient quantization and on analog communication
strategies exploiting the additive nature of the wireless channel, while Chen
et al. (2020) studied a framework to jointly optimize learning and
communication for FL in wireless networks.

Second-order methods. Newton-type (NT) methods exploit second-order
information of the cost function to provide accelerated optimization, and are
therefore appealing candidates to speed up FL. NT methods have been widely
investigated for distributed learning purposes: GIANT (Wang et al., 2018) is
an NT approach exploiting the harmonic mean of local Hessian matrices in
distributed settings. Other related techniques are LocalNewton (Gupta et al.,
2021), DANE (Shamir et al., 2014), AIDE (Reddi et al., 2016), DiSCO (Zhang and
Lin, 2015), DINGO (Crane and Roosta, 2019) and DANLA (Zhang et al., 2022).
DONE (T. Dinh et al., 2022) is another technique inspired by GIANT and
specifically designed to tackle federated edge learning scenarios.
Communication efficient NT methods like GIANT and DONE exploit an extra
communication round to obtain estimates of the global Hessian from the
harmonic mean of local Hessian matrices. These algorithms, however, were all
designed assuming that data is i.i.d. distributed across agents and, as we
empirically show in this paper, under-perform if such assumption does not
hold. A recent study, FedNL (Safaryan et al., 2021), proposed algorithms based
on matrix compression which use theory developed in Islamov et al. (2021) to
perform distributed training, by iteratively learning the Hessian matrix at
the optimum. A similar approach based on matrix compression techniques has
been considered in Basis Learn (Qian et al., 2021), proposing variants of
FedNL that exploit change of basis in the space of matrices to improve the
quality of the compressed Hessian. However, these solutions do not consider
heterogeneity in the CRs, and require the computation of the local Hessian at
each iteration. Another NT technique has been recently proposed in Liu et al.
(2021), based on the popular L-BFGS (Liu and Nocedal, 1989) quasi-Newton
method. Other recent related NT approaches have been proposed in FLECS
(Agafonov et al., 2022), FedNew (Elgabli et al., 2022) and Quantized Newton
(Alimisis et al., 2021). Although some preliminary NT approaches have been
proposed for the FL framework, there is still large space for improvements
especially in ill-conditioned setups with heterogeneous CRs and non i.i.d.
data distributions.

###  1.2 Contribution

In this work, we design SHED (Sharing Hessian Eigenvectors for Distributed
learning), a communication-constrained NT optimization algorithm for FL that
is suitable for networks with heterogeneous CRs and non i.i.d. data. SHED
requires only sporadic Hessian computation and it is shown to accelerate
convergence of convex federated learning problems, reducing communication
rounds and communication load in networks with heterogeneous CRs and non
i.i.d. data distributions.

In SHED, local machines, the _agents_ , compute the local Hessian and its
eigendecomposition. Then, according to their available communication resources
(CRs), agents share some of the most _relevant_ Hessian eigenvector-eigenvalue
pairs (EEPs) with the server (the Master). Together with the EEPs, agents send
to the Master a scalar quantity - computed locally - that allows a full-rank
Hessian approximation. The obtained approximations are averaged at the Master
in order to get an approximated global Hessian, to be used to perform global
NT descent steps. The approximation of the Hessian matrix at the Master is
then improved by _incrementally_ sending additional EEPs. In particular, for a
general convex problem, SHED is designed to make use of _outdated_ Hessian
approximations. Thanks to this, local Hessian computation is required only
sporadically. We analyze the convergence rate via the analysis of the dominant
Lyapunov exponent of the objective parameter estimation error. In this novel
framework, we show that the use of outdated Hessian approximations provides an
improvement in the convergence rate. Furthermore, we show that the EEPs of the
outdated local Hessian approximations can be shared heterogeneously with the
Master, so the algorithm is versatile with respect to (per-iteration)
heterogeneous CRs. We also show that, thanks to the proposed incremental
strategy, SHED enjoys a super-linear convergence rate.

###  1.3 Organization of the paper

The rest of the paper is organized as follows: in Section 2, we detail the
problem formulation and the general idea behind SHED design. In Section 3, we
present SHED and the theoretical results in the linear regression (least
squares) case, describing the convergence rate with respect to the algorithm
parameters. Section 4 is instrumental to extending the algorithm to a general
strongly convex problem, which is done later on in Section 5. At the end of
Section 5, based on the theoretical results, we propose some heuristic choices
for tuning SHED parameters. Finally, in Section 6, empirical performance of
SHED are shown on real datasets. In particular, the theoretical results are
validated and compared against those from state-of-the art approaches. Longer
proofs and some additional experiments are reported in the Appendices, for
completeness.

_Notation:_ Vectors and matrices are written as lower and upper case bold
letters, respectively (e.g., vector $\mathbf{v}$ and matrix $\mathbf{V}$). The
operator $\|\cdot\|$ denotes the $2$-norm for vectors and the spectral norm
for matrices. $diag(\mathbf{v})$ denotes a diagonal matrix with the components
of vector $\mathbf{v}$ as diagonal entries. $\mathbf{I}$ denotes the identity
matrix.

##  2 Problem formulation

In Fig. 1 we illustrate the typical framework of a federated learning
scenario. In the rest of the paper, we denote the optimization parameter by
$\boldsymbol{\theta}\in\mathbb{R}^{n}$, and a generic cost function by
$f:\mathbb{R}^{n}\rightarrow\mathbb{R}$, where $n$ is the dimension of the
feature data vectors $\mathbf{x}\in\mathbb{R}^{n}$.
Let us denote a dataset by $\mathcal{D}=\\{\mathbf{x}_{j},y_{j}\\}_{j=1}^{N}$,
where $N$ is the global number of data samples.
$\mathbf{x}_{j}\in\mathbb{R}^{n}$ denotes the $j$-th data sample and
$y_{j}\in\mathbb{R}$ the response of sample $j$. In the case of classification
problems, $y_{j}$ would be an integer specifying the class to which sample
$\mathbf{x}_{j}$ belongs.

Figure 1: The typical FL optimization framework: agents $A_{1},...,A_{M}$
cooperate to solve a common learning problem. After receiving the global
parameter $\operatorname*{{\boldsymbol{\theta}}}^{t}$, at the current
iteration $t$, they share their optimization sets
$U_{t}^{(1)},...,U_{t}^{(M)}$ so that a new global update can be performed at
the Master.

We consider the problem of regularized empirical risk minimization of the
form:

$\min_{\boldsymbol{\theta}}f(\boldsymbol{\theta}),$ (1)

with

$f(\boldsymbol{\theta}):=\frac{1}{N}\sum_{j=1}^{N}l_{j}(\boldsymbol{\theta})+\frac{\mu}{2}\|\boldsymbol{\theta}\|_{2}^{2},$

where each $l_{j}(\boldsymbol{\theta})$ is a convex function related to the
$i$-th element of $\mathcal{D}$ and $\mu$ a regularization parameter. In
particular,

$l_{j}(\boldsymbol{\theta})=l(\mathbf{x}_{j}^{T}\boldsymbol{\theta},y_{j}),$

where $l:\mathbb{R}\rightarrow\mathbb{R}$ is a convex function. Examples of
convex cost functions considered in this work are:

  * •

linear regression with quadratic cost (least squares):

$l(\mathbf{x}_{j}^{T}\boldsymbol{\theta},y_{j})=\frac{1}{2}(\mathbf{x}_{j}^{T}\boldsymbol{\theta}-y_{j})^{2}$

  * •

logistic regression:

$l(\mathbf{x}_{j}^{T}\boldsymbol{\theta},y_{j})=\log(1+e^{-y_{j}(\mathbf{x}_{j}^{T}\boldsymbol{\theta})})$

Logistic regression belongs to the class of generalized linear models. In
general, in this paper we consider optimization problems in which the
following assumptions on the agents’ cost functions $f^{(i)}$ hold:

###### Assumption 1

Let
$\mathbf{H}^{(i)}(\boldsymbol{\theta}):=\nabla^{2}f^{(i)}(\boldsymbol{\theta})$
the local Hessian matrix of the cost
$f^{(i)}(\operatorname*{{\boldsymbol{\theta}}})$ of agent $i$.
$f^{(i)}(\operatorname*{{\boldsymbol{\theta}}})$ is twice continuously
differentiable, $K_{i}$-smooth, $\kappa_{i}$-strongly convex and
$\mathbf{H}^{(i)}(\boldsymbol{\theta})$ is $L_{i}$-Lipschitz continuous for
all $i=1,...,M$.

The above assumptions imply that

$\kappa_{i}\mathbf{I}\leq\mathbf{H}^{(i)}(\boldsymbol{\theta})\leq
K_{i}\mathbf{I},\ \forall{\boldsymbol{\theta}}$
$\|\mathbf{H}^{(i)}(\boldsymbol{\theta})-\mathbf{H}^{(i)}(\boldsymbol{\theta}^{\prime})\|\leq
L_{i}\|\boldsymbol{\theta}-\boldsymbol{\theta}^{\prime}\|,\
\forall{\boldsymbol{\theta},\boldsymbol{\theta}^{\prime}}$

Note that $K_{i}$-smoothness is also a consequence of the $K$-smoothness of
any global cost function $f$ for some constant $K$, while strong convexity of
agents’ cost functions is always guaranteed if the functions
$l_{j}(\operatorname*{{\boldsymbol{\theta}}})$ are convex in the presence of a
regularization term $\mu>0$.

In the paper, we denote the data matrix by
$\mathbf{X}=[\mathbf{x}_{1},...,\mathbf{x}_{N}]\in\mathbb{R}^{n\times N}$, and
the label vector as $\mathbf{Y}=[y_{1},...,y_{N}]\in\mathbb{R}^{1\times N}$.
We refer to the dataset $\mathcal{D}$ as the tuple
$\mathcal{D}=(\mathbf{X},\mathbf{Y})$. We denote by $M$ the number of agents
involved in the optimization algorithm, and we can write
$\mathbf{X}=[\mathbf{X}_{1},...,\mathbf{X}_{M}]$,
$\mathbf{Y}=[\mathbf{Y}_{1},...,\mathbf{Y}_{M}]$, where
$\mathbf{X}_{i}\in\mathbb{R}^{n\times N_{i}}$ and
$\mathbf{Y}_{i}\in\mathbb{R}^{1\times N_{i}}$, where $N_{i}$ is the number of
data samples of the $i$-th agent. We denote the local dataset of agent $i$ by
$\mathcal{D}_{i}=(\mathbf{X}_{i},\mathbf{Y}_{i})$.

###  2.1 A Newton-type method based on eigendecomposition

The Newton method to solve (1) works as follows (Newton update):

$\boldsymbol{\theta}^{t+1}=\boldsymbol{\theta}^{t}-\eta_{t}\mathbf{H}_{t}^{-1}\mathbf{g}_{t},$

where $t$ denotes the $t$-th iteration,
$\mathbf{g}_{t}=\mathbf{g}(\boldsymbol{\theta}_{t})=\nabla
f(\boldsymbol{\theta}^{t})$ is the gradient at iteration $t$ and $\eta_{t}$ is
the step size (also called learning rate) at iteration $t$.
$\mathbf{H}_{t}=\nabla^{2}f(\boldsymbol{\theta}^{t})$ denotes the Hessian
matrix at iteration $t$. Compared to gradient descent, the Newton method
exploits the curvature information provided by the Hessian matrix to improve
the descent direction. We define
$\mathbf{p}_{t}:=\mathbf{H}_{t}^{-1}\mathbf{g}_{t}$. In general, Newton-type
(NT) methods try to get an approximation of $\mathbf{p}_{t}$. In an FL
scenario, assuming for simplicity that all $M$ agents have the same amount of
data, we have that:

$\mathbf{H}_{t}=\frac{1}{M}\sum_{i=1}^{M}\mathbf{H}_{t}^{(i)},\ \
\mathbf{g}_{t}=\frac{1}{M}\sum_{i=1}^{M}\mathbf{g}_{t}^{(i)},$ (2)

where $\mathbf{H}_{t}^{(i)}=\nabla^{2}f^{(i)}(\boldsymbol{\theta}^{t})$ and
$\mathbf{g}_{t}^{(i)}=\nabla f^{(i)}(\boldsymbol{\theta}^{t})$ denote local
Hessian and gradient of the local cost $f^{(i)}(\boldsymbol{\theta}^{t})$ of
agent $i$, respectively. To get a Newton update at the master, in an FL
setting one would need each agent to transfer the whole matrix
$\mathbf{H}_{t}^{(i)}$ of size $O(n^{2})$ to the master at each iteration,
that is considered a prohibitive communication complexity in a federated
learning setting, especially when the size of the feature data vectors, $n$,
increases. Because of this, in this work we propose an algorithm in which the
Newton-type update takes the form:

$\boldsymbol{\theta}^{t+1}=\boldsymbol{\theta}^{t}-\eta_{t}\hat{\mathbf{H}}_{t}^{-1}\mathbf{g}_{t},$
(3)

where $\hat{\mathbf{H}}_{t}$ is an approximation of $\mathbf{H}_{t}$. In some
previous contributions, like Wang et al. (2018); T. Dinh et al. (2022),
$\hat{\mathbf{H}}_{t}$ is the harmonic mean of local Hessians, obtained at the
master through an intermediate additional communication round to provide each
agent with the global gradient $\mathbf{g}_{t}$. In the recent FedNL algorithm
(Safaryan et al., 2021), each agent at each iteration sends a compressed
version of a Hessian-related matrix. In the best performing version of the
algorithm, the rank-1 approximation using the first eigenvector of the matrix
is sent. Authors show that thanks to this procedure they can eventually get
the Hessian matrix at the optimum, which provides super-linear convergence
rate.

Instead, our proposal is to incrementally obtain at the master an average of
full-rank approximations of local Hessians through the communication of the
local Hessian most relevant eigenvectors together with a carefully computed
local approximation parameter. In particular, we approximate the Hessian
matrix $\mathbf{H}_{t}$ exploiting eigendecomposition in the following way:
the symmetric positive definite Hessian $\mathbf{H}_{t}$ can be diagonalized
as $\mathbf{H}_{t}=\mathbf{V}_{t}\boldsymbol{\Lambda}_{t}\mathbf{V}_{t}^{T}$,
with $\boldsymbol{\Lambda}_{t}=diag(\lambda_{1,t},...,\lambda_{n,t})$, where
$\lambda_{k,t}$ is the eigenvalue corresponding to the $k$-th eigenvector,
$\mathbf{v}_{k,t}$. $\mathbf{H}_{t}$ can be approximated as

$\hat{\mathbf{H}}_{t}:=\hat{\mathbf{H}}_{t}(\rho_{t},q_{t})=\mathbf{V}_{t}\hat{\boldsymbol{\Lambda}}_{t}\mathbf{V}_{t}^{T}=\sum_{k=1}^{q_{t}}(\lambda_{k,t}-\rho_{t})\mathbf{v}_{k,t}\mathbf{v}_{k,t}^{T}+\rho_{t}\mathbf{I},$
(4)

with
$\hat{\boldsymbol{\Lambda}}_{t}:=\hat{\boldsymbol{\Lambda}}_{t}(\rho_{t},q_{t})=diag(\lambda_{1,t},...,\lambda_{q_{t}},\rho_{t},...,\rho_{t})$.
The scalar $\rho_{t}>0$ is the approximation parameter (if $\rho_{t}=0$ this
becomes a low-rank approximation). The integer $q_{t}=1,...,n$ denotes the
number of eigenvalue-eigenvector pairs (EEPs)
$\\{\lambda_{k,t},\mathbf{v}_{k,t}\\}_{k=1}^{q_{t}}$ being used to approximate
the Hessian matrix. We always consider eigenvalues ordered so that
$\lambda_{1,t}\geq\lambda_{2,t}\geq...\geq\lambda_{n,t}$. The approximation of
Eqn. (4) was used in Erdogdu and Montanari (2015) for a subsampled centralized
optimization problem. There, the parameter $\rho_{t}$ was chosen to be equal
to $\lambda_{q_{t}+1}$. In the FL setting, we use approximation shown in Eqn.
(4) to approximate the local Hessian matrices of the agents. In particular,
letting $\hat{\mathbf{H}}_{t}^{(i)}$ be the approximated local Hessian of
agent $i$, the approximated global Hessian is the average of the local
approximated Hessian matrices:

$\hat{\mathbf{H}}_{t}=\sum_{i=1}^{M}{p_{i}\hat{\mathbf{H}}_{t}^{(i)}},\ \
p_{i}=N_{i}/N$

where
$\hat{\mathbf{H}}_{t}^{(i)}:=\hat{\mathbf{H}}_{t}(\rho_{t}^{(i)},q_{t}^{(i)})$
is a function of the local approximation parameter $\rho_{t}^{(i)}$ and of the
number of eigenvalue-eigenvector pairs $q_{t}^{(i)}$ shared by agent $i$,
which are denoted by
$\\{\lambda_{k,t}^{(i)},\mathbf{v}_{k,t}^{(i)}\\}_{k=1}^{q_{t}^{(i)}}$.

###  2.2 The algorithm in a nutshell

The idea of the SHED algorithm is that agents share with the Master, together
with the gradient, some of their local Hessian EEPs, according to the
available CRs. They share the EEPs in a decreasing order dictated by the value
of the positive eigenvalues corresponding to the eigenvectors. At each
iteration, they incrementally add new EEPs to the information they have sent
to the Master. In a linear regression problem, in which the Hessian does not
depend on the current parameter, agents would share their EEPs incrementally
up to the $n$-th. When the $n$-th EEP is shared, the Master has the full
Hessian available and no further second order information needs to be
transmitted. In a general convex problem, in which the Hessian matrix changes
at each iteration, as it is a function of the current parameter, SHED is
designed in a way in which agents perform a _renewal_ operation at certain
iterations, i.e., they re-compute the Hessian matrix and re-start sharing the
EEPs from the most relevant ones of the new matrix.

##  3 Linear regression (least squares)

In this section we illustrate our algorithm and present the convergence
analysis considering the problem of solving (1) via Newton-type updates (3) in
the least squares (LSs) case, i.e., in the case of linear regression with
quadratic cost. Before moving to the FL case, we prove some results for the
convergence rate of the centralized iterative least squares problem.
First, we provide some important definitions. In the case of linear regression
with quadratic cost, the Hessian is such that
$\mathbf{H}_{LS}:=\mathbf{H}(\boldsymbol{\theta}),\
\forall{\boldsymbol{\theta}}$, i.e., the Hessian does not depend on the
parameter $\boldsymbol{\theta}$. Because of this, when considering LSs, we
write the eigendecomposition as
$\mathbf{H}_{LS}=\mathbf{V}\boldsymbol{\Lambda}\mathbf{V}^{T}$, with
$\boldsymbol{\Lambda}=diag(\lambda_{1},...,\lambda_{n})$ without specifying
the iteration $t$ when not needed. Let $\boldsymbol{\theta}^{*}$ denote the
solution to (1). In the following, we use the fact that the cost for parameter
$\boldsymbol{\theta}^{t}$ can be written as
$f(\boldsymbol{\theta}^{t})=f(\boldsymbol{\theta}^{*})+\bar{f}(\boldsymbol{\theta}^{t})$,
with
$\bar{f}(\boldsymbol{\theta}^{t}):=\frac{1}{2}(\boldsymbol{\theta}^{t}-\boldsymbol{\theta}^{*})^{T}\mathbf{H}_{LS}(\boldsymbol{\theta}^{t}-\boldsymbol{\theta}^{*})$.
Similarly, the gradient can be written as
$\mathbf{g}_{t}=\mathbf{H}_{LS}(\boldsymbol{\theta}^{t}-\boldsymbol{\theta}^{*})$.
In this setup, the update rule of Eqn. (3) can be written as a time-varying
linear discrete-time system:

${\operatorname*{{\boldsymbol{\theta}}}}^{t+1}-{\operatorname*{{\boldsymbol{\theta}}}}^{*}=\mathbf{A}_{t}({\operatorname*{{\boldsymbol{\theta}}}}^{t}-{\operatorname*{{\boldsymbol{\theta}}}}^{*})$
(5)

where

$\mathbf{A}_{t}:=\mathbf{A}(\rho_{t},\eta_{t},q_{t})=\mathbf{I}-\eta_{t}\hat{\mathbf{H}}_{t}^{-1}\mathbf{H}_{LS}.$

Indeed,

$\displaystyle{\operatorname*{{\boldsymbol{\theta}}}}^{t+1}-{\operatorname*{{\boldsymbol{\theta}}}}^{*}$
$\displaystyle={\operatorname*{{\boldsymbol{\theta}}}}^{t}-{\operatorname*{{\boldsymbol{\theta}}}}^{*}-\eta_{t}\hat{\mathbf{H}}_{t}^{-1}\mathbf{g}_{t}$
$\displaystyle=\boldsymbol{\theta}^{t}-\boldsymbol{\theta}^{*}-\eta_{t}\hat{\mathbf{H}}_{t}^{-1}\mathbf{H}_{LS}(\boldsymbol{\theta}^{t}-\boldsymbol{\theta}^{*})$
$\displaystyle=(\mathbf{I}-\eta_{t}\hat{\mathbf{H}}_{t}^{-1}\mathbf{H}_{LS})(\boldsymbol{\theta}^{t}-\boldsymbol{\theta}^{*}).$

###  3.1 Centralized iterative least squares

In this sub-section, we study the optimization problem in the centralized
case, so when all the data is kept in a single machine. We provide a range of
choices for the approximation parameter $\rho_{t}$ that are optimal in the
convergence rate sense. We denote the convergence factor of the descent
algorithm described by Eqn. (3) by

$r_{t}:=r(\rho_{t},\eta_{t},q_{t}),$

making its dependence on the tuple $(\rho_{t},\eta_{t},q_{t})$ explicit.

###### Theorem 1

Consider solving problem (1) via Newton-type updates (3) in the least squares
case. At iteration $t$, let the Hessian matrix $\mathbf{H}_{LS}$ be
approximated as in Eqn. (4) (centralized case). The convergence rate is
described by

$\|\boldsymbol{\theta}^{t+1}-\boldsymbol{\theta}^{*}\|\leq
r_{t}\|\boldsymbol{\theta}^{t}-\boldsymbol{\theta}^{*}\|.$ (6)

For a given $q_{t}\in\\{0,1,...,n\\}$ the best achievable convergence factor
is

$r_{t}^{*}=r^{*}(q_{t}):=\min_{(\rho_{t},\eta_{t})}{r(\rho_{t},\eta_{t},q_{t})}=(1-\frac{\lambda_{n}}{\rho^{*}_{t}}),$
(7)

where

$\rho_{t}^{*}:=(\lambda_{q_{t}+1}+\lambda_{n})/2.$ (8)

$r_{t}^{*}$ is achievable if and only if
$(\rho_{t},\eta_{t})\in\mathcal{S}^{*}$, with

$\mathcal{S}^{*}=\\{(\rho_{t},\eta_{t}):\rho_{t}\in[\lambda_{n},\lambda_{q_{t}+1}],\eta_{t}^{*}=\frac{2\rho_{t}}{\lambda_{q_{t}+1}+\lambda_{n}}\\}$
(9)

Proof  From (4), writing
$\mathbf{H}_{LS}=\mathbf{V}\boldsymbol{\Lambda}\mathbf{V}^{T}$ and
$\hat{\mathbf{H}}_{t}=\mathbf{V}\hat{\boldsymbol{\Lambda}}_{t}\mathbf{V}^{T}$,
with
$\hat{\boldsymbol{\Lambda}}_{t}=diag(\lambda_{1},...,\lambda_{q_{t}},\rho_{t},...,\rho_{t})$,
define
$\boldsymbol{\theta}_{\rho_{t},\eta_{t}}^{t+1}:=\boldsymbol{\theta}^{t}-\eta_{t}{\hat{\mathbf{H}}_{t}}^{-1}\mathbf{g}_{t}$.
Recalling that
$\mathbf{g}_{t}=\mathbf{H}_{LS}(\boldsymbol{\theta}^{t}-\boldsymbol{\theta}^{*})$,
we have:

$\displaystyle\boldsymbol{\theta}_{\rho_{t},\eta_{t}}^{t+1}-\boldsymbol{\theta}^{*}=\mathbf{A}_{t}(\boldsymbol{\theta}^{t}-\boldsymbol{\theta}^{*})$
$\displaystyle=(\mathbf{I}-\eta_{t}\hat{\mathbf{H}}_{t}^{-1}\mathbf{H}_{LS})(\boldsymbol{\theta}^{t}-\boldsymbol{\theta}^{*})$
(10)
$\displaystyle=\mathbf{V}(\mathbf{I}-\eta_{t}\hat{\boldsymbol{\Lambda}}_{t}^{-1}\boldsymbol{\Lambda})\mathbf{V}^{T}(\boldsymbol{\theta}^{t}-\boldsymbol{\theta}^{*}).$

For some given $q_{t}\in\\{1,...,n\\}$, $r_{t}$ is a function of two tunable
parameters, i.e., the tuple $(\eta_{t},\rho_{t})$. We now prove that
$r_{t}^{*}$ can be achieved if and only if
$\rho_{t}\in[\lambda_{n},\lambda_{q_{t}+1}]$. The convergence rate is
determined by the eigenvalue of
$(\mathbf{I}-\eta_{t}\hat{\boldsymbol{\Lambda}}_{t}^{-1}\boldsymbol{\Lambda})$
with the greatest absolute value. First, we show that
$\rho_{t}\notin[\lambda_{n},\lambda_{q_{t}+1}]$ implies $r_{t}>r_{t}^{*}$,
then we show that, if $\rho_{t}\in[\lambda_{n},\lambda_{q_{t}+1}]$, there
exists an optimal $\eta_{t}^{*}$ for which $r_{t}^{*}$ is achieved. If
$\rho_{t}<\lambda_{n}$, the choice of $\eta_{t}$ minimizing the maximum
absolute value of
$(\mathbf{I}-\eta_{t}\hat{\boldsymbol{\Lambda}}_{t}^{-1}\boldsymbol{\Lambda})$
is the solution of $|1-\eta_{t}|=|1-\eta_{t}\lambda_{q_{t+1}}/\rho_{t}|$,
which is $\eta_{t}^{*}=2\rho_{t}/(\rho_{t}+\lambda_{q_{t}+1})$. The
corresponding convergence factor is $1-\eta_{t}^{*}>r_{t}^{*}$. Similarly, if
$\rho_{t}>\lambda_{q_{t}+1}$, one gets
$\eta_{t}^{*}=2\rho_{t}/(\rho_{t}+\lambda_{n})$ and convergence factor equal
to $1-2\lambda_{n}/(\rho_{t}+\lambda_{n})>r_{t}^{*}$. If
$\rho_{t}\in[\lambda_{n},\lambda_{q_{t}+1}]$, the best $\eta_{t}$ is such that
$|1-\eta_{t}\lambda_{n}/\rho_{t}|=|1-\eta_{t}\lambda_{q_{t+1}}/\rho_{t}|$,
whose solution is

$\eta_{t}^{*}=\frac{2\rho_{t}}{\lambda_{q_{t}+1}+\lambda_{n}}$ (11)

and the achieved factor is
$r_{t}=1-\eta_{t}^{*}\lambda_{n}/\rho_{t}=1-\lambda_{n}/\rho_{t}^{*}=r_{t}^{*}$.
We see that the definition of the set $\mathcal{S}^{*}$ immediately follows.

In the above Theorem, we have shown that the best convergence rate is
achievable, by tuning the step size, as long as
$\rho_{t}\in[\lambda_{n},\lambda_{q_{t}+1}]$. In the following Corollary, we
provide an optimal choice for the tuple $(\rho_{t},\eta_{t})$ with respect to
the estimation error.

###### Corollary 2

Among the tuples $(\rho_{t},\eta_{t})\in\mathcal{S}^{*}$, the choice of the
tuple $(\rho_{t}^{*},1)$, with $\rho_{t}^{*}$ defined in (8), is optimal with
respect to the estimation error
$\|\boldsymbol{\theta}^{t+1}_{\rho_{t},\eta_{t}}-\boldsymbol{\theta}^{*}\|$,
for any $\operatorname*{{\boldsymbol{\theta}}}^{t}$ and for any $t$, in the
sense that

$\|\boldsymbol{\theta}^{t+1}_{\rho_{t}^{*},1}-\boldsymbol{\theta}^{*}\|\leq\|\boldsymbol{\theta}^{t+1}_{\rho_{t},\eta_{t}}-\boldsymbol{\theta}^{*}\|,\
\ \forall{(\eta_{t},\rho_{t})}\in\mathcal{S}^{*}$

Proof  Define
$\mathbf{B}_{t}^{*}:=(\mathbf{I}-(\hat{\boldsymbol{\Lambda}}_{t}(\rho_{t}^{*},q_{t}))^{-1}\boldsymbol{\Lambda})=diag(0,...,0,1-\lambda_{q_{t}+1}/\rho_{t}^{*},...,1-\lambda_{n}/\rho_{t}^{*})$.
For $\rho_{t}\neq\rho_{t}^{*}$, with
$\rho_{t}\in[\lambda_{n},\lambda_{q_{t}+1}]$, define
$\mathbf{B}_{t}:=(\mathbf{I}-\eta_{t}^{*}(\hat{\boldsymbol{\Lambda}}_{t}(\rho_{t},q_{t}))^{-1}\boldsymbol{\Lambda})=diag(1-\eta_{t}^{*},...,1-\eta_{t}^{*},1-\lambda_{q_{t}+1}/\rho_{t}^{*},...,1-\lambda_{n}/\rho_{t}^{*})=\mathbf{B}^{*}_{t}+\delta\mathbf{B}_{t}$,
with $\delta\mathbf{B}_{t}=diag(1-\eta_{t}^{*},...,1-\eta_{t}^{*},0,...,0)$,
where $\eta_{t}^{*}$ is defined in (11) and $\rho_{t}^{*}$ in (8). Now define
$\mathbf{z}^{t}:=\mathbf{V}^{T}(\boldsymbol{\theta}^{t}-\boldsymbol{\theta}^{*})$
and
$\mathbf{z}^{t+1}_{\rho_{t},\eta_{t}}:=\mathbf{V}^{T}(\boldsymbol{\theta}^{t+1}_{\rho_{t},\eta_{t}}-\boldsymbol{\theta}^{*})$,
where
$(\boldsymbol{\theta}^{t+1}_{\rho_{t},\eta_{t}}-\boldsymbol{\theta}^{*})$ is
defined in (10). We have

$\displaystyle\|\boldsymbol{\theta}^{t+1}_{\rho_{t}^{*},1}-\boldsymbol{\theta}^{*}\|^{2}=\|\mathbf{z}^{t+1}_{\rho_{t}^{*},1}\|^{2}=\|\mathbf{B}^{*}_{t}\mathbf{z}^{t}\|^{2},$
$\|\boldsymbol{\theta}^{t+1}_{\rho_{t},\eta_{t}^{*}}-\boldsymbol{\theta}^{*}\|^{2}=\|\mathbf{B}_{t}\mathbf{z}^{t}\|^{2}=\|\mathbf{B}^{*}_{t}\mathbf{z}^{t}\|^{2}+\|\delta\mathbf{B}_{t}\mathbf{z}^{t}\|^{2},$

because the cross term is
$2(\mathbf{B}^{*}_{t}\mathbf{z}^{t})^{T}(\delta\mathbf{B}_{t}\mathbf{z}^{t})=0$.
We see that, for any $t$ and for any $\boldsymbol{\theta}^{t}$,

$\|\boldsymbol{\theta}^{t+1}_{\rho_{t}^{*},1}-\boldsymbol{\theta}^{*}\|\leq\|\boldsymbol{\theta}^{t+1}_{\rho_{t},\eta_{t}^{*}}-\boldsymbol{\theta}^{*}\|$

We remark that the bound in (6) is tight for $r_{t}=r_{t}^{*}$. If $q_{t}$
increases, the convergence factor $r_{t}^{*}$ decreases until it becomes zero,
when $q_{t}=n-1$, thus we can have convergence in a finite number of steps.

Figure 2: Performance comparison for different values of $\rho_{t}$. The best
choice in terms of estimation error from Corollary 2 is compared against the
choice that was proposed in Erdogdu and Montanari (2015), that is
$\rho_{t}=\lambda_{q_{t}+1}$. In (b), we show an example of the set
$\mathcal{S}^{*}$, and outline two points: LSs in red is the choice of the
tuple in $\mathcal{S}^{*}$ providing the optimal convergence factor in LSs,
while the purple CVX which does not belong to $\mathcal{S}^{*}$, and it is the
choice which we do in the scenario of FL with convex cost

###  3.2 Federated least squares

We now consider the FL scenario described in section 2, in which $M$ agents
keep their local data and share optimization parameters to contribute to the
learning algorithm.

###### Assumption 2

In the rest of the paper, we assume that each agent has the same amount of
data samples, $N_{i}=\frac{N}{M}$. This allows us to express global functions,
such as the gradient, as the arithmetic mean of local functions (e.g.,
$\mathbf{g}_{t}=(1/M)\sum_{i=1}^{M}\mathbf{g}_{t}^{(i)}$).

This assumption is made only for notation convenience. It is straightforward
to show that all the results are valid also for $N_{i}$ different for each
$i$. To show it, it is sufficient to replace the arithmetic mean of local
functions with the weighted average, weighting each local function with
$p_{i}=N_{i}/N$. As an example, the global gradient would be written in the
following way:

$\mathbf{g}_{t}=\sum_{i=1}^{M}p_{i}\mathbf{g}_{t}^{(i)},p_{i}=N_{i}/N$

In this subsection we introduce the algorithm for the LSs case (Algorithm 1),
that is a special case of Algorithm 4, described in Sec. 5, which is designed
for a general convex cost. We refer to Algorithm 1 as SHED-LS and it works as
follows: at iteration $t$, each device shares with the Master some of its
local Hessian eigenvectors scaled by quantities related to the corresponding
eigenvalues. The eigenvectors are shared _incrementally_ , where the order in
which they are shared is given by the corresponding eigenvalues. For example,
at $t=1$ agent $i$ will start by sharing its first Hessian eigenvectors
$\mathbf{v}_{1}^{(i)},...,\mathbf{v}_{q_{1}^{(i)}}^{(i)}$ according to its
communication resources (CRs) and then will incrementally send up to
$\mathbf{v}_{n-1}^{(i)}$ in the following iterations. To enable the
approximation of the local Hessian via a limited number of eigenvectors using
(4), each eigenvector is sent to the Master together with $\lambda_{j}^{(i)}$
and the parameter $\rho_{t}^{(i)}$.

Input:
$\\{\mathcal{D}_{i}\\}_{i=1}^{M}=\\{(\mathbf{X}_{i},\mathbf{Y}_{i})\\}_{i=1}^{M}$,
$t_{max}$, $\boldsymbol{\theta}^{0}$, $\mathcal{A}=\\{agents\\}$,
$\mathcal{I}=\\{1\\}$

Output: $\operatorname*{{\boldsymbol{\theta}}}^{t_{max}}$

1 for  _$t\leftarrow 1$ to $t_{max}$_ do

2 for  _$agent\ i\in\mathcal{A}$_ do

3 when  _Received $\operatorname*{{\boldsymbol{\theta}}}^{t}$ from the Master_
do

4 if  _$t\in\mathcal{I}$_ then

5 compute
$\mathbf{H}_{LS}^{(i)}=\nabla^{2}f^{(i)}(\boldsymbol{\theta}^{t})=\mathbf{X}_{i}\mathbf{X}_{i}^{T}$

$\\{({\lambda}_{j}^{(i)},{\mathbf{v}}_{j}^{(i)})\\}_{j=1}^{n}\leftarrow
eigendec(\mathbf{H}_{LS}^{(i)})$

// eigendecomposition

6 $q_{0}^{(i)}\leftarrow 0$

7 compute $\mathbf{g}_{t}^{(i)}=\nabla f^{(i)}(\boldsymbol{\theta}^{t})$

set $d_{t}^{(i)}$

// according to CRs

$q_{t}^{(i)}\leftarrow q_{t-1}^{(i)}+d_{t}^{(i)}$

// increment

${\rho}_{t}^{(i)}\leftarrow({\lambda}_{q_{t}^{(i)}+1}^{(i)}+{\lambda}_{n}^{(i)})/2$

// approximation parameter

8
$U_{t}^{(i)}\leftarrow\\{\\{{\mathbf{v}}_{j}^{(i)},\lambda_{j}^{(i)}\\}_{j=q^{(i)}_{t-1}+1}^{q^{(i)}_{t}},\mathbf{g}_{t}^{(i)},{\rho}_{t}^{(i)}\\}$

9 Send $U_{t}^{(i)}$ to the Master.

10

11

12 At the Master:

13 when  _Received $U_{t}^{(i)}$ from all agents_ do

compute $\hat{\mathbf{H}}_{t}^{(i)},\ \forall{i}$.

// see (13)

14 Compute $\hat{\mathbf{H}}_{t}$ (as in eq. (12)) and $\mathbf{g}_{t}$.

15 Perform Newton-type update (3) with $\eta_{t}=1$.

16 Broadcast $\boldsymbol{\theta}^{t+1}$ to all agents.

17

Algorithm 1 Federated Least Squares - SHED-LS

The Master averages the received information to obtain an estimate of the
global Hessian as follows:

${\hat{\mathbf{H}}}_{t}=\frac{1}{M}\sum_{i=1}^{M}{\hat{\mathbf{H}}}_{t}^{(i)},$
(12)

where $\hat{\mathbf{H}}_{t}^{(i)}$ is

$\hat{\mathbf{H}}_{t}^{(i)}:=\hat{\mathbf{H}}^{(i)}(\rho_{t}^{(i)},q_{t}^{(i)})=\sum_{j=1}^{q_{t}^{(i)}}(\lambda_{j}^{(i)}-\rho_{t}^{(i)}){\mathbf{v}}_{j}^{(i)}{{\mathbf{v}}_{j}^{(i)T}}+\rho_{t}^{(i)}\mathbf{I},$
(13)

in which $q_{t}^{(i)}$ is the number of the local Hessian eigenvectors-
eigenvalues pairs that agent $i$ has already sent to the Master at iteration
$t$. We denote by $d_{t}^{(i)}$ the _increment_ , meaning the number of
eigenvectors that agent $i$ can send to the Master at iteration $t$. Given the
results shown in Section 3.1, in this Section we fix the local approximation
parameter to be
$\rho_{t}^{(i)}=\rho_{t}^{(i)*}=({\lambda}_{q_{t}^{(i)}+1}^{(i)}+{\lambda}_{n}^{(i)})/2$
and the step size to be $\eta_{t}=1$. The following results related to the
convergence rate allow the value $q_{t}^{(i)}$ to be different for each agent
$i$, so we define

$\mathbf{q}_{t}=[q_{t}^{(1)},...,q_{t}^{(M)}]^{T},\
{q}_{t}^{(i)}\in\\{0,...,n\\}\ \ \forall{i\in\\{1,...,M\\}}$

By construction, the matrix $\hat{\mathbf{H}}_{t}$ is positive definite, being
the sum of positive definite matrices, implying that
$-\mathbf{p}_{t}=-\hat{\mathbf{H}}_{t}^{-1}\mathbf{g}_{t}$ is a descent
direction.

###### Theorem 3

Consider the problem in (1) in the least squares case. Given
$\hat{\mathbf{H}}_{t}$ defined in (12), the update rule defined in (3) is such
that, for $\rho_{t}^{(i)}=(\lambda_{q_{t}^{(i)}+1}^{(i)}+\lambda_{n}^{(i)})/2$
and $\eta_{t}=1$:

$\|\boldsymbol{\theta}^{t+1}-\boldsymbol{\theta}^{*}\|\leq
c_{t}\|\boldsymbol{\theta}^{t}-\boldsymbol{\theta}^{*}\|,$ (14)

with $c_{t}=(1-\bar{\lambda}_{n}/\bar{\rho}_{t})$ and

$\bar{\rho}_{t}:=\bar{\rho}(\mathbf{q}_{t})=\frac{1}{M}\sum_{i=1}^{M}\rho_{t}^{(i)},\
\ \bar{\lambda}_{n}=\frac{1}{M}\sum_{i=1}^{M}\lambda_{n}^{(i)}.$ (15)

If Algorithm 1 is applied, $c_{t+1}\leq c_{t}\ \forall{t}$, and, if for all
$i$ it holds that $q_{t^{\prime}}^{(i)}={n-1}$, at some iteration $t^{\prime}$
$c_{t^{\prime}}=0$.

Proof  Fix $\eta_{t}=1$ in (5),

$\boldsymbol{\theta}^{t+1}-\boldsymbol{\theta}^{*}=\mathbf{A}_{t}({\operatorname*{{\boldsymbol{\theta}}}}^{t}-{\operatorname*{{\boldsymbol{\theta}}}}^{*})=(\mathbf{I}-\hat{\mathbf{H}}_{t}^{-1}\mathbf{H}_{LS})({\operatorname*{{\boldsymbol{\theta}}}}^{t}-{\operatorname*{{\boldsymbol{\theta}}}}^{*}).$

We have that

$\displaystyle\|{\operatorname*{{\boldsymbol{\theta}}}}^{t+1}-{\operatorname*{{\boldsymbol{\theta}}}}^{*}\|$
$\displaystyle\leq\|\mathbf{I}-\hat{\mathbf{H}}_{t}^{-1}\mathbf{H}_{LS}\|\|{\operatorname*{{\boldsymbol{\theta}}}}^{t}-{\operatorname*{{\boldsymbol{\theta}}}}^{*}\|$
$\displaystyle\leq\|\hat{\mathbf{H}}_{t}^{-1}\|\|\hat{\mathbf{H}}_{t}-\mathbf{H}_{LS}\|\|{\operatorname*{{\boldsymbol{\theta}}}}^{t}-{\operatorname*{{\boldsymbol{\theta}}}}^{*}\|$
$\displaystyle\leq\frac{(\bar{\rho}_{t}-\bar{\lambda}_{n})}{\bar{\rho}_{t}}\|{\operatorname*{{\boldsymbol{\theta}}}}^{t}-{\operatorname*{{\boldsymbol{\theta}}}}^{*}\|$

The last inequality follows from two inequalities: (i)
$\|\hat{\mathbf{H}}_{t}^{-1}\|\leq 1/\bar{\rho}_{t}$ and (ii)
$\|\hat{\mathbf{H}}_{t}-\mathbf{H}_{LS}\|\leq\bar{\rho}_{t}-\bar{\lambda}_{n}$.
(i) holds because
$\|\hat{\mathbf{H}}_{t}^{-1}\|=(\lambda_{min}(\hat{\mathbf{H}}_{t}))^{-1}$,
and $\lambda_{min}(\hat{\mathbf{H}}_{t})\geq\bar{\rho}_{t}$, thus implying
$\|\hat{\mathbf{H}}_{t}^{-1}\|\leq 1/\bar{\rho}_{t}$.
(ii) follows recalling that
$\mathbf{H}_{LS}=\frac{1}{M}\sum_{i=1}^{M}\mathbf{H}^{(i)}_{LS}$, with
$\mathbf{H}^{(i)}_{LS}$ the local Hessian at agent $i$. We have

$\|\hat{\mathbf{H}}_{t}-\mathbf{H}_{LS}\|=\frac{1}{M}\|\sum_{i=1}^{M}(\hat{\mathbf{H}}_{t}^{(i)}-\mathbf{H}^{(i)}_{LS})\|\leq\frac{1}{M}\sum_{i=1}^{M}\|\hat{\mathbf{H}}_{t}^{(i)}-\mathbf{H}^{(i)}_{LS}\|.$

Being $\hat{\mathbf{H}}_{t}^{(i)}-\mathbf{H}^{(i)}_{LS}$ symmetric, it holds
that

$\|\hat{\mathbf{H}}_{t}^{(i)}-\mathbf{H}^{(i)}_{LS}\|=\max_{j}|\lambda_{j}(\hat{\mathbf{H}}_{t}^{(i)}-\mathbf{H}^{(i)}_{LS})|=\rho_{t}^{(i)}-\lambda_{n}^{(i)},$

where the last equality holds because

$\hat{\mathbf{H}}_{t}^{(i)}-\mathbf{H}^{(i)}_{LS}=\mathbf{V}^{(i)}(\hat{\boldsymbol{\Lambda}}_{t}^{(i)}-\boldsymbol{\Lambda}_{t}^{(i)})\mathbf{V}^{(i)T}$
(16)

where

$\displaystyle\hat{\boldsymbol{\Lambda}}_{t}^{(i)}=diag(\lambda_{1}^{(i)},...,\lambda_{q_{t}^{(i)}}^{(i)},\rho_{t}^{(i)},...,\rho_{t}^{(i)}),$
(17)
$\displaystyle\boldsymbol{\Lambda}_{t}^{(i)}=diag(\lambda_{1}^{(i)},...,\lambda_{q_{t}^{(i)}}^{(i)},\lambda_{q_{t}^{(i)}+1}^{(i)},...\lambda_{n}^{(i)}),$

and because
$\rho_{t}^{(i)}=(\lambda_{q_{t}^{(i)}+1}^{(i)}+\lambda_{n}^{(i)})/2$.

This theorem shows that our approach in the LSs case provides convergence in a
finite number of iterations, if $q_{t}^{(i)}$ keeps increasing through time
for each agent $i$. Indeed, as in the centralized case, if $q_{t}^{(i)}$
increases for all $i$, the factor $c_{t}$ decreases until it becomes zero.
Furthermore, each agent is free to send at each iterations an arbitrary number
of eigenvector-eigenvalue pairs, according to its CRs, and by doing so it can
improve the convergence rate.

##  4 From least squares to strongly convex cost

We want to extend the analysis and algorithm presented in the previous
sections of the paper to a general convex cost
$f(\operatorname*{{\boldsymbol{\theta}}}^{t})$. With respect to the proposed
approach, the general convex case requires special attention for two main
reasons: (i) the update rule defined in (3) requires tuning of the step-size
$\eta_{t}$, usually via backtracking line search, and (ii) the Hessian matrix
is in general a function of the parameter
$\operatorname*{{\boldsymbol{\theta}}}$. In this section, still focusing on
the least squares case, we provide some results that are instrumental to the
analysis of the general convex case.

###  4.1 Backtracking line search for step size tuning

We recall the well-known Armijo-Goldstein condition for accepting a step size
$\eta_{t}$ via backtracking line search:

$f(\boldsymbol{\theta}^{t}-\eta_{t}\mathbf{p})\leq
f(\boldsymbol{\theta}^{t})-\alpha\eta_{t}\mathbf{p}^{T}\mathbf{g}_{t},$ (18)

where $\alpha\in(0,1/2)$. The corresponding line search algorithm is the
following:

Input: $\alpha\in(0,1/2)$, $\beta\in(0,1)$,
$\operatorname*{{\boldsymbol{\theta}}}^{t}$, $\mathbf{p}_{t}$,
$\mathbf{g}_{t}$, $f$

Output: $\bar{\eta}_{t}$

1 $\eta_{t}^{0}\leftarrow 1$, $k\leftarrow 0$

2 while  _$f(\boldsymbol{\theta}^{t}-\eta_{t}^{(k)}\mathbf{p})
>f(\boldsymbol{\theta}^{t})-\alpha\eta_{t}^{(k)}\mathbf{p}^{T}\mathbf{g}_{t}$_
do

3 $k\leftarrow k+1$

4 $\eta_{t}^{(k)}=\beta\eta_{t}^{(k-1)}$

$\bar{\eta}_{t}=\eta_{t}^{(k)}$

Algorithm 2 Backtracking line search algorithm

###### Lemma 4

Consider the problem in (1) in the least squares case. Let
$\mathbf{p}_{t}=\hat{\mathbf{H}}_{t}^{-1}\mathbf{g}_{t}$, with
$\hat{\mathbf{H}}_{t}$ defined in (12). A sufficient condition for a step size
$\eta_{t}$ to satisfy Armijo-Goldstein condition (18), for any
$\alpha\in(0,1/2)$, is

$\eta_{t}=\min_{i=1,...,M}{\frac{\rho_{t}^{(i)}}{\lambda_{q^{(i)}_{t}+1}^{(i)}}}$
(19)

Proof  The quadratic cost in $\boldsymbol{\theta}^{t}$ can be written as

$f(\boldsymbol{\theta}^{t})=f(\boldsymbol{\theta}^{*})+\bar{f}(\boldsymbol{\theta}^{t}).$

with
$\bar{f}(\boldsymbol{\theta}^{t})=\frac{1}{2}(\boldsymbol{\theta}^{t}-\boldsymbol{\theta}^{*})^{T}\mathbf{H}_{LS}(\boldsymbol{\theta}^{t}-\boldsymbol{\theta}^{*})$.
Given that $f(\boldsymbol{\theta}^{*})$ does not depend on
$\boldsymbol{\theta}^{t}$, we can focus on $\bar{f}(\boldsymbol{\theta}^{t})$.
We have that

$\displaystyle\bar{f}(\boldsymbol{\theta}^{t}-\eta_{t}\mathbf{p}_{t})$
$\displaystyle=\frac{1}{2}(\boldsymbol{\theta}^{t}-\eta_{t}\mathbf{p}_{t}-\boldsymbol{\theta}^{*})^{T}\mathbf{H}_{LS}(\boldsymbol{\theta}^{t}-\eta_{t}\mathbf{p}_{t}-\boldsymbol{\theta}^{*})$
(20)
$\displaystyle\stackrel{{\scriptstyle(1)}}{{=}}\bar{f}(\boldsymbol{\theta}^{t})+\frac{1}{2}\eta_{t}^{2}\mathbf{p}_{t}^{T}\mathbf{H}_{LS}\mathbf{p}_{t}-\eta_{t}\mathbf{p}_{t}^{T}\mathbf{g}_{t}$
$\displaystyle\stackrel{{\scriptstyle(2)}}{{=}}\bar{f}(\boldsymbol{\theta}^{t})-\eta_{t}\mathbf{p}_{t}^{T}(\hat{\mathbf{H}}_{t}-\eta_{t}\frac{\mathbf{H}_{LS}}{2})\mathbf{p}_{t}$

where we have used identity
$\mathbf{H}_{LS}(\boldsymbol{\theta}_{t}-\boldsymbol{\theta}^{*})=\mathbf{g}_{t}$
and the fact that
$\mathbf{p}_{t}^{T}\mathbf{g}_{t}=\mathbf{p}_{t}^{T}\hat{\mathbf{H}}_{t}\hat{\mathbf{H}}_{t}^{-1}\mathbf{g}_{t}=\mathbf{p}_{t}^{T}\hat{\mathbf{H}}_{t}\mathbf{p}_{t}$
to get equality (1) and (2), respectively. We see that if

$\hat{\mathbf{H}}_{t}-\eta_{t}\mathbf{H}_{LS}/2\geq\hat{\mathbf{H}}_{t}/2,$
(21)

then Armijo-Goldstein condition (18) is satisfied. Indeed, in that case,

$\displaystyle\bar{f}(\boldsymbol{\theta}^{t})-\eta_{t}\mathbf{p}_{t}^{T}(\hat{\mathbf{H}}_{t}-\eta_{t}\frac{\mathbf{H}_{LS}}{2})\mathbf{p}_{t}$
$\displaystyle\leq\bar{f}(\boldsymbol{\theta}^{t})-\frac{1}{2}\eta_{t}\mathbf{p}_{t}^{T}\hat{\mathbf{H}}_{t}\mathbf{p}_{t}$
$\displaystyle\leq\bar{f}(\boldsymbol{\theta}^{t})-\alpha\eta_{t}\mathbf{p}_{t}^{T}\mathbf{g}_{t}.$

So, we need to find a sufficient condition on $\eta_{t}$ for (21) to be true.
We see that (21) is equivalent to
$\hat{\mathbf{H}}_{t}-\eta_{t}\mathbf{H}_{LS}\geq 0$. We have
$\hat{\mathbf{H}}_{t}-\eta_{t}\mathbf{H}_{LS}=\frac{1}{M}\sum_{i=1}^{M}\mathbf{V}^{(i)}(\hat{\boldsymbol{\Lambda}}_{t}^{(i)}-\eta_{t}\boldsymbol{\Lambda}_{t}^{(i)})\mathbf{V}^{(i)T}$,
and, $\forall{i}$, all the elements of the diagonal matrix
$\hat{\boldsymbol{\Lambda}}_{t}^{(i)}-\eta_{t}\boldsymbol{\Lambda}_{t}^{(i)}$
are positive if $\eta_{t}\leq\rho^{(i)}_{t}/\lambda_{q_{t}^{(i)}+1}^{(i)}$
(see Eq. (17)), and we see that the choice (19) provides a sufficient
condition to satisfy (21), from which we can conclude.

###### Corollary 5

In the least squares case (Algorithm 1), choosing
$\rho_{t}^{(i)}=\rho_{t}^{(i)*},\ \forall{i}$, Armijo backtracking line search
(Algorithm 2) would choose a step size $\eta_{t}\geq\frac{1}{2}$. Choosing
$\rho_{t}^{(i)}=\lambda_{q_{t}^{(i)}+1},\ \forall{i}$, Armijo backtracking
line search would choose a step size $\eta_{t}=1$.

Proof  The proof is straightforward from Eqn. (19) of Lemma 4.

###### Remark 6

For the choice $\rho^{(i)}=\rho^{(i)*}$, the Armijo-backtracking might not
choose a step size $\eta_{t}=1$ even for arbitrarily small $\alpha<0.5$.
Indeed, we can easily build a counter-example in the centralized case
considering
$\hat{\mathbf{H}}_{t}=\hat{\boldsymbol{\Lambda}}_{t}=diag(\lambda_{1},...,\lambda_{q},\rho^{*},...,\rho^{*})$,
a gradient $\mathbf{g}_{t}$ such that
$\mathbf{p}_{t}={\mathbf{\hat{H}}}_{t}^{-1}\mathbf{g}_{t}=[0,...,0,1,0,...,0]^{\top}$
(e.g., $\mathbf{g}_{t}=[0,...0,\lambda_{q_{t}+1},0,...,0]^{\top}$). We see
that with $\eta_{t}=1$, Eq. (20) becomes
$f({\operatorname*{{\boldsymbol{\theta}}}}^{t+1})=f({\operatorname*{{\boldsymbol{\theta}}}}^{t})-\lambda_{n}/2$
and, in order to be satisfied, the Armijo condition would require
$\alpha\leq\frac{\lambda_{n}}{\lambda_{q+1}+\lambda_{n}}$, where the right
hand side can become arbitrarily small depending on the eigenspectrum.

The results of Lemma 4 and of Corollary 5 and the counter-example of the above
Remark are important for the design of the algorithm in the general convex
case. Indeed, as illustrated in the next Section (Section 5), a requirement
for the theoretical results on the convergence rate is that the step size
becomes equal to one, which is not guaranteed by the Armijo backtracking line
search, even when considering the least squares case, if
$\rho_{t}^{(i)}<\lambda_{q_{t}^{(i)}+1}^{(i)}$. For this reason, the algorithm
in the general convex case is designed with
$\rho_{t}^{(i)}=\lambda_{q_{t}^{(i)}+1}^{(i)}$

###  4.2 Algorithm with periodic renewals

Now, we introduce a variant of Algorithm 1 that is instrumental to study the
convergence rate of the proposed algorithm in the general convex case. The
variant is Algorithm 3. The definition of $\mathcal{I}=\\{1,T,2T,...\\}$
implies that every $T$ iterations the incremental strategy is restarted from
the first EEPs of $\mathbf{H}_{LS}$, in what we call a periodic _renewal_.

In Algorithm 1, substitute $\mathcal{I}=\\{1\\}$ with
$\mathcal{I}=\\{1,T,2T,...\\}$, for some input parameter $T<n$.

Algorithm 3 Variant of Algorithm 1, SHED-LS-periodic

Differently from Algorithm 1, Algorithm 3 can not guarantee convergence in a
finite number of steps, because $T<n$ and thus it could be that $c_{t}>0,\
\forall{t}$ (see Theorem 3). We study the convergence rate of the algorithm by
focusing on upper bounds on the Lyapunov exponent (Lyapunov, 1992) of the
discrete-time dynamical system ruled by the descent algorithm. The Lyapunov
exponent characterizes the rate of exponential (linear) convergence and it is
defined as the positive constant $a_{*}>0$ such that, considering
$h(\operatorname*{{\boldsymbol{\theta}}}^{t}):=(\operatorname*{{\boldsymbol{\theta}}}^{t}-\operatorname*{{\boldsymbol{\theta}}}^{*})a^{-t}$,
if $a>a_{*}$ then $h(\operatorname*{{\boldsymbol{\theta}}}^{t})$ vanishes with
$t$, while if $a<a_{*}$, for some initial condition,
$h(\operatorname*{{\boldsymbol{\theta}}}^{t})$ diverges. The usual definition
of Lyapunov exponent for discrete-time linear systems (see Czornik et al.,
2012) is, considering the system defined in (5),

$a_{*}:=\limsup_{t\rightarrow\infty}{\|\boldsymbol{\Psi}_{t}\|^{1/t}},\
\boldsymbol{\Psi}_{t}=\mathbf{A}_{1}\cdots\mathbf{A}_{t}.$ (22)

From (14), we have that, for each $k$, $\|\mathbf{A}_{k}\|\leq
c_{k}=(1-\bar{\lambda}_{n}/\bar{\rho}_{k})$. This implies that, defining

$\displaystyle a_{t}:$ $\displaystyle=(\prod_{k=1}^{t}{c_{k}})^{1/t},$ (23)
$\displaystyle\bar{a}:$ $\displaystyle=\limsup_{t\rightarrow+\infty}{a_{t}},$

it is $a_{*}\leq\bar{a}$. The following Lemma formalizes this bound and
provides an upper bound on the Lyapunov exponent obtained by applying
Algorithm 3.

###### Lemma 7

Let $c_{k}=(1-\bar{\lambda}_{n}/\bar{\rho}_{k})$ (see Eq. (15)). Applying
Algorithm 3, the Lyapunov exponent of system (5) is such that

$a_{*}\leq\bar{a}=\limsup_{t\rightarrow+\infty}(\prod_{k=1}^{t}{c_{k}})^{1/t}.$

If $q_{t}^{(i)}=q_{t-1}^{(i)}+1,\ \forall{i,t},$

$a_{*}\leq\bar{a}_{T}:=(\prod_{k=1}^{T}c_{k})^{1/T}.$ (24)

where $\bar{a}_{T}$ is such that $\bar{a}_{T+1}\leq\bar{a}_{T}$ and
$\bar{a}_{n}=0$.

Proof  We see that $a_{*}\leq\bar{a}$ from definition (22), because for each
$k$, $\|\mathbf{A}_{k}\|\leq c_{k}=(1-\bar{\lambda}_{n}/\bar{\rho}_{k})$. Now,
we show that, if $q_{t}^{(i)}=q_{t-1}^{(i)}+1,\ \forall{i,t}$, then
$\bar{a}=\bar{a}_{T}$. To show this latter inequality, let us consider the
logarithm of the considered values, so we prove
$\log{\bar{a}}=\log{\bar{a}_{T}}$. Indeed, if $q_{t}^{(i)}=q_{t-1}^{(i)}+1$
and applying Algorithm 3, $c_{k}$ is a periodic sequence of the index $k$
($c_{k+T}=c_{k}$):

$\displaystyle\log{\bar{a}}$
$\displaystyle=\limsup_{t}{\frac{1}{t}\sum_{k=1}^{t}\log{c_{k}}}$ (25)
$\displaystyle=\limsup_{R}{\frac{1}{RT+T^{\prime}}(R\sum_{k=1}^{T}\log{c_{k}}+\sum_{k=1}^{T^{\prime}}\log{c_{k}}})$
$\displaystyle=\frac{1}{T}\sum_{k=1}^{T}\log{c_{k}}=\log{\bar{a}_{T}}$

where $R=\lfloor t/T\rfloor$ and $T^{\prime}=t-RT$.

##  5 Federated learning with convex cost

Given the previous analysis and theoretical results for linear regression with
quadratic cost, we are now ready to illustrate our Newton-type algorithm
(Algorithm 4) for general convex FL problems, of which Algorithm 1 is a
special case. We refer to the this general version of the algorithm simply as
SHED.

Input: $\\{\mathcal{D}_{i}\\}_{i=1}^{M}$, $t_{max}$, $\mathcal{I}$,
$\boldsymbol{\theta}^{0}$, $\mathcal{A}=\\{agents\\}$

Output: $\operatorname*{{\boldsymbol{\theta}}}^{t_{max}}$

1 for  _$t\leftarrow 0$ to $t_{max}$_ do

2 for  _$agent\ i\in\mathcal{A}$_ do

3 when  _Received $\operatorname*{{\boldsymbol{\theta}}}^{t}$ from the Master_
do

4 if  _$t\in\mathcal{I}$_ then

5 $k_{t}\leftarrow t$

compute $\mathbf{H}_{t}^{(i)}=\nabla^{2}f^{(i)}(\boldsymbol{\theta}^{t})$

// renewal

$\\{(\hat{\lambda}_{j,t}^{(i)},\hat{\mathbf{v}}_{j,t}^{(i)})\\}_{j=1}^{n}\leftarrow
eigendec(\mathbf{H}_{t}^{(i)})$

// eigendecomposition

6 $q_{t-1}^{(i)}\leftarrow 0$

7 compute
$\mathbf{g}_{t}^{(i)}=\mathbf{g}^{(i)}(\boldsymbol{\theta}^{t})=\nabla
f^{(i)}(\boldsymbol{\theta}^{t})$

set $d_{t}^{(i)}$

// according to CRs

$q_{t}^{(i)}\leftarrow q_{t-1}^{(i)}+d_{t}^{(i)}$

// increment

$\hat{\rho}_{t}^{(i)}\leftarrow\hat{\lambda}_{q_{t}^{(i)}+1,t}^{(i)}$

// approximation parameter

$U_{t}^{(i)}\leftarrow\\{\\{{\hat{\mathbf{v}}}_{j,t}^{(i)},\hat{\lambda}_{j,t}^{(i)}\\}_{j=q^{(i)}_{t-1}+1}^{q^{(i)}_{t}},\mathbf{g}_{t}^{(i)},{\hat{\rho}}_{t}^{(i)}\\}$

// see (28) and (30)

8 Send $U_{t}^{(i)}$ to the Master.

9

10

11 At the Master:

12 when  _Received $U_{t}^{(i)}$ from all agents_ do

compute $\hat{\mathbf{H}}_{t}^{(i)},\ \forall{i}$.

// see (30)

13 Compute $\hat{\mathbf{H}}_{t}$ (as in eq. (26)) and $\mathbf{g}_{t}$.

14 Get $\eta_{t}$ via federated Armijo backtracking line search.

15 Perform Newton-type update (3).

16 Broadcast $\boldsymbol{\theta}^{t+1}$ to all agents.

17

Algorithm 4 FL with convex cost - SHED

Since in a general convex problem the Hessian depends on the current
parameter, $\boldsymbol{\theta}^{t}$, we denote by
$\mathbf{H}(\boldsymbol{\theta}^{t})$ the global Hessian at the current
iterate, while we denote by $\hat{\mathbf{H}}_{t}$ the global approximation,
defined similarly to (12), with the difference that now eigenvalues and
eigenvectors depend on the parameter for which the Hessian was computed.
The expression of $\hat{\mathbf{H}}_{t}$ thus becomes:

$\hat{\mathbf{H}}_{t}=\frac{1}{M}\sum_{i=1}^{M}\hat{\mathbf{H}}_{t}^{(i)}{(\boldsymbol{\theta}^{k_{t}^{(i)}})},$
(26)

where $k_{t}^{(i)}\leq t$ denotes the iteration in which the local Hessian of
agent $i$ was computed. The parameter $\boldsymbol{\theta}^{k_{t}^{(i)}}$ is
the parameter for which agent $i$ computed the local Hessian, that in turn is
being used for the update at iteration $t$.
The idea of the algorithm is to use previous versions of the Hessian rather
than always recomputing it. This is motivated by the fact that as we approach
the solution of the optimization problem, the second order approximation
becomes more accurate and the Hessian changes more slowly. Hence, recomputing
the Hessian and restarting the incremental approach provides less and less
advantages as we proceed. From time to time, however, we need to re-compute
the Hessian corresponding to the current parameter $\boldsymbol{\theta}^{t}$,
because $\mathbf{H}(\boldsymbol{\theta}^{k_{t}})$ could have become too
different from $\mathbf{H}(\boldsymbol{\theta}^{t})$. As in Section 4.2, we
call this operation a _renewal_. We denote by $\mathcal{I}$ the set of
iteration indices at which a renewal takes place. In principle, each agent
could have its own set of renewal indices, and decide to recompute the Hessian
matrix independently. In this work, we consider for simplicity that the set
$\mathcal{I}$ is the same for all agents, meaning that all agents use the same
parameter for the local Hessian computation, i.e., $k_{t}^{(i)}=k_{t},\
\forall{i}$. At the end of this section we describe heuristic strategies to
choose $\mathcal{I}$ with respect to the theoretical analysis. We remark that
in the case of a quadratic cost, in which the Hessian is constant, one chooses
$\mathcal{I}=\\{1\\}$, and so Algorithm 1 is a special case of Algorithm 4.
The eigendecomposition can be applied to the local Hessian as before, we
define

$\hat{\mathbf{v}}_{j,t}^{(i)}=\mathbf{v}_{j}^{(i)}(\boldsymbol{\theta}^{k_{t}}),\
\ \hat{\lambda}_{j,t}^{(i)}=\lambda_{j}^{(i)}(\boldsymbol{\theta}^{k_{t}}).$
(27)

For notation convenience we also define

$\tilde{\mathbf{v}}_{j,t}^{(i)}=(\hat{\lambda}_{j,t}^{(i)}-\hat{\rho}_{t}^{(i)})^{1/2}\hat{\mathbf{v}}_{j,t}^{(i)},$
(28)

The theoretical results on the convergence rate in this Section require that
$\hat{\mathbf{H}}^{(i)}({\operatorname*{{\boldsymbol{\theta}}}})\geq\mathbf{H}^{(i)}(\operatorname*{{\boldsymbol{\theta}}}),\forall{\operatorname*{{\boldsymbol{\theta}}}}$,
which in turn requires, defining
$\hat{\rho}_{t}^{(i)}={\rho}_{t}^{(i)}(\operatorname*{{\boldsymbol{\theta}}}^{k_{t}})$,
that $\hat{\rho}_{t}^{(i)}\geq\hat{\lambda}^{(i)}_{q_{t}^{(i)}+1,t}$. Given
the results of Theorem 1 on the range of the approximation parameter
$\rho_{t}$ achieving the best convergence rate in the centralized least
squares case, we set:

$\hat{\rho}_{t}^{(i)}=\hat{\lambda}_{q_{t}+1,t}^{(i)}.$ (29)

The local Hessian can be approximated as

$\hat{\mathbf{H}}_{t}^{(i)}(\boldsymbol{\theta}^{k_{t}})=\sum_{j=1}^{q_{t}^{(i)}}{\tilde{\mathbf{v}}_{j,t}^{(i)}\tilde{\mathbf{v}}_{j,t}^{(i)T}}+\hat{\rho}_{t}^{(i)}\mathbf{I}.$
(30)

Clearly, it still holds that
$\hat{\mathbf{H}}_{t}\geq\bar{\rho}_{t}\mathbf{I},\ \forall{t}$, where
$\bar{\rho}_{t}=\frac{1}{M}\sum_{i=1}^{M}\hat{\rho}_{t}^{(i)}$. Furthermore,
it is easy to see that $\hat{\mathbf{H}}_{t}\leq K\mathbf{I}$, with $K$ the
smoothness constant of $f$. We assume to use Armijo backtracking condition,
that is recalled in Algorithm 2.

###### Theorem 8

For any initial condition, Algorithm 2 ensures convergence to the optimum,
i.e.,

$\lim_{t\rightarrow+\infty}\|\boldsymbol{\theta}_{t}-\boldsymbol{\theta}^{*}\|=0.$

Proof  See Appendix A.

Now, we provide results related to the convergence rate of the algorithm. In
order to prove the following results, we need to impose some constraints on
the renewal indices set $\mathcal{I}$. Specifically,

###### Assumption 3

Denoting $\mathcal{I}=\\{C_{j}\\}_{j\in\mathbb{N}}$, there exists a finite
positive integer $\bar{l}$ such that $C_{j}\leq C_{j-1}+\bar{l},\ \forall{j}$.
This is equivalent to state that, writing $k_{t}=t-\tau$, the ‘delay’ $\tau$
is bounded.

The next Theorem provides a bound describing the relation between the
convergence rate and the increments of outdated Hessians.

###### Theorem 9

Applying SHED (Algorithm 4), for any iteration $t$, it holds that:

$\|\boldsymbol{\theta}^{t+1}-\boldsymbol{\theta}^{*}\|\leq
c_{1,t}\|\boldsymbol{\theta}^{t}-\boldsymbol{\theta}^{*}\|+c_{2,t}\|\boldsymbol{\theta}^{t}-\boldsymbol{\theta}^{*}\|^{2}$
(31)

where, defining
$\bar{\lambda}_{n,t}=\frac{1}{M}\sum_{i=1}^{M}\hat{\lambda}_{{n},t}^{(i)}$ and
$\bar{\rho}_{t}=\frac{1}{M}\sum_{i=1}^{M}\hat{\rho}_{t}^{(i)}$ (see (27) and
(29)),

$\displaystyle c_{1,t}$
$\displaystyle=(1-\frac{\bar{\lambda}_{n,t}}{\bar{\rho}_{t}})+\frac{L}{\bar{\rho}_{t}}\|\boldsymbol{\theta}^{t}-\boldsymbol{\theta}^{k_{t}}\|+(1-\eta_{t})\frac{\|{\mathbf{H}}(\boldsymbol{\theta}^{t})\|}{\bar{\rho}_{t}},\
\ $ (32) $\displaystyle c_{2,t}$
$\displaystyle=\frac{\eta_{t}L}{2\bar{\rho}_{t}}.$

Proof  The beginning of the proof follows from the proof of Lemma 3.1 in
Erdogdu and Montanari (2015), (see page 18). In particular, we can get the
same inequality as (A.1) in Erdogdu and Montanari (2015) (the $\mathbf{Q}^{t}$
here is $\hat{\mathbf{H}}_{t}$) with the difference that since we are not sub-
sampling, we have (using the notation of Erdogdu and Montanari (2015))
$S=[n]$. The following inequality holds:

$\displaystyle\|\boldsymbol{\theta}^{t+1}-\boldsymbol{\theta}^{*}\|$
$\displaystyle\leq\|\boldsymbol{\theta}^{t}-\boldsymbol{\theta}^{*}\|\|I-\eta_{t}\hat{\mathbf{H}}_{t}^{-1}\mathbf{H}(\boldsymbol{\theta}^{t})\|$
$\displaystyle+\eta_{t}L\frac{\|\hat{\mathbf{H}}_{t}^{-1}\|}{2}\|\boldsymbol{\theta}^{t}-\boldsymbol{\theta}^{*}\|^{2}.$

Note that, as we have shown in the proof of Theorem 3, it holds
$\|\hat{\mathbf{H}}_{t}^{-1}\|\leq 1/\bar{\rho}_{t}$. We now focus on the
first part of the right hand side of the inequality:

$\displaystyle\|I-\eta_{t}\hat{\mathbf{H}}_{t}^{-1}\mathbf{H}(\boldsymbol{\theta}^{t})\|$
$\displaystyle\leq\|\hat{\mathbf{H}}_{t}^{-1}\|(\|\hat{\mathbf{H}}_{t}-\mathbf{H}(\boldsymbol{\theta}^{t})\|$
$\displaystyle+(1-\eta_{t}){\|{\mathbf{H}}(\boldsymbol{\theta}^{t})\|})$
$\displaystyle\leq\frac{1}{\bar{\rho}_{t}}\|\hat{\mathbf{H}}_{t}-\mathbf{H}(\boldsymbol{\theta}^{t})\|+\frac{(1-\eta_{t})}{\bar{\rho}_{t}}{\|{\mathbf{H}}(\boldsymbol{\theta}^{t})\|}$

and focusing now on the first term of the right hand side of the last
inequality

$\displaystyle\frac{1}{\bar{\rho}_{t}}(\|\hat{\mathbf{H}}_{t}-\mathbf{H}(\boldsymbol{\theta}^{k_{t}})\|+\|\mathbf{H}(\boldsymbol{\theta}^{t})-\mathbf{H}(\boldsymbol{\theta}^{k_{t}})\|)$
$\displaystyle\leq$
$\displaystyle\frac{1}{\bar{\rho}_{t}}(\frac{1}{M}\sum_{i=1}^{M}\|\hat{\mathbf{H}}_{t}^{(i)}(\boldsymbol{\theta}^{k_{t}})-\mathbf{H}^{(i)}(\boldsymbol{\theta}^{k_{t}})\|+L\|\boldsymbol{\theta}^{k_{t}}-\boldsymbol{\theta}^{t}\|)$
$\displaystyle=$ $\displaystyle
1-\frac{\bar{\lambda}_{n,t}}{\bar{\rho}_{t}}+\frac{L}{\bar{\rho}_{t}}\|\boldsymbol{\theta}^{k_{t}}-\boldsymbol{\theta}^{t}\|,$

where the last equality holds being
$\|\hat{\mathbf{H}}_{t}^{(i)}(\boldsymbol{\theta}^{k_{t}})-\mathbf{H}^{(i)}(\boldsymbol{\theta}^{k_{t}})\|=\hat{\rho}_{t}^{(i)}-\lambda_{n,t}^{(i)}$,
which in turn is true given that
$\hat{\rho}_{t}^{(i)}=\hat{\lambda}_{q_{t}^{(i)}+1,t}^{(i)}$.

The above theorem is a generalization of Lemma 3.1 in Erdogdu and Montanari
(2015) (without sub-sampling). In particular, the difference is that (i) the
dataset is distributed and (ii) an outdated Hessian matrix is used.

###### Theorem 10

Recall the definition of the average strong convexity constant
$\bar{\kappa}=(1/M)\sum_{i=1}^{M}\kappa_{i}$, with $\kappa_{i}$ the strong
convexity constant of agent $i$. Let $K$ be the smoothness constant of $f$.
The following results hold:

  1. 1.

If

$3\bar{\kappa}(M(t)+\|{\operatorname*{{\boldsymbol{\theta}}}}^{t}-{\operatorname*{{\boldsymbol{\theta}}}}^{*}\|)+K\|{\operatorname*{{\boldsymbol{\theta}}}}^{t}-{\operatorname*{{\boldsymbol{\theta}}}}^{*}\|\leq\frac{3\bar{\kappa}^{2}}{L}(1-2\alpha),$
(33)

and

$\frac{3}{2}L\|{\operatorname*{{\boldsymbol{\theta}}}}^{t}-{\operatorname*{{\boldsymbol{\theta}}}}^{*}\|+LM(t)\leq\bar{\kappa}$
(34)

with
$M(t)=\max\\{\|{\operatorname*{{\boldsymbol{\theta}}}}^{t}-{\operatorname*{{\boldsymbol{\theta}}}}^{*}\|,\|{\operatorname*{{\boldsymbol{\theta}}}}^{k_{t}}-{\operatorname*{{\boldsymbol{\theta}}}}^{*}\|\\}$
then SHED (Algorithm 4) enjoys at least linear convergence.

  2. 2.

Define $\mathcal{X}_{t}:=\\{k\leq t:\bar{\rho}_{k}=\bar{\lambda}_{n,k}\\}$.
Let $|\mathcal{X}_{t}|$ denote the cardinality of $\mathcal{X}_{t}$. If
$|\mathcal{X}_{t}|=0$ then SHED enjoys linear convergence and the Lyapunov
exponent of the estimation error can be upper bounded as

${a_{*}}\leq\limsup_{t}(\prod_{k=1}^{t}{1-\frac{\bar{\lambda}_{n}^{o}}{\bar{\rho}_{k}^{o}}})^{1/t}$
(35)

where $\bar{\lambda}_{n}^{o}$ and $\bar{\rho}_{k}^{o}$ are the average of the
$n$-th eigenvalues and approximation parameters, respectively, computed at the
optimum:

$\bar{\lambda}_{n}^{o}=\frac{1}{M}\sum_{i=1}^{M}\lambda_{n}^{(i)}(\boldsymbol{\theta}^{*}),\
\ \bar{\rho}_{k}^{o}=\frac{1}{M}\sum_{i=1}^{M}\rho_{k}^{(i)o},$

with $\rho_{k}^{(i)o}=\lambda_{q_{k}^{(i)}+1}^{(i)}(\boldsymbol{\theta}^{*})$.

  3. 3.

Let $|\mathcal{X}_{t}|$ denote the cardinality of $\mathcal{X}_{t}$. If
$|\mathcal{X}_{t}|\geq t^{1/2}h(t)$, with $h(t)$ any function such that
$h(t)\rightarrow\infty$ as $t\rightarrow\infty$, then the Lyapunov exponent is
$a_{*}=0$ and thus SHED enjoys super-linear convergence.

Proof [Sketch of proof] For 1), we first show that when condition (33) holds,
the step size is chosen equal to one by the Armijo backtracking line search.
We then show that when also (34) holds, then the cost converges at least
linearly for any subsequent iteration. For 2), we upper bound the Lyapunov
exponent and exploit local Lipschitz continuity to provide the result. For 3),
we exploit at least linear convergence proved in 1) together with the
assumption on the cardinality of the set $\mathcal{X}_{t}$. For the complete
proof see Appendix A.

In the above theorem we have shown that SHED enjoys at least linear
convergence and provided a sufficient condition on the choice of renewals
indices set to guarantee super-linear convergence. The sufficient condition
could be easily guaranteed for example with a choice of periodic renewals with
a period such that the cardinality of $\mathcal{X}_{t}$ is big enough. Note
that in the case in which all agents send one EEP at each iteration we can get
an explicit expression for the Lyapunov exponent also in the convex cost case,
getting

$a_{*}\leq\bar{a}_{T}:=\prod_{k=1}^{T}(1-\frac{\bar{\lambda}_{n}^{o}}{{\bar{\rho}_{k}^{o}}})^{1/T}.$
(36)

###  5.1 Heuristics for the choice of $\mathcal{I}$

From the theoretical results provided above we can do a heuristic design for
the renewal indices set $\mathcal{I}$. In particular, we see from the bound
(31) in Theorem 9 that when we are at the first iterations of the optimization
we would frequently do the renewal operation, given that the Hessian matrix
changes much faster and that the term
$\|{\operatorname*{{\boldsymbol{\theta}}}}^{k_{t}}-{\operatorname*{{\boldsymbol{\theta}}}}^{t}\|$
is big. As we converge, instead, we would like to reduce the number of renewal
operations, to improve the convergence rate, and this is strongly suggested by
the result 2) related to the Lyapunov exponent in Theorem 10. Furthermore, the
super-linear convergence that follows from 3) in the same Theorem suggests to
keep performing renewals in order to let the cardinality of $\mathcal{X}_{t}$
to grow sufficiently fast with $t$.
To evaluate SHED performance, we obtain the results of the next section using
two different renewal strategies. In the first we choose the distance between
renewals to be determined by the _Fibonacci_ sequence, so
$\mathcal{I}=\\{C_{j}\\}$, where $C_{j}=\sum_{k=1}^{j}F_{k}$, $F_{k}$ being
the _Fibonacci_ sequence, with $F_{0}=0,F_{1}=1$. When the sequence $C_{j}$
reaches $n-1$, the next values of the sequence are chosen so that
$C_{j+1}=C_{j}+n-1$. We will call this method Fib-SHED. The second strategy is
based on the inspection of the value of the gradient norm,
$\|\mathbf{g}_{t}\|$, which is directly related to
$\|{\operatorname*{{\boldsymbol{\theta}}}}^{t}-{\operatorname*{{\boldsymbol{\theta}}}}^{*}\|$.
In particular, we will make a decision concerning renewals at each iteration
by evaluating the empirically observed decrease in the gradient norm. If
$\|\mathbf{g}_{t}\|-\|\mathbf{g}_{t-1}\|<b(\|\mathbf{g}_{t-1}\|-\|\mathbf{g}_{t-2}\|)$
for some constant $b$, this strategy triggers a renewal. To guarantee at least
linear convergence, we impose a renewal after $n$ iterations in which no
renewal has been triggered. In the rest of the paper we call this strategy GN-
SHED (Gradient Norm-based SHED). See Fig. 3 for an illustration of some
different possible choices for the renewal indices sets.

Figure 3: Illustration of possible choices of renewal indices set. The set
$\mathcal{I}=\\{C_{j}\\}$ specifies the iterations at which a renewal takes
place. (a) illustrates the least squares case in which renewal is performed
only once, see Algorithm 1, (b) the periodic renewals case with $T=5$, so
$\mathcal{I}=\\{1,5,10,...\\}$, see Algorithm 3, and (c) the set in which the
distance between renewals increases according to the Fibonacci sequence.

##  6 Empirical Results

In this section we present empirical results obtained with real datasets. In
particular, to illustrate Algorithm 1, thus in the case of regression via
least squares, we consider the popular Million Song (1M Songs) Dataset Bertin-
Mahieux et al. (2011) and the Online News Popularity dataset Fernandes et al.
(2015), both taken from the UCI Machine Learning Repository Dua and Graff
(2017). For the more general non linear convex case, we apply logistic
regression to two image datasets, in particular the FMNIST Xiao et al. (2017)
and EMNIST digits Cohen et al. (2017) datasets, and to the ‘w8a’ web page
dataset available from the libSVM library Chang and Lin (2011). First, we show
performance assessments related to the theoretical results of the previous
sections. In particular, we show how the versatility of the algorithm applies
effectively to the FL framework, showing results related to parameters choice
and heterogeneity of CRs. Finally, we compare our algorithm against state-of-
the-art approaches in both i.i.d. and non i.i.d. data distributions, showing
the advantage that can be provided by our approach, especially in the case of
heterogeneous CRs. When comparing with other algorithms, we provide results
for both i.i.d. and non i.i.d. partitions to show that our algorithm is
competitive with state-of-the-art approaches also in the i.i.d. configuration.

In the general convex case, we follow the heuristics discussed at the end of
Section 5, where the version using the Fibonacci sequence to define
$\mathcal{I}$ is called Fib-SHED while the event-triggered one based on the
gradient norm inspection is called GN-SHED. As done previously in the paper,
Algorithm 1 is referred to as SHED-LS, and Algorithm 4 as SHED, while their
periodic variants (Algorithm 3), are referred to as SHED-LS-periodic and SHED-
periodic, respectively. If the value of $d_{t}^{(i)}$ is not specified, SHED-
LS and SHED stand for the proposed approach when $d_{t}^{(i)}=1$ for each
iteration $t$ and agent $i$, $i=1,...,M$.

In this section, we show results for $\mu=10^{-5},10^{-6},10^{-8}$, and we
have obtained similar results with $\mu=10^{-4}$. Note that
$10^{-4},10^{-6},10^{-8}$ were the values considered in Wang et al. (2018).

###  6.1 Datasets

We consider the following datasets. For each dataset, we partition the dataset
for the FL framework in both an i.i.d. and non i.i.d. way. In particular:

####  6.1.1 Least squares

We consider the popular Milion Song Dataset Bertin-Mahieux et al. (2011). The
objective of the regression is the prediction of the year in which a song has
been recorded. The dataset consists in 515K dated tracks starting from 1922
recorded songs. We take a subset of 300K tracks as training set. The feature
size of each data sample is $n=90$. As preprocessing of the dataset we apply
rescaling of the features and labels to have them comprised between $0$ and
$1$. To get a number of samples per agent in line with typical federated
learning settings, we spread the dataset to $M=200$ agents, so that each agent
has around $800$ data samples.

As a non i.i.d configuration, we consider a scenario with _label distribution
skew_ and unbalancedness (see Kairouz et al. (2021)). In particular, we
partition the dataset so that each agent has songs from only one year, and
different number of data samples. In this configuration, some agent has a very
small number of data samples, while other agents have up to 11K data samples.

We also provide some results using the Online News Popularity Fernandes et al.
(2015), made of $40$K online news data samples. The target of the prediction
is the number of shares given processed information about the online news.
Data is processed analogously to 1M Songs. $M=30$ agents are considered with
$500$ data sample each.

####  6.1.2 Logistic regression

To validate our algorithm with a non linear convex cost we use two popular and
widely adopted image datasets and a web page prediction dataset.
Image datasets: We consider the FMNIST (Fashion-MNIST) dataset Xiao et al.
(2017), and the EMNIST (Extended-MNIST) digits dataset Cohen et al. (2017).
For both, we consider a training set of $60K$ data samples. We apply standard
pre-processing to the datasets: we rescale the data to be between $0$ and $1$
and apply PCA Sheikh et al. (2020), getting a parameter dimension of $n=300$.
Given that the number of considered data samples is smaller with respect to
1M-Songs, we distribute the dataset to $M=28$ agents.

We focus on one-vs-all binary classification via logistic regression, which is
the basic building block of multinomial logistic regression for multiclass
classification. In a one-vs-all setup, one class (that we call target class),
needs to be distinguished from all the others. The i.i.d. configuration is
obtained by uniformly assigning to each agent samples from the target class
and mixed samples of the other classes. When the agents perform one-vs all, we
provide them with a balanced set where half of the samples belong the target
class, and the other half to the remaining classes, equally mixed. In this
way, each agent has around $400$ data samples.

To partition the FMNIST and EMNIST datasets in a non i.i.d. way, we provide
each agent with samples from only two classes: one is the target class and the
other is one of the other classes. This same approach was considered in Li et
al. (2020b). In this section we show results obtained with a target class
corresponding to label ’one’ of the dataset, but we have seen that equivalent
results are obtained choosing one of the other classes as target class.
Web page dataset: We consider the web page category prediction dataset ‘w8a’,
available from libSVM Chang and Lin (2011). The dataset comprises 50K data
samples and the objective of the prediction is the distinction between two
classes, depending on the website category. The dataset is already pre-
processed. Being a binary and strongly unbalanced dataset, we do not consider
an i.i.d. configuration, and spread the dataset to the agents similarly to the
case of FMNIST, with the difference that some agent only has data samples from
one class. The number of data samples per agent is between $700$ and $900$.

###  6.2 Federated backtracking

To tune the step size when there are no guarantees that a step size equal to
one decreases the cost, we adopt the same strategy adopted in Wang et al.
(2018): an additional communication round takes place in which each agent
shares with the master the loss obtained when the parameter is updated via the
new descent direction for different values of the step size. In this way, we
can apply a distributed version of the popular Armijo backtracking line search
(see Algorithm 2 or Boyd and Vandenberghe (2004), page 464). When showing the
results with respect to communication rounds, we always include also the
additional communication round due to backtracking.

###  6.3 Heterogeneous channels model

As previously outlined, our algorithm is suitable for FL frameworks in which
agents have heterogeneous transmission resources, that is for instance the
case in Federated Edge Learning. As an example we consider the case of
wireless channels, that are characterized by strong variability. In
particular, in the wireless channel subject to fading conditions and
interference, at each communication round agents could have very different CRs
available Pase et al. (2021); Wadu et al. (2020); Chen et al. (2020). To show
the effectiveness of our algorithm in such a scenario, we consider the
Rayleigh fading model. In a scenario with perfect channel knowledge, in which
agent $i$ allocates a certain bandwidth $B^{(i)}=B$ (we consider it to be the
same for all agents) to the FL learning task, we have that the achievable rate
of transmission for the $i$-th user is:

$R^{(i)}=B\log_{2}(1+\gamma\Gamma^{(i)})$ (37)

where $\Gamma^{(i)}$ is a value related to transmission power and
environmental attenuation for user $i$. For simplicity we fix
$\Gamma^{(i)}=\Gamma=5$ for all users (in Pase et al. (2021), for instance,
$\Gamma=1$ and $\Gamma=10$ were considered). The only source of variability is
then $\gamma\sim\mathit{Exp}(\nu)$, modeling the Rayleigh fading effect. We
fix $\nu=1$.

With respect to our algorithm, for illustration purposes, we compare scenarios
in which the number of vectors in $\mathbb{R}^{n}$ related to second order
information, $d_{t}^{(i)}$, is deterministically chosen and fixed for all
users, with a scenario in which $d_{t}^{(i)}$ is randomly chosen for each
round for each user. In particular, we fix:

$d_{t}^{(i)}=d_{\gamma}:=\lfloor d_{0}\log_{2}(1+\gamma\Gamma)\rfloor.$ (38)

This would correspond to each user allocating a bandwidth of $B=d_{0}B_{0}$ to
the transmission of second order information, where $B_{0}$ is the
transmission speed needed to transmit a vector in $\mathbb{R}^{n}$ per channel
use. In the rest of the paper, we fix $d_{0}=2$. This implies that in average
the number of vectors that an agent is able to transmit is $4$, and the actual
number can vary from $0$ (only the gradient is sent) to $10$. For simplicity,
when adopting this framework, we assume that agents are always able to send at
least a vector in $\mathbb{R}^{n}$ per communication round (so the gradient is
always sent). Further simulations are left as future work.

###  6.4 Choice of $\rho_{t}$

In Fig. 4 we analyze the optimization performance with respect to the choice
of the approximation parameter $\rho_{t}$, following the theoretical results
of Section 3.1, in particular Theorem 1 and Corollary 2. We provide
performance analysis in the distributed setting for least squares with the 1M
Song and Online News Popularity datasets in (a) and (b), respectively. We show
the error $\|\theta^{t}-\theta^{*}\|$ versus the number of communication
rounds. We adopt the notation of Section 5 as it includes least squares as a
special case. In this subsection, the same parameter choice is done for all
the agents, thus we omit to specify that the parameter is specific of the
$i$-th agent. We compare two choices for the approximation parameter: (i) the
one providing the best estimation error according to Corollary 2,
$\rho_{t}^{*}=(\hat{\lambda}_{q_{t}+1,t}+\hat{\lambda}_{n,t})/2$ with the step
size $\eta_{t}=1$, and (ii) the one that was proposed in Erdogdu and Montanari
(2015), where $\rho_{t}=\hat{\lambda}_{q_{t}+1,t}$. In the latter case,
following the result of Theorem 1 related to the best convergence rate we pick
the step size to be equal to the $\eta_{t}^{*}$ in Eq. (11). From the plots,
we can see the improvement provided by the choice of $\rho_{t}^{*}$ against
$\hat{\lambda}_{q_{t}+1,t}$, although the performance is very similar for both
values of the approximation parameter.

Figure 4: Performance comparison for different values of $\rho_{t}$. The best
choice in terms of estimation error from Corollary 2 is compared against the
choice that was proposed in Erdogdu and Montanari (2015), that is
$\rho_{t}=\lambda_{q_{t}+1}$

###  6.5 Lyapunov exponent convergence bound

In Fig. 5, we illustrate how the Lyapunov exponent bounds derived in Sec. 4.2
and Theorem 10-2) characterize the linear convergence rate of the algorithms.
In particular, we consider the cases of periodic renewals in both the least
squares and logistic regression case, considering the case in which
$q_{t}^{(i)}=q_{t-1}^{(i)}+1$, $\forall i,t$ and renewals are periodic with
period $T$. The plots show how the linear convergence rate is dominated by the
Lyapunov exponent bound characterized by $\bar{a}_{T}$. For illustrative
purposes, we show the results for the choice of $T=35$ and of $T=25$ for least
squares on 1M songs and logistic regression on FMNIST, respectively.

Figure 5: Linear convergence with periodic renewals illustrated via the study
of the upper bound on the dominant Lyapunov exponent of the estimation error
from equation (23) and point 2) of Theorem 10. $\bar{a}_{T}$ is as in eqs.
(23) and (36) for (a) and (b), respectively. Note that the upper bound is on
the slope of the decreasing cost.

###  6.6 Role of $d_{t}$

In Figure 6, we show the impact of each agent transmitting more eigenvalue-
eigenvector pairs per communication rounds, showing the optimization results
when the increment, $d_{t}^{(i)}$, takes different values. In this subsection,
for illustration purposes, renewals of SHED are determined by the Fibonacci
sequence (see Sec. 5.1). We consider the case when the number is the same and
fixed (specifically we consider $d_{t}^{(i)}\in\\{1,3,6,30\\}$) for all the
agents and the case $d_{t}^{(i)}=d_{\gamma}$, where $d_{\gamma}$ is as defined
in Eq. (38), with $d_{0}=2$ and $\Gamma=5$, that means that in average
$d_{\gamma}$ is equal to $4$. In this latter case each agent is able to
transmit a different random number of eigenvalue-eigenvector pairs. This
configuration is relevant as our algorithm allows agents to contribute to the
optimization according to their specific communication resources (CRs). In
Figures 6-(a)-(c) we show the results for the least squares on 1M Songs, while
in Figures 6-(b)-(d) we show the results in the convex case of logistic
regression on the FMNIST dataset. We can see from Figures 6-(a)-(b) how the
global number of communication rounds needed for convergence can be
significantly reduced by increasing the amount of information transmitted at
each round. In particular, at each round, the number of vectors in
$\mathbb{R}^{n}$ being transmitted is $d_{t}+1$, since together with the
$d_{t}$ scaled eigenvectors,
$\\{\tilde{\mathbf{v}}_{j,t}^{(i)}\\}_{j=q^{(i)}_{t-1}+1}^{q^{(i)}_{t}}$ (see
Algorithm 2), agents need to transmit also the gradient.

Figure 6: Performance comparison for different values of $d_{t}$ for (i)
linear regression on 1M Songs ((a) and (c)) and (ii) logistic regression on
FMNIST ((b) and (d)). $d_{\gamma}$ is as defined in Eq. (38), with $d_{0}=2$,
so $d_{\gamma}$ is in average equal to $4$. In (b), we emphasize three points
where the renewal operation (see Algorithm 4) takes place.

In the least squares case, the number of iterations needed for convergence is
surely smaller than the number of EEPs sent: when the $n-1$-th EEP has been
shared, convergence surely occurs. When the number of EEPs is random
($d_{\gamma}$), but equal to $4$ in average, we see that we can still get a
significant improvement, which is between the choices $d_{t}=3$ and $d_{t}=6$.

Let us now analyze the case of logistic regression. In Figures 6-(b)-(d) we
emphasize the role of the _renewal_ operation, showing also how incrementally
adding eigenvectors of the outdated Hessian improves the convergence, as
formalized and shown in Theorems 9 and 10. In particular, from Figure 6-(b) it
is possible to appreciate the impact of increasing the interval between
renewals.

Quantifying the improvement provided by the usage of greater increments
analyzing the empirical results, we see that when $d_{t}=3$, the number of
vectors transmitted per round is $4$, against the $2$ transmitted when
$d_{t}=1$, so the communication load per round is twice as much. We see that,
despite this increase in data transmitted per round, in the case $d_{t}=3$ the
overall number of communication rounds is halved with respect to the case
$d_{t}=1$. An equivalent result occurs for the case $d_{t}=6$. Notably, for
the case $d_{t}=d_{\gamma}$, with an average increment equal to $4$, we see
that the convergence speed is in between the cases $d_{t}=3$ and $d_{t}=6$,
which shows the effectiveness of the algorithm under heterogeneous channels.
For $d_{t}=30$, even if we increase the communication by a factor $15$, we get
a convergence that is only $10$ times faster. This is of interest because it
shows that if we increase arbitrarily the increment $d_{t}$, we pay in terms
of overall communication load.

In Figures 6-(c)-(d) we plot the error as a function of the amount of data
transmitted per agent, where the number of vectors in $\mathbb{R}^{n}$ is the
unit of measure. These plots show that, for small values of $d_{t}$ (in
particular, $d_{t}\in\\{1,3,6\\}$), the overall data transmitted does not
increase for the considered values of $d_{t}$, meaning that we can
significantly reduce the number of global communication rounds by transmitting
more data per round without increasing the overall communication load. This is
true in particular also for the case $d_{t}=d_{\gamma}$, thus when the agents’
channels availability is heterogeneous at each round, showing that our
algorithm works even in this relevant scenario without increasing the
communication load. On the other hand, in the case of $d_{t}=30$, if we
increase the amount of information by an order of magnitude, even if we get to
faster convergence, we have as a consequence a significant increase in the
communication load that the network has to take care of.

###  6.7 Comparison against other algorithms

In this subsection, SHED-LS+ and SHED+ signifies that the number of Hessian
EEPs is chosen randomly for each agent $i$, according to the model for
simulating heterogeneous CRs illustrated in Sec. 6.3, so
$d_{t}^{(i)}=d_{\gamma}$, with an average increment equal to $4$. To
distinguish the different heuristic strategies for the renewals indices set,
we use Fib-SHED to denote the approaches where $\mathcal{I}$ is pre-determined
by the Fibonacci sequence, while GN-SHED to denote the approach in which
renewals are event-triggered by the gradient norm inspection (see Sec. 5.1)

We compare the performance of our algorithm with a state-of-the art first-
order method, accelerated gradient descent (AGD, with the same implementation
of Wang et al. (2018)).

As benchmark second-order methods we consider the following:

  * •

a distributed version of the Newton-type method proposed in Erdogdu and
Montanari (2015), to which we refer as Mont-Dec, which is the same as
Algorithm 4 with the difference that the renewal occurs at each communication
round, so the Hessian is always recomputed and the second-order information is
never outdated. In this way, the eigenvectors are always the first
($\\{{v_{j,t}}\\}_{j=1}^{d_{t}}$), so this algorithm is neither incremental
nor exploits outdated second-order information. We fix the amount of second-
order information sent by the agents to be the same as the previously
described SHED+.

  * •

the distributed optimization state-of-the-art NT approach exploiting the
harmonic mean, GIANT Wang et al. (2018). While in GIANT the local Newton
direction is computed via conjugate gradient method, in DONE T. Dinh et al.
(2022), a similar approach has been proposed but using Richardson iterations.
To get our results, we provide the compared algorithm with the actual exact
harmonic mean, so that we include both the algorithms at their best. Referring
to this approach, we write GIANT in the figures. We remark that GIANT requires
an extra communication round, because the algorithm requires the agents to get
the global gradient. We also implemented the determinantal averaging approach
Dereziński and Mahoney (2019) to compensate for the inversion bias of the
harmonic mean, but we did not see notable improvements so we do not show its
results for the sake of plot readability.

  * •

the very recently proposed FedNL method presented in Safaryan et al. (2021).
We adopted the best performing variant (FedNL-LS), with rank-1 matrix
compression, which initializes the Hessian at the master with the complete
Hessian computed at the initial parameter value, so the transmission of the
Hessian matrix $\mathbf{H}(\operatorname*{{\boldsymbol{\theta}}}^{0})$ at the
first round is required. This method requires the computation of the Hessian
and of an $n\times n$ matrix SVD at each iteration. FedNL uses the exact same
amount of CRs as SHED (so with $d_{t}^{(i)}=1,\forall i$), as it sends the
gradient and one eigenvector at each iteration. For the Hessian learning step
size, we adopt the same choice of Safaryan et al. (2021) (step size equal to
1), and we observed that different values of this parameter did not provide
improvements.

We do not compare Mont-Dec and FedNL in the least squares case, as these
algorithms are specifically designed for a time-varying Hessian and it does
not make sense to use them in a least squares problem.

Figure 7: Performance comparison for least squares on 1M Songs. Relative cost
is
$f({\operatorname*{{\boldsymbol{\theta}}}}^{t})-f({\operatorname*{{\boldsymbol{\theta}}}}^{*})$.
For SHED-LS and SHED-LS+, see Algorithm 1

In the following, we show the details of results for the least squares on 1M
Songs dataset and then for logistic regression on the FMNIST dataset. The
detailed Figures of convergence results obtained with the EMNIST digits and
‘w8a’ are provided in Appendix B, while a summary of the results obtained on
all the considered datasets is provided in Section 6.7.3.

####  6.7.1 Least squares on 1M Songs

In Fig. 7 we show the performance of the algorithm in the case of least
squares applied on the 1M Songs dataset (thus the procedure is the one
described in Algorithm 1). We plot the relative cost
$f(\theta^{t})-f(\theta^{*})$ versus the communication rounds. We see that AGD
is characterized by very slow convergence due to the high condition number.
Together with AGD, in Fig. 7 we show three Newton-type approaches: the state-
of-the-art GIANT Wang et al. (2018), SHED and SHED+. We see how in the i.i.d.
case GIANT performs well even if SHED+ is the fastest, but all the three
algorithms have similar performance compared to the first-order approach
(AGD). In the non i.i.d. case, whose framework has been described in Sec.
6.1.1, we can see how GIANT has a major performance degradation, even if still
obtaining a better performance with respect to first-order methods. On the
other hand, both SHED and SHED+ do not show performance degradation in the non
i.i.d. case.

####  6.7.2 Logistic regression on FMNIST

Figure 8: Performance comparison of logistic regression on FMNIST when
$\mu=10^{-5}$. Relative cost is
$f({\operatorname*{{\boldsymbol{\theta}}}}^{t})-f({\operatorname*{{\boldsymbol{\theta}}}}^{*})$.

In Figs. 8, 9 we show the results obtained with the FMNIST dataset with the
setup described in Sec. 6.1.2. We show the performance for two values of the
regularization parameter, $\mu$, in particular $\mu=10^{-5}$ and
$\mu=10^{-6}$. We compare our approach (Algorithm 4, SHED) also against FedNL
and Mont-Dec. We show the performance of Fib-SHED, Fib-SHED+ and GN-SHED. As
we did for Fig. 6, to provide a complete comparison, in Figs. 8-9-(c)-(d) we
show also the relative cost versus the overall amount of data transmitted, in
terms of overall number of vectors in $\mathbb{R}^{n}$ transmitted.

Figure 9: Performance comparison of logistic regression on FMNIST when
$\mu=10^{-6}$. Relative cost is
$f({\operatorname*{{\boldsymbol{\theta}}}}^{t})-f({\operatorname*{{\boldsymbol{\theta}}}}^{*})$.

From Fig. 8-(a)-(b), we can see how also in the FMNIST case the non i.i.d.
configuration causes a performance degradation for GIANT, while SHED and FedNL
are not impacted.

In the i.i.d. case, for $\mu=10^{-5}$ (Fig. 8), we see that GN-SHED, Fib-SHED,
Fib-SHED+ and GIANT require a similar number of communication rounds to
converge, while the amount of overall data transmitted per agent is smaller
for GIANT. When $\mu=10^{-6}$, instead, we see from Figure 9 that GIANT
requires more communication rounds to converge while the same amount of
information needs to be transmitted. Hence, in this case GIANT is more
impacted by a lower condition number when compared to our approach. On the
other hand, FedNL in both cases shows a much slower convergence speed with
respect to Fib-SHED and GN-SHED (approaches requiring the same per-iteration
communication load as FedNL). Notice, for FedNL, in Figures 8-(c)-(d) and
9-(c)-(d), the impact that the transmission of the full Hessian matrices at
the first round has on the overall communication load.

In the considered non i.i.d. case, the large advantage that our approach can
provide with respect to the other considered algorithms is strongly evident in
both data transmitted and communication rounds.

Comparing the SHED approaches against Mont-Dec we see the key role that the
_incremental_ strategy exploiting _outdated_ second order information has on
the convergence speed of our approach. Indeed, even though in the first
iterations the usage of the current Hessian information provides the same
performance of the SHED methods, the performance becomes largely inferior in
the following rounds.
From a computational point of view, both the Mont-Dec and the FedNL approach
are much heavier than SHED as they require that each agent recomputes the
Hessian and SVD at each round. SHED, instead, in the more challenging case of
$\mu=10^{-6}$, requires the agents to compute the Hessian matrix only 12 times
out of the 450 rounds needed for convergence. For more details on this, see
Figure 11.

In Fig. 10, we show how SHED is much more resilient to ill-conditioning with
respect to competing algorithms, by comparing the convergence performance when
the regularization parameter is $\mu=10^{-5}$ and $\mu=10^{-8}$. With respect
to FedNL, note how Fib-SHED worsens its performance with $\mu=10^{-8}$ by
being around 2.5 times slower compared to $\mu=10^{-5}$, while FedNL worsens
much more, being more than 8 times slower with $10^{-8}$ compared to
$\mu=10^{-5}$.

Figure 10: Performance of logistic regression on FMNIST, comparing
$\mu=10^{-5}$ and $\mu=10^{-8}$. Relative cost is
$f({\operatorname*{{\boldsymbol{\theta}}}}^{t})-f({\operatorname*{{\boldsymbol{\theta}}}}^{*})$.

####  6.7.3 Comparison on computational complexity

Figure 11: In this plots, we show, for three datasets (FMNIST, EMNIST and
w8a), the number of times an agent is required to compute the local Hessian
matrix in order for an algorithm to converge, comparing the proposed Fib-SHED
and the FedNL Safaryan et al. (2021) algorithms.

In Appendix B, we show the results obtained with the EMNIST and w8a datasets.
We omit the comparison with GN-SHED given that the results are similar to the
ones obtained with Fib-SHED. In the case of EMNIST, we obtain results similar
to FMNIST, except that, when $\mu=10^{-5}$, GIANT is not much impacted by the
considered non i.i.d. configuration. Even if in that case GIANT seems to be
the best choice, all the other results show that GIANT and related approaches
based on the harmonic mean (like DONE) are strongly sensitive to non i.i.d.
data distributions. With image datasets (EMNIST and FMNIST), FedNL is largely
outperformed by our approach (by both Fib-SHED and Fib-SHED+), while, with the
‘w8a’ dataset, FedNL is more competitive. However, the Fib-SHED and Fib-SHED+
approaches have the very appealing feature that they require agents to compute
the local Hessian matrices only _sporadically_. The FedNL approach, instead,
requires that the Hessian is recomputed by agents at each round, implying a
much heavier computational demand. To better illustrate and quantify this
advantage, we show, in Figure 11, the number of times that an agent is
required to compute the local Hessian matrix in order to obtain convergence,
comparing Fib-SHED and FedNL, in the cases of the three datasets. In the case
of EMNIST and FMNIST, we are showing the non i.i.d. configurations, but
similar results can be obtained with the i.i.d. ones. The results show that,
compared with Fib-SHED, the number of times agents are required to compute the
Hessian is always at least ten times greater for FedNL to converge.

##  7 Conclusions

In this work, we have proposed SHED, a Newton-type algorithm to perform FL in
heterogeneous communication networks. SHED is versatile with respect to
agents’ (per-iteration) communication resources and operates effectively in
the presence of non i.i.d. data distributions, outperforming state-of-the-art
techniques. It achieves better performance with respect to the competing FedNL
approach, while involving sporadic Hessian computations. In the case of i.i.d.
data statistics, SHED is also competitive with GIANT, even though the latter
may perform better under certain conditions. We stress that the key advantage
of SHED lies in its robustness under any data distribution, and its
effectiveness and versatility when communication resources differ across nodes
and links.

Future work includes the use and extension of the algorithm for more specific
scenarios and applications, like, for example, wireless networks, and also the
study of new heuristics for the renewal operation in the general convex case.
Furthermore, future research directions may deal with compression techniques
for the Hessian eigenvectors, as well as the development of heuristic
algorithms combining different second-order methods, like FedNL and SHED.
Lastly, the spectral characteristics of the data could also be exploited to
tailor the proposed algorithm to specific FL problems, leading to further
benefits in terms of convergence speed.

Acknowledgments

This work has been supported, in part, by the Italian Ministry of Education,
University and Research, through the PRIN project no. 2017NS9FEY entitled
“Realtime Control of 5G Wireless Networks: Taming the Complexity of Future
Transmission and Computation Challenges”.

## Appendix A.

## Proof of Theorem 8

First, we need to prove the following Lemma:

###### Lemma 11

If $\|\mathbf{g}_{t}\|>\omega>0$, for $\hat{\mathbf{H}}_{t}$ defined in (26),
$\mathbf{p}_{t}=\hat{\mathbf{H}}_{t}^{-1}\mathbf{g}_{t}$ there are
$\gamma_{t},\eta_{t}>0$ such that

$f(\boldsymbol{\theta}^{t}-\eta_{t}\mathbf{p}_{t})\leq
f(\boldsymbol{\theta}^{t})-\gamma_{t},$ (39)

in particular, for a backtracking line search with parameters
$\alpha\in(0,0.5),\beta\in(0,1)$, it holds:

$\gamma_{t}=\alpha\beta\frac{\bar{\rho}_{t}}{K^{2}}\omega^{2}$ (40)

Proof  The proof is the same as the one provided in Boyd and Vandenberghe
(2004) for the damped Newton phase (page 489-490), with the difference that
the ”Newton decrement”, that we denote by $\sigma_{t}$, here is
$\sigma_{t}^{2}:=\mathbf{g}_{t}^{T}\mathbf{p}_{t}=\mathbf{g}_{t}^{T}\hat{\mathbf{H}}_{t}^{-1}\mathbf{g}_{t}=\mathbf{p}_{t}^{T}\hat{\mathbf{H}}_{t}\mathbf{p}_{t}$.
Furthermore, the property $\hat{\mathbf{H}}_{t}\geq\bar{\rho}_{t}\mathbf{I},\
\forall{t}$ is used in place of strong convexity.

Note that, because of Assumption 1, it always holds that $\rho_{t}^{(i)}>0,\ \
\forall{i},\forall{t}$, and this implies $\bar{\rho}_{t}>0,\forall{t}$. When
$\bar{\rho}_{t}>0,\forall{t}$, Lemma 11 implies $\|\mathbf{g}_{t}\|\rightarrow
0$. Indeed, if not, there would be some $\epsilon>0$ such that
$\|\mathbf{g}_{t}\|>\epsilon,\ \forall{t}$. But then (39) immediately implies
that $f(\boldsymbol{\theta}^{t})\rightarrow-\infty$, that contradicts the
strong convexity hypothesis. By strong convexity and differentiability of $f$,
$g(\boldsymbol{\theta})=\nabla
f(\boldsymbol{\theta})=0\implies\boldsymbol{\theta}=\boldsymbol{\theta}^{*}$.

## Proof of Theorem 10

1) We first need the following Lemma:

###### Lemma 12

Let $\bar{\kappa}=\sum_{i=1}^{M}\kappa_{i}$, with $\kappa_{i}$ the strong
convexity constant of the cost $f^{(i)}$ of agent $i$, let $K$ be the
smoothness constant of $f$ and
$M(t)=\max\\{\|{\operatorname*{{\boldsymbol{\theta}}}}^{t}-{\operatorname*{{\boldsymbol{\theta}}}}^{*}\|,\|{\operatorname*{{\boldsymbol{\theta}}}}^{k_{t}}-{\operatorname*{{\boldsymbol{\theta}}}}^{*}\|\\}$.
Applying Algorithm 4, if

$3\bar{\kappa}(M(t)+\|{\operatorname*{{\boldsymbol{\theta}}}}^{t}-{\operatorname*{{\boldsymbol{\theta}}}}^{*}\|)+K\|{\operatorname*{{\boldsymbol{\theta}}}}^{t}-{\operatorname*{{\boldsymbol{\theta}}}}^{*}\|\leq\frac{3\bar{\kappa}^{2}}{L}(1-2\alpha),$
(41)

then, for any $\alpha\in(0,1/2)$, the backtracking algorithm (2) chooses
$\eta_{t}=1$.

Proof [Sketch of proof] The proof uses a procedure similar to the one used in
Boyd and Vandenberghe (2004), page 490-491 to prove the beginning of the
quadratically convergent phase in the centralized Newton method. Then,
leveraging Lipschitz continuity and the fact that, thanks to the choice
$\rho_{t}^{(i)}=\lambda_{q_{t}^{(i)}}^{(i)}$, it is
$\hat{\mathbf{H}}_{t}\geq\mathbf{H}({\operatorname*{{\boldsymbol{\theta}}}}^{k_{t}})$,
the result is derived. For the complete proof, see the Proof of Lemma 12 at
the end of this Appendix.

Let condition (41) be satisfied. Then, if

$\vspace{-0.3cm}(3/2)L\|{\operatorname*{{\boldsymbol{\theta}}}}^{t}-{\operatorname*{{\boldsymbol{\theta}}}}^{*}\|+LM(t)\leq\bar{\kappa},$
(42)

the convergence of SHED is at least linear. Indeed, if (41) is satisfied,
then, from Lemma 12, the step size is $\eta_{t}=1$ and the convergence bound
(see Theorem 9) becomes
$\|\boldsymbol{\theta}^{t+1}-\boldsymbol{\theta}^{*}\|\leq
c_{t}\|\boldsymbol{\theta}^{t}-\boldsymbol{\theta}^{*}\|$, with

$\displaystyle c_{t}$
$\displaystyle=(1-\frac{\bar{\lambda}_{n,t}}{\bar{\rho}_{t}}+\frac{L}{\bar{\rho}_{t}}\|\boldsymbol{\theta}^{t}-\boldsymbol{\theta}^{k_{t}}\|+\frac{L}{2\bar{\rho}_{t}}\|\boldsymbol{\theta}^{t}-\boldsymbol{\theta}^{*}\|)$
(43)
$\displaystyle\leq(1-\frac{\bar{\lambda}_{n,t}}{\bar{\rho}_{t}}+\frac{L}{\bar{\rho}_{t}}\|\boldsymbol{\theta}^{k_{t}}-\boldsymbol{\theta}^{*}\|+\frac{3L}{2\bar{\rho}_{t}}\|\boldsymbol{\theta}^{t}-\boldsymbol{\theta}^{*}\|)$

and it is easy to see that condition (42) implies that $c_{t}<1$ and thus we
get a contraction in
$\|{\operatorname*{{\boldsymbol{\theta}}}}^{t}-{\operatorname*{{\boldsymbol{\theta}}}}^{*}\|$.
Furthermore, when conditions (41) and (42) are both satisfied at some
iteration $\bar{t}$, they are then satisfied for all $t\geq\bar{t}$ and thus
$c_{t}<1$ for all $t\geq\bar{t}$. Indeed, $c_{t}<1$ implies that
$\|{\operatorname*{{\boldsymbol{\theta}}}}^{t+1}-{\operatorname*{{\boldsymbol{\theta}}}}^{*}\|<\|{\operatorname*{{\boldsymbol{\theta}}}}^{t}-{\operatorname*{{\boldsymbol{\theta}}}}^{*}\|$
and $M(t+1)\leq M(t)$ because either $k_{t+1}=k_{t}$ or $k_{t+1}=t+1$. Note
that Assumption 3 is needed to guarantee that (41) and (42) are eventually
satisfied.

2) From 1), we can write
$\|\operatorname*{{\boldsymbol{\theta}}}^{t}-\operatorname*{{\boldsymbol{\theta}}}^{*}\|\leq
C{a}^{t}$ for some $a\in(0,1)$ and some $C>0$. Considering $t\geq\bar{t}$,
with $\bar{t}$ the first iteration for which both conditions (41) and (42) are
satisfied, we consider ${c}_{t}$ as in (43), and let $T=t-k_{t}$

$\displaystyle{c}_{t}$ $\displaystyle\leq
1-\frac{\bar{\lambda}_{n,t}}{\bar{\rho}_{t}}+\frac{L}{\bar{\rho}_{t}}\|\boldsymbol{\theta}^{k_{t}}-\boldsymbol{\theta}^{*}\|+\frac{3L}{2\bar{\rho}_{t}}\|\boldsymbol{\theta}^{t}-\boldsymbol{\theta}^{*}\|$
(44) $\displaystyle\leq
1-\frac{\bar{\lambda}_{n,t}}{\bar{\rho}_{t}}+\frac{L}{\bar{\rho}_{t}}C_{1}{a_{*}}^{t-T}+\frac{3L}{2\bar{\rho}_{t}}C_{2}{a_{*}}^{t}$
$\displaystyle=1-\frac{\bar{\lambda}_{n,t}}{\bar{\rho}_{t}}+Ba_{*}^{t},$

where
$B=\frac{L}{\bar{\rho}_{t}}C_{1}a_{*}^{-T}+\frac{3L}{2\bar{\rho}_{t}}C_{2}$,
and $C_{1},C_{2}$ are some bounded positive constants. For any iteration $t$,
we can write
$\|\operatorname*{{\boldsymbol{\theta}}}^{t+1}-\operatorname*{{\boldsymbol{\theta}}}^{*}\|\leq\bar{c}_{t}\|\operatorname*{{\boldsymbol{\theta}}}^{t}-\operatorname*{{\boldsymbol{\theta}}}^{*}\|$
for some $\bar{c}_{t}$ that could also be greater than one, if $t<\bar{t}$,
but it is easy to see that $\bar{c}_{t}$ is always bounded. Now, we consider
$\bar{a}:=\limsup_{t}{(\prod_{k=1}^{t}{\bar{c}_{k}})^{1/t}}$. It is
straightforward to see that, as in the least squares case, it is
$a_{*}\leq\bar{a}$. We can write
$\log{\bar{a}}=\limsup_{t}{\frac{1}{t}}(\sum_{k=\bar{t}+1}^{t}\log{{c}_{k}})$.
We get

$\displaystyle\log{\bar{a}}$
$\displaystyle\leq\limsup_{t}{\frac{1}{t}}(\sum_{k=\bar{t}+1}^{t}\log{{c}_{k}})$
$\displaystyle\leq\limsup_{t}\frac{1}{t}\sum_{k=\bar{t}+1}^{t}\log{(1-\frac{\bar{\lambda}_{n,k}}{\bar{\rho}_{k}}+B{a_{*}}^{k})}$
$\displaystyle=\limsup_{t}\frac{1}{t}\sum_{k=\bar{t}+1}^{t}\log{(1-\frac{\bar{\lambda}_{n,k}}{\bar{\rho}_{k}})}+\log{(1+\frac{B{a_{*}}^{k}}{1-\frac{\bar{\lambda}_{n,k}}{\bar{\rho}_{k}}})}.$

We see that the last term is

$\vspace{-0.2cm}\log{(1+\frac{B{a_{*}}^{k}}{1-\frac{\bar{\lambda}_{n,k}}{\bar{\rho}_{k}}})}\leq\frac{B{a_{*}}^{k}}{1-\frac{\bar{\lambda}_{n,k}}{\bar{\rho}_{k}}}\leq{\bar{B}{a_{*}}^{k}}$

that comes from the identity $\log{(1+x)}\leq x$, and where
$\bar{B}=\max_{k}\frac{1}{{1-\frac{\bar{\lambda}_{n,k}}{\bar{\rho}_{k}}}}$,
bounded because $|\mathcal{X}_{t}=0|,\ \ \forall{t}$, and thus
$\bar{\rho}_{k}>\bar{\lambda}_{n,k},\forall k$. Now, we see that

$\vspace{-0.3cm}\limsup_{t}{\frac{1}{t}}\sum_{k=1}^{t}{\bar{B}{a_{*}}^{k}}=\limsup_{t}\frac{1}{t}{\bar{B}}(\frac{1-{a_{*}}^{t+1}}{1-a_{*}}-1)=0.$

Hence, we get, using also the finiteness of $\bar{t}$,

$\vspace{-0.2cm}\log{\bar{a}}=\limsup_{t}\frac{1}{t}\sum_{k=1}^{t}\log{(1-\frac{\bar{\lambda}_{n,k}}{\bar{\rho}_{k}})}.$
(45)

Now we use local Lipschitz continuity (Assumption 1) to conclude the proof.
Lipschitz continuity of $f^{(i)}(\boldsymbol{\theta})$ implies that, for any
$k\in\\{1,...,n\\}$,
$|\lambda_{k}^{(i)}(\operatorname*{{\boldsymbol{\theta}}})-\lambda_{k}^{(i)}(\operatorname*{{\boldsymbol{\theta}}}^{*})|\leq\bar{L}\|\operatorname*{{\boldsymbol{\theta}}}-\operatorname*{{\boldsymbol{\theta}}}^{*}\|$,
which in turn implies
$\lambda_{k}^{(i)}(\operatorname*{{\boldsymbol{\theta}}})\geq\lambda_{k}^{(i)}(\operatorname*{{\boldsymbol{\theta}}}^{*})-\bar{L}\|\operatorname*{{\boldsymbol{\theta}}}-\operatorname*{{\boldsymbol{\theta}}}^{*}\|$
and
$\lambda_{k}^{(i)}(\operatorname*{{\boldsymbol{\theta}}})\leq\lambda_{k}^{(i)}(\operatorname*{{\boldsymbol{\theta}}}^{*})+\bar{L}\|\operatorname*{{\boldsymbol{\theta}}}-\operatorname*{{\boldsymbol{\theta}}}^{*}\|$.
For a proof of this result, see also Bisgard (2020), page 116, Theorem 4.25.
It follows that

$\displaystyle\log{\bar{a}}$
$\displaystyle\leq\limsup_{t}\frac{1}{t}\sum_{k=1}^{t}\log{(1-\frac{\bar{\lambda}_{n}^{o}-\bar{L}\|\operatorname*{{\boldsymbol{\theta}}}^{k_{t}}-\operatorname*{{\boldsymbol{\theta}}}^{*}\|}{\bar{\rho}_{k}^{o}+\bar{L}\|\operatorname*{{\boldsymbol{\theta}}}^{k_{t}}-\operatorname*{{\boldsymbol{\theta}}}^{*}\|})}$
$\displaystyle\stackrel{{\scriptstyle(1)}}{{=}}\limsup_{t}\frac{1}{t}\sum_{k=1}^{t}\log{(1-\frac{\bar{\lambda}_{n}^{o}}{\bar{\rho}_{k}^{o}+\bar{L}\|\operatorname*{{\boldsymbol{\theta}}}^{k_{t}}-\operatorname*{{\boldsymbol{\theta}}}^{*}\|})}$
$\displaystyle=\limsup_{t}\frac{1}{t}\sum_{k=1}^{t}\log{(\frac{\bar{\rho}_{k}^{o}-\bar{\lambda}_{n}^{o}+\bar{L}\|\operatorname*{{\boldsymbol{\theta}}}^{k_{t}}-\operatorname*{{\boldsymbol{\theta}}}^{*}\|}{\bar{\rho}_{k}^{o}+\bar{L}\|\operatorname*{{\boldsymbol{\theta}}}^{k_{t}}-\operatorname*{{\boldsymbol{\theta}}}^{*}\|})}$
$\displaystyle\leq\limsup_{t}\frac{1}{t}\sum_{k=1}^{t}\log{(1-\frac{\bar{\lambda}_{n}^{o}}{{\bar{\rho}_{k}^{o}}}+\frac{\bar{L}\|\operatorname*{{\boldsymbol{\theta}}}^{k_{t}}-\operatorname*{{\boldsymbol{\theta}}}^{*}\|}{\bar{\rho}_{k}^{o}})}$
$\displaystyle\stackrel{{\scriptstyle(2)}}{{=}}\limsup_{t}\frac{1}{t}\sum_{k=1}^{t}\log{(1-\frac{\bar{\lambda}_{n}^{o}}{{\bar{\rho}_{k}^{o}}})}$

where equalities (1) and (2) follow from calculations equivalent to the ones
used to obtain (45).
3) Consider $c_{t}$ as it was defined and bounded in (44), $C_{3}>0$ a
constant such that $\log{c_{k}}\leq C_{3},\ \forall k$, and such that
$\log{B}\leq C_{3}$. Let $\bar{a}$ be defined as before. Let $a\in(0,1)$. We
have

$\displaystyle\log{\bar{a}}$
$\displaystyle=\limsup_{t}\frac{1}{t}\sum_{k=1}^{t}\log{c_{k}}$ (46)
$\displaystyle\leq\limsup_{t}\frac{1}{t}\sum_{k=1}^{t}\log{(1-\frac{\bar{\lambda}_{n,k}}{\bar{\rho}_{k}}+Ba^{k})}$
$\displaystyle=\limsup_{t}\frac{1}{t}(\sum_{k\notin\mathcal{X}_{t}}\log{c_{k}}+\sum_{k\in\mathcal{X}_{t}}\log{Ba^{k}})$
$\displaystyle\leq
2C_{3}+\limsup_{t}\frac{1}{t}\log{a}\sum_{k\in\mathcal{X}_{t}}k$
$\displaystyle\leq 2C_{3}+\limsup_{t}\frac{1}{t}\log
aC_{4}|\mathcal{X}_{t}|^{2}$ $\displaystyle=2C_{3}+\limsup_{t}\log
aC_{4}{h(t)}^{2}=-\infty$

where $C_{4}>0$ is some positive constant and the last equality follows
because $\log a<0$ and $\lim_{t}h(t)=\infty$. We see that $0\leq
a_{*}\leq\bar{a}\leq 0$, which implies $a_{*}=0$..

## Proof of Lemma 12

The following proof is similar to the proof for the beginning of the
quadratically convergent phase in centralized Newton method by Boyd and
Vandenberghe (2004), page 490-491. We start with some definitions: at
iteration $t$, let
$f({\operatorname*{{\boldsymbol{\theta}}}}^{t}),\mathbf{g}_{t}=\nabla
f({\operatorname*{{\boldsymbol{\theta}}}}^{t}),\mathbf{H}_{t}=\nabla^{2}f({\operatorname*{{\boldsymbol{\theta}}}}^{t})$
be the cost, the gradient and the Hessian, respectively, computed at
${\operatorname*{{\boldsymbol{\theta}}}}^{t}$. Let
$\mathbf{H}_{k_{t}}=\nabla^{2}f({\operatorname*{{\boldsymbol{\theta}}}}^{k_{t}})$
be the Hessian at ${\operatorname*{{\boldsymbol{\theta}}}}^{k_{t}}$, and
$\hat{\mathbf{H}}_{t}$ the global Hessian approximation of Algorithm 4. Let
the NT descent direction be
$\mathbf{p}_{t}=\hat{\mathbf{H}}_{t}^{-1}\mathbf{g}_{t}$. Define

$\displaystyle\bar{\sigma}_{t}^{2}$
$\displaystyle:=\mathbf{p}_{t}^{T}\mathbf{g}_{t}=\mathbf{g}_{t}^{T}\hat{\mathbf{H}}_{t}^{-1}\mathbf{g}_{t},$
(47) $\displaystyle{\sigma}_{t}^{2}$
$\displaystyle:=\mathbf{p}_{t}^{T}{\nabla^{2}f({\operatorname*{{\boldsymbol{\theta}}}}^{t})}\mathbf{p}_{t}=\mathbf{p}_{t}^{T}{\mathbf{H}}_{t}\mathbf{p}_{t},$
$\displaystyle\tilde{f}(\eta)$
$\displaystyle:=f({\operatorname*{{\boldsymbol{\theta}}}}^{t}-\eta\mathbf{p}_{t}),\
\ \tilde{f}(0)=f({\operatorname*{{\boldsymbol{\theta}}}}^{t}),$
$\displaystyle\tilde{f}^{\prime}(\eta)$
$\displaystyle:=\frac{\partial\tilde{f}(\eta)}{\partial\eta}=-\nabla{f}({\operatorname*{{\boldsymbol{\theta}}}}^{t}-\eta\mathbf{p}_{t})^{T}\mathbf{p}_{t},$
$\displaystyle\tilde{f}^{\prime\prime}(\eta)$
$\displaystyle:=\frac{\partial^{2}\tilde{f}(\eta)}{\partial\eta^{2}}=\mathbf{p}_{t}^{T}\nabla^{2}f({\operatorname*{{\boldsymbol{\theta}}}}^{t}-\eta\mathbf{p}_{t})\mathbf{p}_{t}.$

Note that

$\displaystyle\tilde{f}^{\prime}(0)$
$\displaystyle=-\mathbf{g}_{t}^{T}\mathbf{p}_{t}=-\bar{\sigma}_{t}^{2},$ (48)
$\displaystyle\tilde{f}^{\prime\prime}(0)$
$\displaystyle=\mathbf{p}_{t}^{T}\mathbf{H}_{t}\mathbf{p}_{t}=\sigma_{t}^{2}.$

Note that $\hat{\mathbf{H}}_{t}\geq\bar{\kappa},\ \ \forall{t}$. Thanks to
$L$-Lipschitz continuity, it holds that
$\|\nabla^{2}f({\operatorname*{{\boldsymbol{\theta}}}}^{t}-\eta\mathbf{p}_{t})-\nabla^{2}f({\operatorname*{{\boldsymbol{\theta}}}}^{t})\|\leq\eta
L\|\mathbf{p}_{t}\|$ and we have that

$\displaystyle|\tilde{f}^{\prime\prime}(\eta)-\tilde{f}^{\prime\prime}(0)|$
$\displaystyle=\mathbf{p}_{t}^{T}(\nabla^{2}f({\operatorname*{{\boldsymbol{\theta}}}}^{t}-\eta\mathbf{p}_{t})-\nabla^{2}f({\operatorname*{{\boldsymbol{\theta}}}}^{t}))\mathbf{p}_{t}$
(49) $\displaystyle\leq\eta L\|\mathbf{p}_{t}\|^{3}\leq\eta
L\frac{\bar{\sigma}^{3}({\operatorname*{{\boldsymbol{\theta}}}}^{t})}{\bar{\kappa}^{3/2}}$

where the last inequality holds because
$\bar{\kappa}\|\mathbf{p}_{t}\|^{2}\leq\mathbf{p}_{t}^{T}\hat{\mathbf{H}}\mathbf{p}_{t}=\bar{\sigma}_{t}^{2}$.
From (49) it follows that

$\tilde{f}^{\prime\prime}(\eta)\leq\tilde{f}^{\prime\prime}(0)+\eta
L\frac{\bar{\sigma}^{3}({\operatorname*{{\boldsymbol{\theta}}}}^{t})}{\bar{\kappa}^{3/2}}=\sigma_{t}^{2}+\eta
L\frac{\bar{\sigma}^{3}({\operatorname*{{\boldsymbol{\theta}}}}^{t})}{\bar{\kappa}^{3/2}}.$
(50)

Similarly to Boyd and Vandenberghe (2004), page 490-491, we can now integrate
both sides of the inequality getting

$\tilde{f}^{\prime}(\eta)\leq\tilde{f}^{\prime}(0)+\eta\sigma_{t}^{2}+\frac{\eta^{2}}{2}L\frac{\bar{\sigma}_{t}^{3}}{\bar{\kappa}^{3/2}}=-\bar{\sigma}_{t}^{2}+\eta\sigma_{t}^{2}+\frac{\eta^{2}}{2}L\frac{\bar{\sigma}_{t}^{3}}{\bar{\kappa}^{3/2}}.$

By integrating again both sides of the inequality, we get, recalling that
$\tilde{f}(0)=f({\operatorname*{{\boldsymbol{\theta}}}}^{t})$,

$\displaystyle\tilde{f}(\eta)$
$\displaystyle=f({\operatorname*{{\boldsymbol{\theta}}}}^{t}-\eta\mathbf{p}_{t})\leq
f({\operatorname*{{\boldsymbol{\theta}}}}^{t})-\eta\bar{\sigma}_{t}^{2}+\frac{\eta^{2}}{2}\sigma_{t}^{2}+\frac{\eta^{3}}{6}L\frac{\bar{\sigma}_{t}^{3}}{\bar{\kappa}^{3/2}}.$
(51)

Now, recalling
$\mathbf{H}({\operatorname*{{\boldsymbol{\theta}}}}^{k_{t}})\leq\hat{\mathbf{H}}_{t}$,
we get that

$\displaystyle\sigma_{t}^{2}$
$\displaystyle=\mathbf{p}_{t}^{T}\mathbf{H}_{t}\mathbf{p}_{t}=\mathbf{p}_{t}^{T}(\mathbf{H}_{k_{t}}+\mathbf{H}_{t}-\mathbf{H}_{k_{t}})\mathbf{p}_{t}$
(52)
$\displaystyle=\mathbf{p}_{t}^{T}\mathbf{H}_{k_{t}}\mathbf{p}_{t}+\mathbf{p}_{t}^{T}(\mathbf{H}_{t}-\mathbf{H}_{k_{t}})\mathbf{p}_{t}$
$\displaystyle\leq\mathbf{p}_{t}^{T}\hat{\mathbf{H}}_{t}\mathbf{p}_{t}+L\|\mathbf{p}_{t}\|^{2}\|{\operatorname*{{\boldsymbol{\theta}}}}^{t}-{\operatorname*{{\boldsymbol{\theta}}}}^{k_{t}}\|\leq\bar{\sigma}_{t}^{2}+\frac{L\bar{\sigma}_{t}^{2}}{\bar{\kappa}}\|{\operatorname*{{\boldsymbol{\theta}}}}^{t}-{\operatorname*{{\boldsymbol{\theta}}}}^{k_{t}}\|,$

where we have used Lipschitz continuity, and the fact that
$\mathbf{p}_{t}^{T}\hat{\mathbf{H}}_{t}\mathbf{p}_{t}=\bar{\sigma}_{t}^{2}$
and $\bar{\kappa}\|\mathbf{p}_{t}\|^{2}\leq{\bar{\sigma}_{t}^{2}}$.
Next, setting $\eta=1$ and plugging (52) in (51), we get

$\displaystyle f({\operatorname*{{\boldsymbol{\theta}}}}^{t}-\mathbf{p}_{t})$
$\displaystyle\leq
f({\operatorname*{{\boldsymbol{\theta}}}}^{t})-\bar{\sigma}_{t}^{2}+\frac{\sigma_{t}^{2}}{2}+\frac{L}{6}\frac{\bar{\sigma}_{t}^{3}}{\bar{\kappa}^{3/2}}$
(53) $\displaystyle\leq
f({\operatorname*{{\boldsymbol{\theta}}}}^{t})-\frac{\bar{\sigma}_{t}^{2}}{2}+\frac{L\bar{\sigma}_{t}^{2}}{2\bar{\kappa}}\|{\operatorname*{{\boldsymbol{\theta}}}}^{t}-{\operatorname*{{\boldsymbol{\theta}}}}^{k_{t}}\|+\frac{L}{6}\frac{\bar{\sigma}_{t}^{3}}{\bar{\kappa}^{3/2}}$
$\displaystyle\leq
f({\operatorname*{{\boldsymbol{\theta}}}}^{t})-{\bar{\sigma}_{t}^{2}}(\frac{1}{2}-\frac{L}{2\bar{\kappa}}\|{\operatorname*{{\boldsymbol{\theta}}}}^{t}-{\operatorname*{{\boldsymbol{\theta}}}}^{k_{t}}\|-\frac{L}{6}\frac{\bar{\sigma}_{t}}{\bar{\kappa}^{3/2}})$
$\displaystyle=f({\operatorname*{{\boldsymbol{\theta}}}}^{t})-\mathbf{p}_{t}^{T}\mathbf{g}_{t}(\frac{1}{2}-\frac{L}{2\bar{\kappa}}\|{\operatorname*{{\boldsymbol{\theta}}}}^{t}-{\operatorname*{{\boldsymbol{\theta}}}}^{k_{t}}\|-\frac{L}{6}\frac{\bar{\sigma}_{t}}{\bar{\kappa}^{3/2}}).$

In order for (53) to satisfy the Armijo-Goldstein condition (18) for any
parameter $\alpha\in(0,1/2)$ we see that

$\frac{1}{2}-\frac{L}{2\bar{\kappa}}\|{\operatorname*{{\boldsymbol{\theta}}}}^{t}-{\operatorname*{{\boldsymbol{\theta}}}}^{k_{t}}\|-\frac{L}{6}\frac{\bar{\sigma}_{t}}{\bar{\kappa}^{3/2}}\geq\alpha$
(54)

provides a sufficient condition. The above inequality can be written as

$3\bar{\kappa}\|{\operatorname*{{\boldsymbol{\theta}}}}^{t}-{\operatorname*{{\boldsymbol{\theta}}}}^{k_{t}}\|+\bar{\kappa}^{1/2}{\bar{\sigma}_{t}}\leq\frac{3\bar{\kappa}^{2}}{L}(1-2\alpha).$
(55)

We have that
$\bar{\sigma}_{t}^{2}=\mathbf{g}_{t}^{T}\hat{\mathbf{H}}^{-1}\mathbf{g}_{t}\leq\|\mathbf{g}_{t}\|^{2}/\bar{\kappa}$,
which implies $\bar{\kappa}^{1/2}{\|}\sigma_{t}{\|}\leq\|\mathbf{g}_{t}\|$.
Furthermore, by triangular inequality, we have
$\|{\operatorname*{{\boldsymbol{\theta}}}}^{t}-{\operatorname*{{\boldsymbol{\theta}}}}^{k_{t}}\|\leq\|{\operatorname*{{\boldsymbol{\theta}}}}^{t}-{\operatorname*{{\boldsymbol{\theta}}}}^{*}\|+\|{\operatorname*{{\boldsymbol{\theta}}}}^{k_{t}}-{\operatorname*{{\boldsymbol{\theta}}}}^{*}\|$.
Therefore, we see that if

$3\bar{\kappa}(\|{\operatorname*{{\boldsymbol{\theta}}}}^{t}-{\operatorname*{{\boldsymbol{\theta}}}}^{*}\|+\|{\operatorname*{{\boldsymbol{\theta}}}}^{k_{t}}-{\operatorname*{{\boldsymbol{\theta}}}}^{*}\|)+\|\mathbf{g}_{t}\|\leq\frac{3\bar{\kappa}^{2}}{L}(1-2\alpha),$
(56)

then the Armijo-Goldstein condition is satisfied and $\eta=1$ is chosen by the
backtracking algorithm, proving the Lemma. Indeed, by $K$-smoothness of the
cost function (see Assumption 1) we have $\|\mathbf{g}_{t}\|\leq
K\|{\operatorname*{{\boldsymbol{\theta}}}}^{t}-{\operatorname*{{\boldsymbol{\theta}}}}^{*}\|$
and so the condition of the Lemma implies (56).

## Appendix B.

## Results on EMNIST digits and w8a

In this appendix, we include the results on the EMNIST digits and ‘w8a’
datasets when comparing the different algorithms. We show results for two
values of the regularization parameter $\mu$, specifically $\mu=10^{-5}$ and
$\mu=10^{-6}$, in Figure 12 and 13, respectively. The results obtained on the
EMNIST digits dataset confirm the results that were obtained with FMNIST, with
the difference that in the case $\mu=10^{-5}$, GIANT is not much impacted by
the considered non i.i.d. configuration. The results obtained with the ‘w8a’
dataset show that, while GIANT performance is largely degraded because of the
non i.i.d. configuration, also in this case Fib-SHED and Fib-SHED+
significantly outperform FedNL in both communication rounds and communication
load required for convergence.

Figure 12: Performance comparison of logistic regression on EMNIST when
$\mu=10^{-5}$. Relative cost is $f(\theta^{t})-f(\theta^{*})$. Figure 13:
Performance comparison of logistic regression on EMNIST when $\mu=10^{-6}$.
Relative cost is $f(\theta^{t})-f(\theta^{*})$. Figure 14: Performance
comparison of logistic regression on w8a when $\mu=10^{-5}$ and $\mu=10^{-6}$.
Relative cost is $f(\theta^{t})-f(\theta^{*})$.

## References

  * Agafonov et al. (2022) Artem Agafonov, Dmitry Kamzolov, Rachael Tappenden, Alexander Gasnikov, and Martin Takáč.  Flecs: A federated learning second-order framework via compression and sketching.  _arXiv preprint arXiv:2206.02009_ , 2022.
  * Alimisis et al. (2021) Foivos Alimisis, Peter Davies, and Dan Alistarh.  Communication-efficient distributed optimization with quantized preconditioners.  In Marina Meila and Tong Zhang, editors, _Proceedings of the 38th International Conference on Machine Learning_ , volume 139 of _Proceedings of Machine Learning Research_ , pages 196–206. PMLR, 18–24 Jul 2021.
  * Amiri and Gündüz (2020) Mohammad Mohammadi Amiri and Deniz Gündüz.  Federated learning over wireless fading channels.  _IEEE Transactions on Wireless Communications_ , 19(5):3546–3557, 2020.
  * Bertin-Mahieux et al. (2011) Thierry Bertin-Mahieux, Daniel P.W. Ellis, Brian Whitman, and Paul Lamere.  The million song dataset.  In _Proceedings of the 12th International Conference on Music Information Retrieval (ISMIR 2011)_ , 2011.
  * Bisgard (2020) James Bisgard.  _Analysis and Linear Algebra: The Singular Value Decomposition and Applications_ , volume 94.  American Mathematical Soc., 2020.
  * Boyd and Vandenberghe (2004) Stephen Boyd and Lieven Vandenberghe.  _Convex optimization_.  Cambridge university press, 2004.
  * Chang and Lin (2011) Chih-Chung Chang and Chih-Jen Lin.  LIBSVM: A library for support vector machines.  _ACM Transactions on Intelligent Systems and Technology_ , 2:27:1–27:27, 2011.  Software available at http://www.csie.ntu.edu.tw/~cjlin/libsvm.
  * Chen et al. (2020) Mingzhe Chen, Zhaohui Yang, Walid Saad, Changchuan Yin, H Vincent Poor, and Shuguang Cui.  A joint learning and communications framework for federated learning over wireless networks.  _IEEE Transactions on Wireless Communications_ , 20(1):269–283, 2020.
  * Chen et al. (2018) Tianyi Chen, Georgios Giannakis, Tao Sun, and Wotao Yin.  Lag: Lazily aggregated gradient for communication-efficient distributed learning.  _Advances in neural information processing systems_ , 2018.
  * Cohen et al. (2017) Gregory Cohen, Saeed Afshar, Jonathan Tapson, and André van Schaik.  Emnist: Extending mnist to handwritten letters.  In _2017 International Joint Conference on Neural Networks (IJCNN)_ , pages 2921–2926, 2017.  doi: 10.1109/IJCNN.2017.7966217.
  * Crane and Roosta (2019) Rixon Crane and Fred Roosta.  DINGO: Distributed newton-type method for gradient-norm optimization.  _Advances in Neural Information Processing Systems_ , 32, 2019.
  * Czornik et al. (2012) Adam Czornik, Aleksander Nawrat, Michal Niezabitowski, and Aneta Szyda.  On the lyapunov and bohl exponent of time-varying discrete linear system.  In _2012 20th Mediterranean Conference on Control Automation (MED)_ , 2012.
  * Dereziński and Mahoney (2019) Michał Dereziński and Michael W. Mahoney.  Distributed estimation of the inverse hessian by determinantal averaging.  In _Proceedings of the 33rd International Conference on Neural Information Processing Systems_ , Red Hook, NY, USA, 2019. Curran Associates Inc.
  * Dua and Graff (2017) Dheeru Dua and Casey Graff.  UCI machine learning repository, 2017.  URL http://archive.ics.uci.edu/ml.
  * Elgabli et al. (2022) Anis Elgabli, Chaouki Ben Issaid, Amrit Singh Bedi, Ketan Rajawat, Mehdi Bennis, and Vaneet Aggarwal.  Fednew: A communication-efficient and privacy-preserving newton-type method for federated learning.  In _International Conference on Machine Learning_ , pages 5861–5877. PMLR, 2022.
  * Erdogdu and Montanari (2015) Murat A Erdogdu and Andrea Montanari.  Convergence rates of sub-sampled newton methods.  _Advances in Neural Information Processing Systems_ , 28, 2015.
  * Fernandes et al. (2015) Kelwin Fernandes, Pedro Vinagre, and Paulo Cortez.  A proactive intelligent decision support system for predicting the popularity of online news.  In _Portuguese conference on artificial intelligence_. Springer, 2015\.
  * Gupta et al. (2021) Vipul Gupta, Avishek Ghosh, Michał Dereziński, Rajiv Khanna, Kannan Ramchandran, and Michael W Mahoney.  LocalNewton: Reducing communication rounds for distributed learning.  In _Uncertainty in Artificial Intelligence_. PMLR, 2021.
  * Hard et al. (2018) Andrew Hard, Kanishka Rao, Rajiv Mathews, Swaroop Ramaswamy, Françoise Beaufays, Sean Augenstein, Hubert Eichner, Chloé Kiddon, and Daniel Ramage.  Federated learning for mobile keyboard prediction.  _arXiv preprint arXiv:1811.03604_ , 2018.
  * Huang et al. (2020) Li Huang, Yifeng Yin, Zeng Fu, Shifa Zhang, Hao Deng, and Dianbo Liu.  Loadaboost: Loss-based adaboost federated machine learning with reduced computational complexity on iid and non-iid intensive care data.  _PLOS ONE_ , 15(4):1–16, 04 2020.
  * Islamov et al. (2021) Rustem Islamov, Xun Qian, and Peter Richtárik.  Distributed second order methods with fast rates and compressed communication.  In _Proceedings of the 38th International Conference on Machine Learning_. PMLR, 2021.
  * Jeong et al. (2018) Eunjeong Jeong, Seungeun Oh, Hyesung Kim, Jihong Park, Mehdi Bennis, and Seong-Lyun Kim.  Communication-efficient on-device machine learning: Federated distillation and augmentation under non-iid private data.  _arXiv preprint arXiv:1811.11479_ , 2018.
  * Kairouz et al. (2021) Peter Kairouz, H. Brendan McMahan, et al.  Advances and open problems in federated learning.  _Foundations and Trends® in Machine Learning_ , 14(1–2):1–210, 2021.
  * Karimireddy et al. (2020) Sai Praneeth Karimireddy, Satyen Kale, Mehryar Mohri, Sashank Reddi, Sebastian Stich, and Ananda Theertha Suresh.  SCAFFOLD: Stochastic controlled averaging for federated learning.  In _Proceedings of the 37th International Conference on Machine Learning_. PMLR, 2020.
  * Li et al. (2020a) Tian Li, Anit Kumar Sahu, Ameet Talwalkar, and Virginia Smith.  Federated learning: Challenges, methods, and future directions.  _IEEE Signal Processing Magazine_ , 37(3):50–60, 2020a.
  * Li et al. (2020b) Tian Li, Anit Kumar Sahu, Manzil Zaheer, Maziar Sanjabi, Ameet Talwalkar, and Virginia Smith.  Federated optimization in heterogeneous networks.  In I. Dhillon, D. Papailiopoulos, and V. Sze, editors, _Proceedings of Machine Learning and Systems_ , 2020b.
  * Liu and Nocedal (1989) Dong C Liu and Jorge Nocedal.  On the limited memory BFGS method for large scale optimization.  _Mathematical programming_ , 45(1):503–528, 1989\.
  * Liu et al. (2020) Yaqiong Liu, Mugen Peng, Guochu Shou, Yudong Chen, and Siyu Chen.  Toward Edge Intelligence: Multiaccess Edge Computing for 5G and Internet of Things.  _IEEE Internet of Things Journal_ , 7, 2020.
  * Liu et al. (2021) Yi Liu, Yuanshao Zhu, and JQ James.  Resource-constrained federated learning with heterogeneous data: Formulation and analysis.  _IEEE Transactions on Network Science and Engineering_ , 2021.
  * Lyapunov (1992) Aleksandr Mikhailovich Lyapunov.  The general problem of the stability of motion.  _International journal of control_ , 55(3):531–534, 1992.
  * McMahan et al. (2017) H. B. McMahan, Eider Moore, D. Ramage, S. Hampson, and B. A. Y. Arcas.  Communication-efficient learning of deep networks from decentralized data.  In _AISTATS_ , 2017.
  * Nguyen et al. (2021) Tuan Dung Nguyen, Amir R. Balef, Canh T. Dinh, Nguyen H. Tran, Duy T. Ngo, Tuan Anh Le, and Phuong L. Vo.  Accelerating federated edge learning.  _IEEE Communications Letters_ , 25(10):3282–3286, 2021.
  * Pase et al. (2021) Francesco Pase, Marco Giordani, and Michele Zorzi.  On the convergence time of federated learning over wireless networks under imperfect CSI.  _IEEE International Conference on Communications Workshops (ICC WKSHPS)_ , pages 1–6, 2021.
  * Pham et al. (2020) Quoc-Viet Pham, Fang Fang, Vu Nguyen Ha, Md. Jalil Piran, Mai Les, Long Bao Le, Won-Joo Hwang, and Zhiguo Ding.  A Survey of Multi-Access Edge Computing in 5G and Beyond: Fundamentals, Technology Integration, and State-of-the-Art.  _IEEE Access_ , 8, 2020.
  * Qian et al. (2021) Xun Qian, Rustem Islamov, Mher Safaryan, and Peter Richtárik.  Basis matters: Better communication-efficient second order methods for federated learning.  _arXiv preprint arXiv:2111.01847_ , 2021.
  * Reddi et al. (2016) Sashank J Reddi, Jakub Konečnỳ, Peter Richtárik, Barnabás Póczós, and Alex Smola.  Aide: Fast and communication efficient distributed optimization.  _arXiv preprint arXiv:1608.06879_ , 2016.
  * Rieke et al. (2020) Nicola Rieke, Jonny Hancox, Wenqi Li, Fausto Milletari, Holger R Roth, Shadi Albarqouni, Spyridon Bakas, Mathieu N Galtier, Bennett A Landman, Klaus Maier-Hein, et al.  The future of digital health with federated learning.  _NPJ digital medicine_ , 3(1):1–7, 2020.
  * Safaryan et al. (2021) Mher Safaryan, Rustem Islamov, Xun Qian, and Peter Richtárik.  FedNL: Making Newton-type methods applicable to federated learning.  _arXiv preprint arXiv:2106.02969_ , 2021.
  * Shamir et al. (2014) Ohad Shamir, Nati Srebro, and Tong Zhang.  Communication-efficient distributed optimization using an approximate newton-type method.  In _Proceedings of the 31st International Conference on Machine Learning_ , pages 1000–1008, 22–24 Jun 2014.
  * Sheikh et al. (2020) Ruksar Sheikh, Mayank Patel, and Amit Sinhal.  Recognizing mnist handwritten data set using pca and lda.  In Garima Mathur, Harish Sharma, Mahesh Bundele, Nilanjan Dey, and Marcin Paprzycki, editors, _International Conference on Artificial Intelligence: Advances and Applications 2019_ , pages 169–177, Singapore, 2020\. Springer Singapore.
  * Shi et al. (2020) Yuanming Shi, Kai Yang, Tao Jiang, Jun Zhang, and Khaled B. Letaief.  Communication-efficient edge AI: algorithms and systems.  _IEEE Communications Surveys Tutorials_ , 22(4):2167–2191, 2020.