# Networked Federated Learning

Yasmin SarcheshmehPour, Yu Tian, Linli Zhang and Alexander Jung

###### Abstract

We develop the theory and algorithmic toolbox for networked federated learning
in decentralized collections of local datasets with an intrinsic network
structure. This network structure arises from domain-specific notions of
similarity between local datasets. Different notions of similarity are induced
by spatio-temporal proximity, statistical dependencies or functional
relations. Our main conceptual contribution is to formulate networked
federated learning using a generalized total variation minimization. This
formulation unifies and considerably extends existing federated multi-task
learning methods. It is highly flexible and can be combined with a broad range
of parametric models including Lasso or deep neural networks. Our main
algorithmic contribution is a novel networked federated learning algorithm
which is well suited for distributed computing environments such as edge
computing over wireless networks. This algorithm is robust against inexact
computations arising from limited computational resources including processing
time or bandwidth. For local models resulting in convex problems, we derive
precise conditions on the local models and their network structure such that
our algorithm learns nearly optimal local models. Our analysis reveals an
interesting interplay between the convex geometry of local models and the
(cluster-) geometry of their network structure.

###### Index Terms:

federated learning, clustering, complex networks, total variation,
regularization

##  I Introduction

Many important application domains generate collections of local datasets that
are related via an intrinsic network structure [1]. Two timely application
domains generating such networked data are (i) the high-precision management
of pandemics and (ii) the Internet of Things (IoT) [2]. Such local datasets
are generated by smartphones, wearables or industrial IoT devices [3]. These
local datasets are related via physical contact networks, social networks [4],
co-morbidity networks [5], or communication networks [6].

Federated learning (FL) is an umbrella term for machine learning (ML)
techniques that collaboratively train models on decentralized collections of
local datasets [7, 8, 9]. These methods carry out computations such as
gradient descent steps during model training at the location of data
generation, rather than first collecting all data at a central location [10].
FL methods are appealing for applications involving sensitive data (such as
healthcare) as they do not require the exchange of raw data but only model
(parameter) updates without leaking sensitive information in local datasets
[8]. Moreover, FL methods can offer robustness against malicious data
perturbation due to its intrinsic averaging or aggregation over large
collections of (mostly benign) datasets [11].

FL applications often face local datasets with different statistical
properties [12]. Each local dataset induces a separate learning task that
consists of learning (or optimizing) the parameters of a local model. This
paper studies an optimization method to train local models that are tailored
(or “personalized”) to the statistical properties of the corresponding local
dataset. This method is an instance of regularized empirical risk minimization
(or structural risk minimization). In particular, it uses a measure (see
Section III) for the variation of local model parameters as regularizer term.
We solve the resulting optimization (or learning) problem using a primal-dual
method that can be implemented as message passing over the network structure
of local datasets (see Section IV).

Clustered FL addresses the heterogeneity of local datasets using various forms
of a clustering assumption [13, 12, 14]. Informally, our clustering assumption
requires local datasets and their associated learning tasks to form a few
disjoint subsets or clusters. As a result, local datasets belonging to the
same cluster have similar statistical properties and, in turn, similar optimal
parameter values for the corresponding local models. Section II makes this
clustering assumption precise via Assumption 1. The main contribution of this
paper is a detailed characterization of the cluster structure and local model
geometry for local datasets that allow our methods to pool local datasets that
form statistically homogeneous clusters of datasets (see Section V).

What sets our approach apart from existing methods for clustered FL [12, 13]
is that we exploit known pairwise similarities between local datasets. These
similarities are encoded by the weighted undirected edges of an _empirical
graph_ [14]. Instead of a trivial combination of clustering methods and
cluster-wise model training, our NFL method (see Algorithm 1) interweaves the
pooling of local datasets with model training. We use the connectivity of the
empirical graph to guide this pooling (see Section IV).

We can interpret NFL methods as enhanced graph clustering methods that combine
the graph structure with the statistical properties of local datasets to
cluster the nodes of the empirical graph [15]. Another related interpretation
of NFL is that of a special form of multi-task learning. Indeed, each cluster
of local datasets gives rise to a separate learning task (finding optimal
parameters for all cluster nodes) [13]. Yet another interpretation of NFL,
that will be revealed in Section III, is that of network (flow) optimization
[16, 17].

Let us emphasize that our approach requires a useful choice for the empirical
graph of networked data. The connectivity of nodes (that represent local
datasets) in the empirical graph must reflect the clustering of local datasets
sharing similar statistical properties. Some application domains offer a
natural choice for the empirical graph, e.g., based on physical models,
functional relations, or the computing infrastructure [4].

If a natural choice for the empirical graph is not available, we can use
principled statistical tests for the similarity between two datasets [18]. We
demonstrate some of these methods in the numerical experiments of Section VI.
However, graph learning is beyond the scope of this paper (see Section VII).

###  I-A Related Work

Similar to [19, 20, 9, 21, 22], we use regularized empirical risk minimization
(RERM) to learn tailored models for local datasets. For each local dataset we
obtain a separate learning task that amounts to finding an (approximately)
optimal choice for the parameters of a local model. These individual learning
tasks are coupled via the undirected weighted edges of an empirical graph (see
Section II). In contrast to [22], which uses a probabilistic model for the
empirical graph, we consider the empirical graph as fixed and known (non-
random).

To capture intrinsic cluster structure of networked data, we use a generalized
total variation (GTV) of the model parameters as the regularizer. GTV unifies
and extends several existing notions of total variation [19, 20, 9, 21]. GTV
is parametrized by a penalty function which is used to measure the difference
of local model parameters at neighbouring nodes in the empirical graph.
Computationally, the main restriction for the choice of penalty function is
that it must allow for efficient computation of the corresponding proximity
operator (3). Some authors refer to such functions as “proximable” [23].

Our analysis reveals conditions on the network structure between local
datasets and their local models such that GTV minimization is able to identify
the cluster structure of the empirical graph. This is relevant for the
application of GTV minimization to clustered FL [13, 12]. In contrast to
existing work on clustered FL, we exploit a known similarity structure between
local datasets. We represent these similarities by the edges of an empirical
graph.

GTV minimization unifies and considerably extends well-known optimization
models for FL including the TV minimization and network Lasso (nLasso) [19,
20, 9, 21]. GTV is an instance of the non-quadratic regularizer put forward in
[21]. In contrast to [21], which uses a combination of gradient descent and
distributed averaging methods, we use a primal-dual method to solve the
resulting GTV minimization. We provide precise conditions such that GTV
minimization captures the inherent cluster structure of networked data.

The methods in [20] are limited to generalized linear models. In contrast, we
consider GTV minimization methods that can be combined with a wide range of
(potentially non-linear) parametrized models including graphical Lasso or deep
neural networks [24, 25].

A substantial body of work studies computational aspects of GTV minimization
[13]. Efficient algorithms for GTV minimization have been proposed for
relevant computational infrastructures such as wireless networks of low-
complexity devices [26, 27]. We would like to highlight the recent study [28]
of intrinsic computational complexity and efficient primal-dual methods for
GTV minimization. While [28] only uses the network diameter, our analysis
involves more fine-grained properties of empirical graph.

Similar to [28], our algorithmic approach to FL is based on an established
primal-dual method for optimization [29]. This primal-dual method is appealing
for FL application in several aspects. First, as we show in Section IV, the
primal-dual method [29, Alg. 6] can be implemented as a message passing
protocol over the empirical graph. Message passing algorithms are scalable to
massive collections of local datasets as long as their empirical graph is
sparse (e.g., a bounded degree network) [30].

The primal-dual method[29, Alg. 6] also offers robustness against limited
computational resources and imperfections. One example for such imperfections
are inexact evaluations of proximity operators. Evaluating the proximity
operator (3) is to solve an optimization problem, e.g., using iterative
methods such as gradient descent [23, 31]. Indeed, each iteration of the
primal-dual method requires to update local model parameters by solving a
separate model training problem for each local dataset. Given finite
computational resources, these local model updates can only be solved up to
some non-zero optimization error. The effect of inexact updates in primal-dual
methods has been analyzed recently [32].

In contrast to its computational aspects, the statistical aspects of GTV
minimization are far less understood [19, 9]. Our main focus is on the
dependency of solutions to GTV minimization on the cluster structure of
networked data. It is possible to frame GTV minimization as the learning or
recovery of group-sparse models [33]. However, it is unclear how the group-
sparse models underlying GTV minimization are related to the fine-grained
properties (such as cluster structure) of the empirical graph.

The closest to our work is the recent analysis of convex clustering [34]. As
convex clustering is a special case of GTV minimization (see Section III-C),
our work significantly extends the results in [34]. Moreover, in contrast to
[34], we characterize the cluster structure of local datasets using network
flows.

###  I-B Contribution

We next enumerate the main contributions of this paper.

  * •

We propose and study GTV minimization as an optimization model for NFL. GTV
minimization is an instance of RERM using the variation of local model
parameters as a regularizer. GTV minimization unifies and extends existing
optimization models for FL, including nLasso and MOCHA [9]. It is flexible as
it can be combined with a range of parametrized local models. The main
restriction is that local models must be parametrized by a common finite-
dimensional Euclidean space $\mathbb{R}^{d}$.

  * •

We explore the duality between GTV minimization and network flow optimization
[16]. This duality is a generalization of some well-known duality results for
network optimization [17]. The dual of GTV minimization generalizes the
minimum cost flow problem of [16] to vector-valued flows.

  * •

Our main methodological contribution is a novel family of NFL algorithms (see
Algorithm 1). This family is obtained by solving different instances of GTV
minimization with an established primal-dual method for structured
optimization [23]. Algorithm 1 alternates between separate local model
parameter updates that are regularized by previously received model parameters
and sharing these updates between neighbouring nodes (see Section IV). The
computational complexity of the algorithm is essentially the same as the
computational complexity of training the local models separately Algorithm 1
can be combined with a wide range of parametric ML models and notions of TV
(obtained for different GTV penalty functions). The main requirement is merely
the existence of an efficient training method for the local models. The choice
for the GTV penalty function is only restricted by requiring an efficient way
to compute its convex conjugate.

  * •

Our main analytical contribution is an upper bound on the estimation error
incurred by GTV minimization. This upper bound reveals sufficient conditions
on the local datasets and network structure such that our method achieves the
performance of an oracle-based method that perfectly knows the true cluster
structure of the data network. We emphasize that our analysis only applies to
GTV minimization (i) using a penalty function being a norm (ii) and local
models characterized by smooth and convex loss functions. Thus, our bounds do
not apply to methods that either use graph Laplacian quadratic form as
regularizer (such as MOCHA) or local models resulting in non-convex loss
functions (deep nets). However, our analysis might provide also insight into
the statistical properties of GTV minimization for non-convex local models
(i.e., resulting in non-convex loss functions).

###  I-C Outline

Section II introduces the concept of an empirical graph to represent
collections of local datasets, the corresponding local models as well as their
similarity structure. Section III introduces GTV as a measure for the
variation of local model parameters across the edges in the empirical graph.
As discussed in Section III-A, GTV minimization balances the variation of
local model parameters over well-connected local datasets (forming a cluster)
and incurring a small loss (training error) for each local dataset. The dual
problem to GTV minimization is then explained in Section III-B. Section III-C
presents several useful interpretations of GTV minimization and its dual.
Section IV applies a well-known primal-dual optimization method to solve GTV
minimization and its dual in a fully distributed fashion via message passing
over the empirical graph (see Algorithm 1). The results of numerical
experiments are discussed in Section VI.

###  I-D Notation

The identity matrix of size $n\\!\times\\!n$ is denoted $\mathbf{I}_{n}$, with
the subscript omitted if the size $n$ is clear from context. We use
$\|\cdot\|$ to denote some norm defined on the Euclidean space
$\mathbb{R}^{d}$ and $\|\cdot\|_{*}$ to denote its dual norm [35, Appx. 1.6.].
Two important examples are the Euclidean norm $\|{\bf
w}\|_{2}\\!:=\\!\sqrt{\sum_{j=1}^{d}w_{j}^{2}}$ and the $\ell_{1}$ norm
$\|{\bf w}\|_{1}\\!:=\\!\sum_{j=1}^{d}|w_{j}|$ of a vector ${\bf
w}\\!=\\!(w_{1},\ldots,w_{d})^{T}\in\mathbb{R}^{d}$ . It will be convenient to
use the notation $(1/2\tau)$ instead of $(1/(2\tau))$. We will need the
(vector-wise) clipping operator

$\vspace{-2mm}\mathcal{T}^{(\gamma)}({\bf w})\\!:=\\!\begin{cases}\gamma{\bf
w}/\|{\bf w}\|_{2}&\mbox{ for }\|{\bf w}\|_{2}\\!\geq\\!\gamma\\\ {\bf
w}&\mbox{ otherwise.}\end{cases}$ (1)

The scalar clipping operator $\mathcal{T}^{(\gamma)}(w)$ is obtained as a
special case of (1) by considering the scalar $w$ as a vector with a single
entry (where $\left\|{\mathbf{w}}\right\|_{2}=|w|$). Given a closed proper
convex function $f({\bf x})$ with domain being a subset of $\mathbb{R}^{d}$,
we define its associated convex conjugate function as [35]

$f^{*}({\bf x}):=\sup\limits_{{\bf z}\in\mathbb{R}^{d}}{\bf x}^{T}{\bf
z}-f({\bf z}).\vspace{-2mm}$ (2)

We will use the proximity operator of a closed proper convex function $f({\bf
x})$, defined as [26]

${\rm\bf prox}_{f,\rho}({\bf x})\\!:=\\!\operatorname*{argmin}_{{\bf
x}^{\prime}}f({\bf x}^{\prime})\\!+\\!(\rho/2)\|{\bf x}-{\bf
x}^{\prime}\|^{2}_{2}\mbox{ with }\rho\\!>\\!0.$ (3)

Note that the minimum in (3) exists and is unique since the objective function
is strongly convex [35].

##  II Problem Formulation

We find it useful to represent networked data by an undirected weighted
_empirical graph_ $\mathcal{G}=(\mathcal{V},\mathcal{E})$. For notational
convenience, we identify the nodes of an empirical graph with natural numbers,
$\mathcal{V}=\\{1,\ldots,n\\}$. Each node $i\\!\in\\!\mathcal{V}$ of the
empirical graph $\mathcal{G}$ carries a separate local dataset
$\mathcal{X}^{(i)}$. It might be instructive to think of a local dataset
$\mathcal{X}^{(i)}$ as a labeled dataset

$\mathcal{X}^{(i)}:=\left\\{\big{(}{\bf
x}^{(i,1)},y^{(i,1)}\big{)},\ldots,\big{(}{\bf
x}^{(i,m_{i})},y^{(i,m_{i})}\big{)}\right\\}.\vspace{-2mm}$ (4)

Here, ${\bf x}^{(r)}$ and $y^{(r)}$ denote, respectively, the feature vector
and true label of the $r$-th data point in the local dataset
$\mathcal{X}^{(i)}$. Note that the size $m_{i}$ of the local dataset might
vary across nodes $i\in\mathcal{V}$. Figure 1 depicts an empirical graph with
$n\\!=\\!11$ nodes $\mathcal{V}=\\{1,\ldots,n\\}$, each carrying a local
dataset $\mathcal{X}^{(i)}$.

We highlight that our NFL method (see Section IV) is not restricted to local
datasets of the form (4). Indeed, Algorithm 1 and its analysis (see Section V)
only requires indirect access to $\mathcal{X}^{(i)}$ via the evaluation of
some local loss function $L_{i}\left({\bf v}\right)$. The value
$L_{i}\left({\bf v}\right)$ measures how well a model with parameters ${\bf
v}$ fits the local dataset $\mathcal{X}^{(i)}$ (see Section II-A). We will
study different choices for the local loss function in Section VI.

Let us point out two particular aspects of our data-access model via the
evaluation of local loss functions. First, it lends naturally to privacy-
friendly methods as they do not need to share raw data $\mathcal{X}^{(i)}$.
Instead, our methods only exchange (local) information about local loss
function $L_{i}\left(\cdot\right)$ such as the gradient $\nabla
L_{i}\left({\bf v}\right)$ or the proximity operator value ${\rm\bf
prox}_{L_{,}\left(\rho\right)}({\bf v})$ (3) for a given choice
$\mathbf{w}^{(i)}={\bf v}$ for the local model parameter vector. This
information is typically obtained from averages over data points and therefore
revealing only little information about individual data points (if the sample
size is no too small).

Beside its privacy-friendliness, our data access model also handles
applications where only a fraction of local datasets are accessible. This is
relevant for wireless sensor networks that consist of battery powered devices
for computation and wireless communication [36]. The lack of access to the
local dataset at some node $i$ can be taken into account by using a trivial
loss function $L_{i}\left({\bf v}\right)=0$ for all parameter vectors ${\bf
v}\in\mathbb{R}^{d}$ (see Section VI).

$\mathcal{X}^{(i)},{\bf w}^{(i)}$$\mathcal{X}^{(j)},{\bf
w}^{(j)}$$\overline{\mathbf{w}}^{(\mathcal{C}_{1})}$$\overline{\mathbf{w}}^{(\mathcal{C}_{2})}$$\overline{\mathbf{w}}^{(\mathcal{C}_{3})}$$A_{i,j}$
Figure 1:  We represent networked data and corresponding models using an
undirected empirical graph
$\mathcal{G}=\big{(}\mathcal{V},\mathcal{E}\big{)}$. Each node
$i\in\mathcal{V}$ of the graph carries a local dataset $\mathcal{X}^{(i)}$ and
model weights $\mathbf{w}^{(i)}$ which are scored using a local loss function
$L_{i}\left(\mathbf{w}^{(i)}\right)$ (that encapsulates the local dataset
$\mathcal{X}^{(i)}$). Two nodes are connected by a weighted edge
$\\{i,i^{\prime}\\}$ if they carry datasets with similar statistical
properties. The amount of similarity is encoded in an edge weight
$A_{i,i^{\prime}}>0$ (indicated by thickness of the links). We rely on a
clustering assumption, requiring optimal parameter vectors for nodes in the
same cluster $\mathcal{C}^{(c)}\subseteq\mathcal{V}$ to be nearly identical.
The empirical graph is partitioned into three disjoint clusters
$\mathcal{C}^{(1)},\mathcal{C}^{(2)},\mathcal{C}^{(3)}$. Note that our FL
method does not require the (typically unknown) partition but rather learns
the partition based on the local datasets and network structure of
$\mathcal{G}$.

An undirected edge $\\{i,i^{\prime}\\}\\!\in\\!\mathcal{E}$ indicates that the
corresponding local datasets $\mathcal{X}^{(i)}$ and
$\mathcal{X}^{(i^{\prime})}$ have similar statistical properties. The strength
of the similarity is quantified by the edge weight $A_{i,i^{\prime}}\\!>\\!0$.
We will also use the notation $A_{e}:=A_{i,i^{\prime}}$ for an edge
$e=\\{i,i^{\prime}\\}$. It will be convenient to indicate absence of an edge
between two nodes $i,i^{\prime}\in\mathcal{V}$ by a zero weight, i.e.,
$A_{i,i^{\prime}}=0$ if and only if $\\{i,i^{\prime}\\}\notin\mathcal{E}$.

The undirected edge $\\{i,i^{\prime}\\}\in\mathcal{E}$ encodes a symmetric
notion of similarity between local datasets. If the local dataset at node $i$
is (statistically) similar to the local dataset at node $i^{\prime}$ then also
vice-versa. The symmetric nature of the similarities between local datasets is
also reflected in the edge weights,

$A_{i,i^{\prime}}=A_{i^{\prime},i}\mbox{ for any two nodes
}i,i^{\prime}\in\mathcal{V}.\vspace{-2mm}$

It will be convenient for the formulation and analysis of our NFL algorithm
(see Algorithm 1) to orient the edges in $\mathcal{E}$. In particular, we
define the head and tail of an undirected edge
$e=\\{i,i^{\prime}\\}\in\mathcal{E}$ as $e_{+}:=\min\\{i,i^{\prime}\\}$ and
$e_{-}:=\max\\{i,i^{\prime}\\}$, respectively. The entire set of directed
edges for an empirical graph is obtained as

$\overrightarrow{\mathcal{E}}:=\big{\\{}\left(i,i^{\prime}\right):i,i^{\prime}\in\mathcal{V},i<i^{\prime}\mbox{
and }\\{i,i^{\prime}\\}\in\mathcal{E}\big{\\}}.$ (5)

We abuse notation and use $\mathcal{E}$ not only to denote the set of
undirected edges but also to denote the set (5) of directed edges in the
empirical graph $\mathcal{G}$.

There are two vector spaces that are naturally associated with an empirical
graph $\mathcal{G}$. The “node space” $\mathcal{W}$ consists of maps
$\mathbf{w}:\mathcal{V}\rightarrow\mathbb{R}^{d}:i\mapsto\mathbf{w}^{(i)}$
that assign each node $i\in\mathcal{V}$ a vector
$\mathbf{w}^{(i)}\in\mathbb{R}^{d}$. The “edge space” $\mathcal{U}$ of all
maps ${\bf u}:\mathcal{E}\rightarrow\mathbb{R}^{d}:e\mapsto\mathbf{u}^{(e)}$
that assign each edge $e\in\mathcal{E}$ a vector
$\mathbf{u}^{(e)}\in\mathbb{R}^{d}$. These two spaces are linked via the
block-incidence matrix $\mathbf{D}$ with entries $D_{e,i}=1$ for $i=e_{+}$,
$D_{e,i}=-1$ for $i=e_{-}$, and $D_{e,i}=0$ otherwise. The block-incidence
matrix $\mathbf{D}$ represents a linear map

${\bf D}:\mathcal{W}\rightarrow\mathcal{U}:\mathbf{w}\mapsto\mathbf{u}\mbox{
with }\mathbf{u}^{(e)}=\mathbf{w}^{(e_{+})}-\mathbf{w}^{(e_{-})}\vspace{-2mm}$
(6)

with the adjoint (transpose) $\mathbf{D}^{T}$ representing another linear map,

$\mathbf{D}^{T}:\mathcal{U}\rightarrow\mathcal{W}:\mathbf{u}\mapsto\mathbf{w}\mbox{
with
}\mathbf{w}^{(i)}=\sum_{e\in\mathcal{E}}\sum_{i=e_{+}}\mathbf{u}^{(e)}-\sum_{i=e_{-}}\mathbf{u}^{(e)}.\vspace{-4mm}$
(7)

###  II-A Networked Models

A networked model consists of a separate local model for each local dataset
$\mathcal{X}^{(i)}$. Our approach to NFL allows for a large variety of design
choices for the local models. We only require all local models to be
parametrized by a common finite-dimensional Euclidean space $\mathbb{R}^{d}$.
This setting covers some widely-used ML models such as (regularized)
generalized linear models or linear time series models [20, 37]. However, our
setting does not cover non-parametric local models such as decision trees.

Networked models are parametrized by a map $\mathbf{w}\in\mathcal{W}$ that
assigns each node $i\in\mathcal{V}$ in the empirical graph $\mathcal{G}$ a
local model parameter vector $\mathbf{w}^{(i)}\in\mathbb{R}^{d}$, 111With a
slight abuse of notation we will refer by $\mathbf{w}^{(i)}$ also to the
entire collection of local model parameters. ${\bf
w}:\mathcal{V}\rightarrow\mathbb{R}^{d}:i\mapsto{\bf w}^{(i)}.$ We measure the
usefulness of a particular choice for the local model parameters
$\mathbf{w}^{(i)}$ by some local loss function
$L_{i}\left(\mathbf{w}^{(i)}\right)$. Unless stated otherwise, we consider
local loss functions that are convex and differentiable. The NFL method
proposed in Section IV allows for different choices for the local loss
functions. These different choices might be obtained, in turn, from different
combinations of ML models and performance metrics [38, Ch. 3].

From a computational perspective, our main requirement on the choice for the
local loss function $L_{i}\left(\mathbf{w}^{(i)}\right)$ is that it allows for
efficient solution of the regularized problem,

$\min_{{\bf w}^{\prime}\in\mathbb{R}^{d}}L_{i}\left({\bf
w}^{\prime}\right)+\lambda\|{\bf w}^{\prime}-{\bf
w}^{\prime\prime}\|^{2}_{2}.\vspace{-2mm}$ (8)

The computational complexity of our approach depends on the ability to
efficiently solve (8) for any given $\lambda\in\mathbb{R}_{+}$ and ${\bf
w}^{\prime\prime}\in\mathbb{R}^{d}$. Note that solving (8) is equivalent to
evaluating the proximity operator ${\rm\bf
prox}_{L_{i}\left(\cdot\right),2\lambda}({\bf w}^{\prime\prime})$.

The NFL method in Section IV applies to parametric models that can be trained
by minimizing a loss function $L_{i}\left(\cdot\right)$ whose proximity
operator can be evaluated efficiently. Convex functions for which the
proximity operator can be computed efficiently are sometimes referred to as
“proximable” or “simple” [23]. Note that the shape of the loss function
typically depends on both, the choice for the local model and the metric used
to measure prediction errors [38, Ch. 4].

Our focus is on applications where the local loss functions
$L_{i}\left(\mathbf{w}^{(i)}\right)$ do not carry sufficient statistical power
to guide learning of model parameters $\mathbf{w}^{(i)}$. As a case in point,
consider a local dataset $\mathcal{X}^{(i)}$ of the form (4), with feature
vectors ${\bf x}^{(r)}\in\mathbb{R}^{d}$ with $m_{i}\ll d$. We would like to
learn the parameter vector $\mathbf{w}^{(i)}$ of a linear hypothesis $h({\bf
x})=\big{(}\mathbf{w}^{(i)}\big{)}^{T}{\bf x}$. Linear regression methods
learn the parameter vector by minimizing the average squared error loss
$L_{i}\left(\mathbf{w}^{(i)}\right)=(1/m_{i})\sum_{r=1}^{m_{i}}\big{(}y^{(r)}-\big{(}\mathbf{w}^{(i)}\big{)}^{T}{\bf
x}^{(r)}\big{)}^{2}$. However, for $m_{i}\ll d$ (the “high-dimensional”
regime) the minimum of $L_{i}\left(\cdot\right)$ is not unique and might also
provide a poor hypothesis incurring large prediction errors on data points
outside $\mathcal{X}^{(i)}$ [38, Ch. 6]. Training linear models in the high-
dimensional regime requires regularization such as in ridge regression or
Lasso [39].

The main theme of this paper is to use the empirical graph $\mathcal{G}$ to
regularize the learning of local model parameters by requiring them to not
vary too much over edges with large weights. Section III introduces the
concept of GTV as a quantitative measure for the variation of local parameter
vectors. Regularization by requiring a small GTV is an instance of the
smoothness assumption used in semi-supervised learning [14].

Our analysis of GTV minimization in Section V relates its underlying
smoothness assumption to a clustering assumption. Section II-B formalizes this
clustering assumption which requires local model parameters to be constant
over subsets (clusters) of nodes in the empirical graph. Theorem 2 then offers
precise conditions on the empirical graph and local loss functions such that
GTV minimization successfully recovers the clusters of nodes.

###  II-B Clustering Assumption

Consider networked data with empirical graph
$\mathcal{G}=\left(\mathcal{V},\mathcal{E}\right)$. Each node $i$ in the graph
carries a local dataset $\mathcal{X}^{(i)}$ and a local model with parameters
$\mathbf{w}^{(i)}$. Our goal is to learn the local model parameters
$\mathbf{w}^{(i)}$ for each node $i\in\mathcal{V}$. The key assumption of
clustered FL is that the local datasets form clusters with local datasets in
the same cluster having similar statistical properties [13]. Given a cluster
$\mathcal{C}$ of nodes, it seems natural to pool their local datasets or,
equivalently, add their local functions to learn a cluster-specific parameter
vector

$\overline{{\bf w}}^{(\mathcal{C})}=\operatorname*{argmin}_{{\bf
v}\in\mathbb{R}^{d}}f^{(\mathcal{C})}\left({\bf v}\right)\mbox{ with
}f^{(\mathcal{C})}\left({\bf v}\right):=\sum_{i\in\mathcal{C}}L_{i}\left({\bf
v}\right).\vspace{-2mm}$ (9)

Note that (9) cannot be implemented in practice since we typically do not know
the cluster $\mathcal{C}$. The main analytical contribution of this paper is
an upper bound for the deviation between solutions of GTV minimization and the
cluster-wise (but impractical) learning problem (9). This bound characterizes
the statistical properties of FL algorithms that are obtained by applying
optimization techniques for solving GTV minimization (see Section IV).

The solution $\overline{{\bf w}}^{(\mathcal{C})}$ of (9) minimizes the
aggregation (sum) of all local loss functions that belong to the same cluster
$\mathcal{C}\subseteq\mathcal{V}$. Thus, $\overline{{\bf w}}^{(\mathcal{C})}$
is the optimal model parameter for a training set obtained by pooling all
local datasets that belong to the cluster $\mathcal{C}$. As indicated by our
notation, we tacitly assume that the solution to (9) is unique. The uniqueness
of the solution in (9) will be ensured by Assumption 2 below.

We now make our assumption of datasets in the same cluster “having similar
statistical properties” precise. In particular, we require the local loss
functions at nodes $i\in\mathcal{C}$ in the same cluster $\mathcal{C}$ to have
nearby minimizers. Thus, we require a small deviation $\left\|{{\bf
v}^{(i)}-\overline{{\bf w}}^{(\mathcal{C})}}\right\|$ between the minimizer
${\bf v}^{(i)}$ of $L_{i}\left(\cdot\right)$ and the corresponding cluster-
wise optimal parameter vector (9). This requirement is, for differentiable and
convex loss functions, equivalent to requiring a small gradient of the local
loss functions at the cluster-wise minimizer (9). It will be convenient for
our analysis to formulate this requirement by upper bounding the dual norm
$\left\|{\nabla L_{i}\left(\overline{{\bf
w}}^{(\mathcal{C})}\right)}\right\|_{*}$ of the local loss gradient.

######  Assumption 1 (Clustering).

Consider some networked data represented by an empirical graph $\mathcal{G}$
whose nodes carry local loss functions $L_{i}\left({\bf v}\right)$, for
$i\in\mathcal{V}$. There is a partition of the nodes $\mathcal{V}$ into
disjoint clusters

$\displaystyle{\mathcal{P}\\!=\\!\\{\mathcal{C}^{(1)},\ldots,\mathcal{C}^{(k)}\\}\mbox{
with }\mathcal{C}^{(c)}\cap\mathcal{C}^{(c^{\prime})}=\emptyset,}$ (10)
$\displaystyle{\quad\mbox{ for }c\\!\neq\\!c^{\prime}\mbox{ and
}\mathcal{V}=\mathcal{C}^{(1)}\cup\ldots\cup\mathcal{C}^{(k)}.}$

Moreover, for each cluster $\mathcal{C}^{(c)}\in\mathcal{P}$,

$\left\|{\nabla L_{i}\left(\overline{{\bf
w}}^{(c)}\right)}\right\|_{*}\leq\delta^{(i)}\mbox{ for all
}i\in\mathcal{C}^{(c)}.\vspace{-2mm}$ (11)

Here, $\overline{{\bf w}}^{(c)}\\!\in\\!\mathbb{R}^{d}$ denotes the solution
of the cluster-wise minimization (9) for cluster $\mathcal{C}^{(c)}$.

The clustering assumption requires the dual norm $\left\|{\nabla L_{i}\left(\
\overline{{\bf w}}^{(c)}\right)}\right\|_{*}$ to be bounded by a constant
$\delta^{(i)}$ for each nodes $i\in\mathcal{V}$ in the empirical graph. We can
interpret this norm as a measure for the discrepancy between the cluster-wise
minimizer $\overline{{\bf w}}^{(c)}$ (see (9)) and the minimizers of the local
loss functions $L_{i}\left(\overline{{\bf w}}^{(c)}\right)$ for each node
$i\in\mathcal{C}^{(c)}$.

Section IV uses the cluster-wise minimization (9) as a theoretical device to
analyze the solutions of GTV minimization (16). It is important to note that
(9) does not inform a practical FL method as it requires knowledge of the
clusters in the partition (10). It might be unrealistic to assume perfect
knowledge of the partition (10) postulated by Assumption 1. Rather, we show
that GTV minimization is able to recover this partition using solely the edges
of the empirical graph $\mathcal{G}$.

Section III formulates NFL as GTV minimization which is an instance of RERM.
GTV minimization is enforces “clusteredness” of local model parameters by
requiring a small variation across edges in the empirical graph. Under
Assumption 1, this regularization strategy will be useful if many (large
weight) edges connect nodes in the same cluster but only few (small weight)
edges connect nodes in different clusters. Section V presents a precise
condition on the network structure such that GTV minimization succeeds in
capturing the true underlying cluster structure of the local loss functions.

The analysis of the NFL method proposed in Section V requires the local loss
functions to be convex and smooth. Moreover, we require their (partial) sums
in the cluster-wise objective $f^{(c)}\left(\cdot\right)$ (9) to be strongly
convex [40, Exercise 1.9].

######  Assumption 2 (Convexity and Smoothness).

For each node $i\in\mathcal{V}$, the local loss function
$L_{i}\left(\mathbf{w}^{(i)}\right)$ is convex and differentiable with
gradient satisfying

$\left\|{\nabla L_{i}\left({\bf v}^{\prime}\right)-\nabla L_{i}\left({\bf
v}\right)}\right\|_{*}\leq\beta^{(i)}\left\|{{\bf v}^{\prime}-{\bf
v}}\right\|.\vspace{-2mm}$ (12)

For each cluster $\mathcal{C}^{(c)}\in\mathcal{P}$ in the partition (10), the
cluster-wise objective $f^{(c)}\left(\cdot\right)$ (9) is strongly convex,

$\displaystyle f^{(c)}\left({\bf v}^{\prime}\right)\geq f^{(c)}\left({\bf
v}\right)+$ $\displaystyle{\big{(}{\bf v}^{\prime}\\!-\\!{\bf
v}\big{)}^{T}\partial f^{(c)}\left({\bf v}\right)+(\alpha^{(c)}/2)\|{{\bf
v}^{\prime}\\!-\\!{\bf v}}\|^{2},}$ (13) $\displaystyle{\mbox{ for any }{\bf
v}^{\prime},{\bf v}\in\mathbb{R}^{d}.}$

Here, $\alpha^{(c)}>0$ is a positive constant that might be different for
different clusters $\mathcal{C}^{(c)}$. The norm $\left\|{\cdot}\right\|$ in
(12), (13) is the dual of the norm $\left\|{\cdot}\right\|_{*}$ used in (12)
and (11).

Assumption 2 is rather standard in FL literature [12, 21]. In particular,
Assumption 2 is satisfied by many important ML models [20]. Assumption 2 does
not hold for many deep learning models that result in non-convex loss
functions [7]. Nevertheless, we expect our theoretical analysis to provide
useful insight also for settings where Assumption 2 is violated.

We emphasize that Assumption 2 does not require strong convexity for each
local loss function $L_{i}\left(\cdot\right)$ individually. Rather, it only
requires their cluster-wise sums (9) to be strongly convex. We also allow for
trivial local loss functions that are constant and might represent
inaccessibility of local datasets due to privacy-constraints or lack of
computational resources.

The NFL method in Section IV can tolerate the presence of non-informative
local loss functions by exploiting the similarities between local datasets as
reflected by the edges in the empirical graph $\mathcal{G}$. Supplementary
material IX-B presents examples for local loss functions, obtained from some
widely-used ML models, that conform to Assumption 2.

##  III Generalized Total Variation Minimization

The clustering Assumption 1 suggests to learn the model parameters
$\mathbf{w}^{(i)}$ via the cluster-wise optimization (9). For each cluster
$\mathcal{C}^{(c)}$ in the partition (10), we use the solution of (9) as the
local model parameters at all nodes $i\in\mathcal{C}^{(c)}$. However, this
approach is not practical since the partition (10) is typically unknown and
therefore we cannot directly implement (9). Instead, we use the empirical
graph $\mathcal{G}$ to penalize variations of local model parameter
$\mathbf{w}^{(i)}$ over well-connected nodes (see Section III-A).

We hope that penalizing their variation (over the edges in the empirical
graph) favours local model parameters that are approximately constant over
nodes in the same cluster (10). For this approach to be successful, the nodes
in the same cluster must be densely connected by many edges (or large weight)
in the empirical graph, which should have only few edges (with low weight)
between nodes in different clusters. We will make this informal assumption
precise in Section V. For now, we use the informal clustering assumption to
motivate GTV as a useful regularizer for learning the local model parameters.

If the cluster structure of $\mathcal{G}$ is reflected by a high density of
edges within clusters and few boundary edges between them, it seems reasonable
to require a small variation of local parameter vectors $\mathbf{w}^{(i)}$
across edges. We measure the variation of local parameter vectors
$\mathbf{w}\in\mathcal{W}$ across the edges in $\mathcal{G}$ via the variation
$\mathbf{u}:e\in\mathcal{E}\mapsto\mathbf{u}^{(e)}:=\mathbf{w}^{(e_{+})}-\mathbf{w}^{(e_{-})}$.
Using the block-incidence matrix (6) we can express the variation of
$\mathbf{w}$ more compactly as $\mathbf{u}=\mathbf{D}\mathbf{w}$.

A quantitative measure for the variation of local model parameters
$\mathbf{w}$ is the GTV

$\displaystyle{\left\|{{\bf w}}\right\|_{\rm
GTV}:=\sum_{\\{i,i^{\prime}\\}\in\mathcal{E}}A_{i,i^{\prime}}\phi\big{(}\mathbf{w}^{(i^{\prime})}-\mathbf{w}^{(i)}\big{)}}$
(14)

with some convex penalty function
$\phi(\cdot):\mathbb{R}^{d}\rightarrow\mathbb{R}$. We also define the GTV for
a subset of edges $\mathcal{S}\subseteq\mathcal{E}$ as

$\|{\bf
w}\|_{\mathcal{S}}:=\sum_{\\{i,i^{\prime}\\}\in\mathcal{S}}A_{i,i^{\prime}}\phi\big{(}\mathbf{w}^{(i^{\prime})}-\mathbf{w}^{(i)}\big{)}.\vspace{-2mm}$
(15)

The GTV (14) provides a whole ensemble of variation measures. This ensemble is
parametrized by a penalty function $\phi({\bf v})\in\mathbb{R}$ which we
tacitly assume to be convex. The penalty function $\phi(\cdot)$ is an
important design choice that determines the computational and statistical
properties of the resulting GTV minimization problem (see Section IV and
Section V). Two popular choices are $\phi({\bf v}):=\|{\bf v}\|_{2}$, which is
used by nLasso [19], and $\phi({\bf v}):=(1/2)\|{\bf v}\|^{2}_{2}$ which is
used by the method MOCHA [9]. Another recent FL method for networked data uses
the choice $\phi({\bf v}):=\|{\bf v}\|_{1}$ [41].

Different choices for the penalty function offer different trade-offs between
computational complexity and statistical properties of the resulting FL
algorithms. As a case in point, the penalty
$\phi(\mathbf{u})=\|\mathbf{u}\|_{2}$ (used in nLasso) is computationally more
challenging than the penalty $\phi({\bf u})=(1/2)\|\mathbf{u}\|_{2}^{2}$ (used
in MOCHA [9]). On the other hand, nLasso is more accurate in learning models
for data with specific network structures (such as chains) that are
challenging for GTV minimization method using the smooth penalty $\phi({\bf
v}):=(1/2)\|{\bf v}\|^{2}_{2}$ [42].

Section IV designs FL methods whose main computational steps include the
computation of the proximity operator ${\rm\bf prox}_{\phi^{*},\rho}(\cdot)$
for the convex conjugate $\phi^{*}$ of the GTV penalty function $\phi(\cdot)$.
Thus, for these methods to be computationally tractable we must choose
$\phi(\cdot)$ such that the proximity operator ${\rm\bf
prox}_{\phi^{*},\rho}(\cdot)$ can be computed (evaluated) efficiently.222The
difficulty of computing the proximity operator ${\rm\bf
prox}_{\phi^{*},\rho}({\bf u})$ is essentially the same as that of computing
the proximity operator ${\rm\bf prox}_{\phi,\rho}({\bf u})$. Indeed, these two
proximity operators are related via the identity ${\bf u}={\rm\bf
prox}_{\phi,1}({\bf u})+{\rm\bf prox}_{\phi^{*},1}({\bf u})$ [43].
Supplementary material IX-C discusses some choices for the GTV penalty
function $\phi(\cdot)$ for which this is the case.

###  III-A The Primal Problem

GTV minimization learns the local model parameters $\mathbf{w}^{(i)}$ by
balancing (the sum of) local loss functions and GTV (14),

$\widehat{\mathbf{w}}\in\underset{\mathbf{w}\in\mathcal{W}}{\mathrm{arg\
min}}\sum_{i\in\mathcal{V}}L_{i}\left(\mathbf{w}^{(i)}\right)+\lambda\|\mathbf{w}\|_{\rm
GTV}.\vspace{-2mm}$ (16)

GTV minimization (16) is an instance of RERM, using the scaled GTV
$\lambda\|\mathbf{w}\|_{\rm GTV}$ as regularizer. The empirical risk incurred
by the local model parameters $\mathbf{w}\in\mathcal{W}$ is measured by the
sum of the local loss functions
$\sum_{i\in\mathcal{V}}L_{i}\left(\mathbf{w}^{(i)}\right)$. Section III
applies the primal-dual method [29, Alg. 6] to compute (approximate) solutions
of (16).

The regularization parameter $\lambda\\!>\\!0$ in (16) steers the preference
for learning parameter vectors $\mathbf{w}^{(i)}$ with small GTV versus
incurring small local loss
$\sum_{i\in\mathcal{V}}L_{i}\left(\mathbf{w}^{(i)}\right)$. The choice of
$\lambda$ can be guided by cross validation [39] or by our analysis of the
solutions of (16) in Section V.

Increasing the value of $\lambda$ results in the solutions of (16) becoming
increasingly clustered, with local model parameters
$\widehat{\mathbf{w}}^{(i)}$ being constant over (increasingly) large subsets
of nodes. Choosing $\lambda$ larger than some critical value, that depends on
the shape of the local loss functions and the network structure of
$\mathcal{G}$, results in $\widehat{\mathbf{w}}^{(i)}$ being constant over all
nodes $i\in\mathcal{V}$. Section V offers precise conditions on the local loss
functions and the empirical graph such that the solutions of (16) capture the
(unknown) underlying partition (10).

GTV minimization (16) unifies and considerably extends some well-known methods
for distributed optimization and learning. In particular, the nLasso [19] is
obtained from (16) for the choice $\phi({\bf v}):=\left\|{{\bf
v}}\right\|_{2}$. The MOCHA method [9] is obtained from (16) for the choice
$\phi({\bf v}):=(1/2)\left\|{{\bf v}}\right\|_{2}^{2}$. Another special case
of (16), obtained for the choice $\phi({\bf v}):=\left\|{{\bf
v}}\right\|_{1}$, has been studied recently [41].

###  III-B The Dual Problem

The solutions of GTV minimization (16) can be conveniently characterized and
computed by introducing another optimization problem that is dual (in a sense
that we make precise promptly) to (16). We obtain this dual problem by using
the convex conjugate $h^{*}$ (see (2)) of a convex function $h({\bf x})$ (see
[35]). The convex conjugate offers an alternative (or dual) representation of
a convex function $h({\bf x})$ via

$h({\bf x}):=\sup\limits_{{\bf z}}{\bf x}^{T}{\bf z}-h^{*}({\bf
z}).\vspace{-3mm}$ (17)

This alternative (or dual) representation of convex functions lends naturally
to a dual problem for GTV minimization (16) .

While the domain of GTV minimization (16) are the local model parameters
$\mathbf{w}^{(i)}\in\mathbb{R}^{d}$, for all nodes $i\in\mathcal{V}$, the
domain of the dual problem will be flow vectors
$\mathbf{u}^{(e)}\in\mathbb{R}^{d}$, for all edges $e\in\mathcal{E}$. To
formulate the dual problem of GTV minimization (16), we first rewrite it more
compactly as (see (6) and (14)),

$\displaystyle{\widehat{\mathbf{w}}\in\underset{{\mathbf{w}}\in\mathcal{W}}{\mathrm{arg\
min}}\ f(\mathbf{w})\\!+\\!g(\mathbf{D}\mathbf{w})}$ (18)
$\displaystyle{\mbox{with
}f(\mathbf{w}):=\sum_{i\in\mathcal{V}}L_{i}\left(\mathbf{w}^{(i)}\right)},$
$\displaystyle{\mbox{and
}g(\mathbf{u})\\!:=\\!\lambda\sum_{e\in\mathcal{E}}A_{e}\phi\big{(}\mathbf{u}^{(e)}\big{)}.}$

The objective function in (18) is the sum of two convex functions
$f(\mathbf{w})$ and $g(\mathbf{u})$ whose arguments are coupled as
$\mathbf{u}=\mathbf{D}\mathbf{w}$. Representing $f(\mathbf{w})$ and
$g(\mathbf{u})$ via their convex conjugates (see (17)) and interchanging the
minimization with the maximization (“taking the supremum”) in (17) results in
the dual problem

$\max_{\mathbf{u}\in\mathcal{U}}-g^{*}(\mathbf{u})-f^{*}(-\mathbf{D}^{T}\mathbf{u}).\vspace{-3mm}$
(19)

The objective function of the dual problem (19) is composed of the convex
conjugates

$\displaystyle g^{*}(\mathbf{u})$ $\displaystyle:=\sup_{{\bf
z}\in\mathcal{U}}\sum_{e\in\mathcal{E}}\big{(}\mathbf{u}^{(e)}\big{)}^{T}{\bf
z}^{(e)}-g({\bf z})$
$\displaystyle\stackrel{{\scriptstyle\eqref{equ_nLasso_compactly}}}{{=}}\sup_{{\bf
z}\in\mathcal{U}}\sum_{e\in\mathcal{E}}\big{(}\mathbf{u}^{(e)}\big{)}^{T}{\bf
z}^{(e)}\\!-\\!\lambda A_{e}\phi\big{(}{\bf z}^{(e)}\big{)}$
$\displaystyle=\sum_{e\in\mathcal{E}}\lambda
A_{e}\phi^{*}\big{(}\mathbf{u}^{(e)}/(\lambda A_{e})\big{)},\vspace{-3mm}$
(20)

and

$\displaystyle f^{*}(\mathbf{w})$ $\displaystyle\\!:=\\!\sup_{{\bf
z}\in\mathcal{W}}\sum_{i\in\mathcal{V}}\big{(}\mathbf{w}^{(i)}\big{)}^{T}{\bf
z}^{(i)}-f({\bf z})$
$\displaystyle\stackrel{{\scriptstyle\eqref{equ_nLasso_compactly}}}{{=}}\sup_{{\bf
z}\in\mathcal{W}}\sum_{i\in\mathcal{V}}\big{(}\mathbf{w}^{(i)}\big{)}^{T}{\bf
z}^{(i)}-\sum_{i\in\mathcal{V}}L_{i}\left(\mathbf{w}^{(i)}\right)$
$\displaystyle=\sum_{i\in\mathcal{V}}L^{*}_{i}\left(\mathbf{w}^{(i)}\right).\vspace{-3mm}$
(21)

The domain of the dual problem (19) is the space $\mathcal{U}$ of maps
$\mathbf{u}:\mathcal{E}\rightarrow\mathbb{R}^{d}$ that assign a flow vector
$\mathbf{u}^{(e)}$ to each edge $e\in\mathcal{E}$ of the empirical graph
$\mathcal{G}$. Note that the convex conjugates (III-B) and (III-B)

The duality between (16) and (19) is made precise in [44, Ch. 31] (see also
[29, Sec. 3.5]). First, the optimal values of both problems coincide [44, Cor.
31.2.1],

$\min_{\mathbf{w}\in\mathcal{W}}f(\mathbf{w})\\!+\\!g(\mathbf{D}\mathbf{w})\\!=\\!\max_{\mathbf{u}\in\mathcal{U}}-g^{*}(\mathbf{u})\\!-\\!f^{*}(-\mathbf{D}^{T}\mathbf{u}).\vspace{-2mm}$
(22)

Moreover, a necessary and sufficient condition for $\widehat{\mathbf{w}}$ to
solve (16) and $\widehat{\mathbf{u}}$ to solve (19) is [44, Thm. 31.3] (see
also [40, Ch. 7])

$-\mathbf{D}^{T}\widehat{\mathbf{u}}\in\partial f(\widehat{\mathbf{w}})\mbox{
, and }\mathbf{D}\widehat{\mathbf{w}}\in\partial
g^{*}(\widehat{\mathbf{u}}).\vspace{-2mm}$ (23)

The identity (22) (which is an instance of “strong duality”) allows to bound
the sub-optimality of some given local model parameters $\mathbf{w}^{(i)}$.
Indeed, for any given dual variable $\widehat{\mathbf{u}}\in\mathcal{U}$, the
objective function value
$-g^{*}(\widehat{\mathbf{u}})\\!-\\!f^{*}(-\mathbf{D}^{T}\widehat{\mathbf{u}})$
is a lower bound for the optimal value of (16). Section IV will discuss how
this fact can be used to constructing stopping criteria for iterative
optimization methods for GTV minimization.

It is instructive to rewrite the dual problem (19) and the optimality
condition (23) in terms of the local parameter vectors $\mathbf{w}^{(i)}$, for
each node $i\in\mathcal{V}$, and the local flow vectors $\mathbf{u}^{(e)}$,
for each $e\in\mathcal{E}$. Indeed, the final expressions in (III-B) and
(III-B) allow us to rewrite the dual problem (19) as

$\displaystyle\max_{\mathbf{u}\in\mathcal{U}}$
$\displaystyle{-\sum_{i\in\mathcal{V}}L^{*}_{i}\left(\mathbf{w}^{(i)}\right)-\lambda\sum_{e\in\mathcal{E}}A_{e}\phi^{*}\big{(}\mathbf{u}^{(e)}/(\lambda
A_{e})\big{)}}$ (24) $\displaystyle{\mbox{ subject to
}-\mathbf{w}^{(i)}=\sum_{e\in\mathcal{E}}\sum_{i=e_{+}}\mathbf{u}^{(e)}-\sum_{i=e_{-}}\mathbf{u}^{(e)}}$
$\displaystyle{\mbox{ for all nodes }i\in\mathcal{V}.}$

Using the definition of the block-incidence matrix (6) and its transpose (7),
we can also rewrite the optimality condition (23) as

$\displaystyle\sum_{e\in\mathcal{E}}\sum_{i=e_{+}}$
$\displaystyle{\widehat{\mathbf{u}}^{(e)}-\sum_{i=e_{-}}\widehat{\mathbf{u}}^{(e)}=-\nabla
L_{i}\left(\widehat{\mathbf{w}}^{(i)}\right)}$ (25) $\displaystyle{\mbox{ for
all nodes }i\in\mathcal{V}}$ $\displaystyle\widehat{\mathbf{w}}^{(e_{+})}-$
$\displaystyle{\widehat{\mathbf{w}}^{(e_{-})}\in(\lambda
A_{e})\partial\phi^{*}(\widehat{\mathbf{u}}^{(e)}/(\lambda A_{e}))}$
$\displaystyle{\mbox{ for every edge }e\in\mathcal{E}.}$

Let us now specialize the dual problem (24) for a GTV penalty function $\phi$
being a norm $\|\cdot\|$ on $\mathbb{R}^{d}$ [45]. For such a penalty function
$\phi\big{(}\mathbf{u}^{(e)}\big{)}=\left\|{\mathbf{u}^{(e)}}\right\|$, the
convex conjugate is the indicator of the dual-norm ball [35, Example 3.26],

$\phi^{*}\big{(}\mathbf{u}^{(e)}\big{)}=\begin{cases}0&\mbox{ for
}\big{\|}\mathbf{u}^{(e)}\big{\|}_{*}\leq 1\\\ \infty&\mbox{
else.}\end{cases}\vspace{-2mm}$ (26)

Inserting (26) into (24),

$\displaystyle\max_{\mathbf{u}\in\mathcal{U}}-\sum_{i\in\mathcal{V}}L^{*}_{i}\left(\mathbf{w}^{(i)}\right)$
$\displaystyle\mbox{ subject to }-\mathbf{w}^{(i)}$
$\displaystyle=\sum_{e\in\mathcal{E}}\sum_{i=e_{+}}\mathbf{u}^{(e)}-\sum_{e\in\mathcal{E}:i=e_{-}}\mathbf{u}^{(e)}$
$\displaystyle\quad\quad\mbox{ for all nodes }i\in\mathcal{V}$
$\displaystyle\|\mathbf{u}^{(e)}\|_{*}$ $\displaystyle\leq\lambda A_{e}\mbox{
for all edges }e\in\mathcal{E}.\vspace{-2mm}$ (27)

Thus, when the GTV penalty function is a norm $\phi(\cdot)=\|\cdot\|$, the
optimality condition (25) becomes (see [44, p. 215])

$\displaystyle\sum_{e\in\mathcal{E}:i\\!=\\!e_{+}}\widehat{\mathbf{u}}^{(e)}-\sum_{e\in\mathcal{E}:i\\!=\\!e_{-}}$
$\displaystyle\widehat{\mathbf{u}}^{(e)}=-\nabla
L_{i}\left(\widehat{\mathbf{w}}^{(i)}\right)\mbox{ for all nodes
}i\in\mathcal{V},$
$\displaystyle\quad\quad\|\widehat{\mathbf{u}}^{(e)}\|_{*}\leq\lambda A_{e}$
$\displaystyle\mbox{ for every edge }e\in\mathcal{E},$ (28)
$\displaystyle\widehat{\mathbf{w}}^{(e_{+})}=\widehat{\mathbf{w}}^{(e_{-})}\mbox{
for}$ $\displaystyle\mbox{ every edge }e\in\mathcal{E}\mbox{ with
}\|\widehat{\mathbf{u}}^{(e)}\|_{*}<\lambda A_{e}.\vspace{-2mm}$

###  III-C Interpretations

One important special case of GTV minimization (16) is convex clustering [34].
Indeed, convex clustering is obtained from (16) using the local loss functions

$L_{i}\left(\mathbf{w}^{(i)}\right)=\|{\mathbf{w}^{(i)}-{\bf
a}^{(i)}}\|^{2},\mbox{ for all nodes }i\in\mathcal{V}\vspace{-1mm}$ (29)

and GTV penalty $\phi({\bf u})=\left\|{{\bf u}}\right\|_{p}$ being a $p$-norm
$\left\|{{\bf u}}\right\|_{p}:=\big{(}\sum_{j=1}^{d}|u_{j}|^{p}\big{)}^{1/p}$
with some $p\geq 1$. The vectors ${\bf a}^{(i)}$ in (29) are the observations
that we wish to cluster.

The dual problem (19) of GTV minimization (16) is closely related to network
flow optimization. Indeed, the dual problem in the form (24) generalizes the
optimal flow problem [17, Sec. 1J] to vector-valued flows. The special case of
the dual problem (III-B), obtained when the GTV penalty function $\phi$ is a
norm, is equivalent to a generalized minimum-cost flow problem [16, Sec.
1.2.1]. Indeed, the maximization problem (III-B) is equivalent to the
minimization

$\displaystyle\min_{\mathbf{u}\in\mathcal{U}}\sum_{i\in\mathcal{V}}L^{*}_{i}\left(\mathbf{w}^{(i)}\right)$
subject to
$\displaystyle-\mathbf{w}^{(i)}=\sum_{e\in\mathcal{E}}\sum_{i=e_{+}}\mathbf{u}^{(e)}-\sum_{i=e_{-}}\mathbf{u}^{(e)}$
$\displaystyle\quad\quad\mbox{ for all nodes }i\in\mathcal{V}$
$\displaystyle\|\mathbf{u}^{(e)}\|_{*}\leq\lambda A_{e}\mbox{ for all edges
}e\in\mathcal{E}.\vspace{-2mm}$ (30)

The optimization problem (III-C) reduces to the minimum-cost flow problem [16,
Eq. (1.3) - (1.5)] for local model parameters of length $d=1$ (i.,e., local
models parametrized by a scalar).

Another interpretation of GTV minimization is as an instance of locally
weighted learning [46]. Indeed, our analysis in Section V reveals that, under
certain conditions on the connectivity of the empirical graph and local loss
functions, GTV minimization effectively implements cluster-wise optimization
(9). In other words, if the node $i$ belongs to the cluster $\mathcal{C}$, the
solution $\widehat{\mathbf{w}}^{(i)}$ of GTV approximates the solution
$\overline{{\bf w}}^{(\mathcal{C})}$ of (9),
$\widehat{\mathbf{w}}^{(i)}\approx\overline{{\bf w}}^{(\mathcal{C})}$. The
deviation between $\widehat{\mathbf{w}}^{(i)}$ and $\overline{{\bf
w}}^{(\mathcal{C})}$ will be bounded by Theorem 2. The cluster-wise
optimization (9) can be re-written as a locally weighted learning problem [46,
Sec. 3.1.2]

$\overline{{\bf w}}^{(\mathcal{C})}=\operatorname*{argmin}_{{\bf
w}\in\mathbb{R}^{d}}\sum_{i\in\mathcal{V}}\rho_{i}L_{i}\left({\bf
w}\right).\vspace{-2mm}$ (31)

with the weights $\rho_{i}=1$ for $i\in\mathcal{C}$ and $\rho_{i}=0$
otherwise. Note that (31) is an adaptive pooling of local datasets into a
cluster $\mathcal{C}$. The pooling of local datasets is driven jointly by the
geometry (connectivity) of the empirical graph $\mathcal{G}$ and the geometry
(shape) of local loss functions (see Theorem 2).

##  IV A Primal Dual Method for GTV Minimization

The GTV minimization problem (16) is a non-smooth convex optimization problem
with a particular structure. In particular, the objective function in (16)
consists of two components that could be optimized easily when considered
separately. Indeed, the first component
$\sum_{i\in\mathcal{V}}L_{i}\left(\mathbf{w}^{(i)}\right)$ could be minimized
trivially by separately minimizing $L_{i}\left(\mathbf{w}^{(i)}\right)$, for
each node $i\in\mathcal{V}$. The second component $\lambda\|\mathbf{w}\|_{\rm
GTV}$ is minimized trivially by using constant networked model parameters.

Primal-dual methods use tools from convex duality to solve problems with a
composite objective function such as (16) [44, 29]. We obtain Algorithm 1 by
applying the primal-dual method [29, Alg. 6] to jointly solve (16) and (19).
The details of this application are provided in Supplementary material IX-D.
Note that Algorithm 1 represents an entire family of NFL algorithms. This
family is parametrized by the local loss function $L_{i}\left(\cdot\right)$
and the GTV penalty function $\phi(\cdot)$ (see (14)).

At its core, Algorithm 1 computes and distributes (via the edges of the
empirical graph $\mathcal{G}$) the node-wise primal and the edge-wise dual
updates in steps (4) and (8), respectively. The primal update (operator) in
step (4) is operators

$\displaystyle\mathcal{PU}^{(i)}\left\\{{\bf v}\right\\}:=$
$\displaystyle\operatorname*{argmin}_{{\bf z}\in\mathbb{R}^{d}}L_{i}\left({\bf
z}\right)\\!+\\!\left(\big{|}\mathcal{N}^{(i)}\big{|}/2\right)\|{\bf
v}\\!-\\!{\bf z}\|^{2}$ (32) $\displaystyle\mbox{ for each node
}i\in\mathcal{V}.$

Here, we used the neighbourhood
$\mathcal{N}^{(i)}:=\\{i^{\prime}\in\mathcal{V}:\\{i,i^{\prime}\\}\in\mathcal{E}\\}$
of a node $i\in\mathcal{V}$. Comparing (32) with (3) reveals that the primal
update $\mathcal{PU}^{(i)}\left\\{\cdot\right\\}$ is exactly the proximity
operator of the local loss function $L_{i}\left(\cdot\right)$.

$\mathcal{PU}^{(i)}\left\\{{\bf v}\right\\}={\rm\bf
prox}_{L_{i}\left(\cdot\right),\rho}({\bf v})\mbox{ with
}\rho=\big{|}\mathcal{N}^{(i)}\big{|}.\vspace{-2mm}$

We can interpret the primal update (32) as regularized minimization of the
local loss $L_{i}\left(\cdot\right)$ [38]. The regularization term in (32) is
the squared Euclidean distance between local model parameters and the argument
${\bf v}$ of the primal update. It enforces similar local model parameters at
well-connected nodes in the same cluster. The amount of regularization is
controlled by the node degree $\big{|}\mathcal{N}^{(i)}\big{|}$. In
particular, the larger the node degree, the smaller the influence of the local
loss function on the resulting of the primal update (32). The primal update
(32) is parametrized by local loss function $L_{i}\left(\cdot\right)$.
Supplementary material IX-B works out the primal update for the local loss
functions used by some popular ML methods [38].

Algorithm 1 alternates between the primal updates (32), for each node
$i\in\mathcal{V}$, and dual updates

$\displaystyle\mathcal{DU}^{(e)}\left\\{{\bf v}\right\\}:=$
$\displaystyle\operatorname*{argmin}_{{\bf z}\in\mathbb{R}^{d}}\lambda
A_{e}\phi^{*}\big{(}{\bf z}/(\lambda A_{e})\big{)}\\!+\\!(1/2\sigma_{e})\|{\bf
v}\\!-\\!{\bf z}\|^{2}$ $\displaystyle\mbox{ for each
}e\in\mathcal{E}.\vspace{-2mm}$ (33)

Here, we used the convex conjugate $\phi^{*}({\bf v}):=\sup_{{\bf
z}\in\mathbb{R}^{d}}{\bf v}^{T}{\bf z}-\phi({\bf z})$ of the GTV penalty
function $\phi({\bf v})$ (see (14)). A comparison of (IV) with (3) reveals
that the dual update $\mathcal{DU}^{(e)}$ is essentially the (scaled)
proximity operator of the convex conjugate $\phi^{*}({\bf v})$,

$\mathcal{DU}^{(e)}\left\\{{\bf v}\right\\}=\lambda A_{e}{\rm\bf
prox}_{\phi^{*},\rho}({\bf v}/(\lambda A_{e}))\mbox{ with }\rho=\lambda
A_{e}/\sigma_{e}.\vspace{-2mm}$

We can interpret the dual update operator (IV) as a regularized minimization
of the convex conjugate $\phi^{*}({\bf v})$ of the GTV penalty function. The
regularization term in (IV) forces the update to not deviate too much from its
argument ${\bf v}$. Note that the dual update (IV) is parametrized by the GTV
penalty function $\phi({\bf v})$ (14) (via its convex conjugate).
Supplementary material IX-C will derive the following dual update operators
for some popular GTV penalty functions that have been studied recently in the
FL literature [9, 19, 20, 41].

The “MOCHA penalty” $\phi({\bf v}):=(1/2)\|{\bf v}\|^{2}_{2}$ [9], lends to
the dual update operator $\mathcal{DU}^{(e)}\big{\\{}{\bf v}\big{\\}}:={\bf
v}/\big{(}1+(\sigma_{e}/(\lambda A_{e}))\big{)}$. An obvious generalization of
the MOCHA penalty is $\phi({\bf v})\\!:=(1/2){\bf v}^{T}\mathbf{Q}{\bf v}$
with a fixed positive semidefinite matrix $\mathbf{Q}$. The corresponding dual
operator is $\mathcal{DU}^{(e)}\big{\\{}{\bf
v}\big{\\}}:=\left((\sigma_{e}/(\lambda
A_{e}))\mathbf{Q}^{-1}+\mathbf{I}\right)^{-1}{\bf v}$, which is a linear
operator.

For the “nLasso penalty” $\phi({\bf v})\\!:=\\!\|{\bf v}\|_{2}$ [19, 20], the
dual update operator becomes the (vector) clipping operator
$\mathcal{DU}^{(e)}\big{\\{}{\bf v}\big{\\}}:=\mathcal{T}^{(\lambda
A_{e})}\\{{\bf v}\\}$ (see (1)). The nLasso penalty is the Euclidean norm
$\left\|{{\bf v}}\right\|_{2}$ of the flow vector ${\bf v}$ across an edge.
Using instead the $\ell_{1}$ norm yields the GTV penalty function $\phi({\bf
v}):=\|{\bf v}\|_{1}$ [41]whose associated dual update operator
$\mathcal{DU}^{(e)}\big{\\{}{\bf v}\big{\\}}:=\big{(}\mathcal{T}^{(\lambda
A_{e})}\big{(}v_{1}),\ldots,\mathcal{T}^{(\lambda
A_{e})}\big{(}v_{d})\big{)}^{T}$ is an element-wise application of the scalar
clipping operator.

Algorithm 1 Primal-Dual Method for Networked FL

Input: empirical graph $\mathcal{G}$; local loss
$\big{\\{}L_{i}\left(\cdot\right)\big{\\}}_{i\in\mathcal{V}}$, for
$i\in\mathcal{V}$, GTV parameter $\lambda$ and penalty $\phi(\cdot)$
Initialize: $k:=0$; $\widehat{\mathbf{w}}^{(i)}_{0}:={\bf
0},\tau_{i}=1/|\mathcal{N}^{(i)}|$ for all nodes $i\in\mathcal{V}$;
$\widehat{\mathbf{u}}^{(e)}_{0}:={\bf 0},\sigma_{e}=1/2$ for all edges
$e\in\mathcal{E}$;

1:while stopping criterion is not satisfied do

2: for all nodes $i\in\mathcal{V}$ do

3: $\widehat{{\bf w}}_{k+1}^{(i)}:=\widehat{{\bf
w}}^{(i)}_{k}-\tau_{i}\sum_{e\in\mathcal{E}}D_{e,i}\widehat{\mathbf{u}}^{(e)}_{k}$

4: $\widehat{{\bf w}}_{k+1}^{(i)}:=\mathcal{PU}^{(i)}\left\\{\widehat{{\bf
w}}_{k+1}^{(i)}\right\\}$

5: end for

6: for all edges $e\in\mathcal{E}$ do

7:
$\widehat{\mathbf{u}}^{(e)}_{k+1}:=\widehat{\mathbf{u}}^{(e)}_{k}\\!+\\!\sigma_{e}\big{(}2\big{(}\widehat{{\bf
w}}^{(e_{+})}_{k\\!+\\!1}\\!-\\!\widehat{{\bf
w}}^{(e_{-})}_{k+1}\big{)}\\!-\\!\big{(}\widehat{{\bf
w}}^{(e_{+})}_{k}\\!-\\!\widehat{{\bf w}}^{(e_{-})}_{k}\big{)}\big{)}$

8:
$\widehat{\mathbf{u}}^{(e)}_{k\\!+\\!1}\\!:=\\!\mathcal{DU}^{(e)}\big{\\{}\widehat{\mathbf{u}}^{(e)}_{k\\!+\\!1}\big{\\}}$

9: end for

10: $k\\!:=\\!k\\!+\\!1$

11:end while

Algorithm 1 can be implemented as message passing over the edges of the
empirical graph $\mathcal{G}$. During each iteration, local computations are
carried out at each node and each edge. The primal update step (4) can be
executed in parallel for each node $i\in\mathcal{V}$. Note that (4) involves
the local loss function $L_{i}\left(\mathbf{w}^{(i)}\right)$ which
encapsulates the local dataset $\mathcal{X}^{(i)}$. The dual update step (8)
can be carried out in parallel for each edge $e\in\mathcal{E}$. Note that (8)
is parametrized by the GTV penalty function $\phi$ (see (14)). The results of
the node-wise primal and edge-wise dual updates in step (4) and (8) are spread
to neighbouring (incident) edges and nodes during steps (3) and (7) of
Algorithm 1. Trivially, the number of basic computational steps required by
Algorithm 1 is directly proportional to the number of edges in the empirical
graph $\mathcal{G}$.

It can be shown that Algorithm 1 is robust against errors occurring during the
updates in steps (4) and (8) [32, 23]. This is important for applications
where the update operators (32) and (IV) can be evaluated only approximately
or these updates have to be transmitted over imperfect wireless links. Let us
denote the perturbed update (e.g., obtained by a numerical optimization method
for solving (32)) by $\widetilde{{{\bf w}}}_{k\\!+\\!1}$ and the exact update
by $\widehat{\mathbf{w}}^{(i)}_{k+1}$, respectively. Then, Algorithm 1 is
still guaranteed to converge to a solution of (16) as long as
$\sum_{k=1}^{\infty}\left\|{\widetilde{{{\bf
w}}}_{k\\!+\\!1}-\widehat{\mathbf{w}}^{(i)}_{k+1}}\right\|<\infty$ (see [23,
Sec. 3]).

Possible stopping criteria in step (1) of Algorithm 1 include the use of a
fixed number of iterations. It might also be useful to stop iterating when the
decrease in the objective function (16) falls below a threshold. The
construction of a stopping criterion can also be based on the primal-dual gap

$\displaystyle{\rm gap}_{k}:=$
$\displaystyle\sum_{i\in\mathcal{V}}L_{i}\left(\widehat{{\bf
w}}^{(i)}_{k+1}\right)+\lambda\|\widehat{{\bf w}}_{k+1}\|_{\rm GTV}-$ (34)
$\displaystyle\left(-\sum_{i\in\mathcal{V}}L^{*}_{i}\left(-\mathbf{s}^{(i)}\right)-\lambda\sum_{e\in\mathcal{E}}A_{e}\phi^{*}\big{(}\widehat{\mathbf{u}}^{(e)}_{k\\!+\\!1}/(\lambda
A_{e})\big{)}\right)$ $\displaystyle\mbox{ with
}\mathbf{s}^{(i)}:=\sum_{e\in\mathcal{E}}\sum_{i=e_{+}}\widehat{\mathbf{u}}^{(e)}_{k\\!+\\!1}-\sum_{i=e_{-}}\widehat{\mathbf{u}}^{(e)}_{k\\!+\\!1}.$

By comparing (34) with (22), we can bound the sub-optimality of the iterate
$\widehat{{\bf w}}^{(i)}_{k+1}$ as

$\displaystyle\sum_{i\in\mathcal{V}}$ $\displaystyle L_{i}\left(\widehat{{\bf
w}}^{(i)}_{k\\!+\\!1}\right)+\lambda\|\widehat{{\bf w}}_{k+1}\|_{\rm GTV}$
$\displaystyle-\min_{\mathbf{w}\in\mathcal{W}}\left[\sum_{i\in\mathcal{V}}L_{i}\left({\bf
w}^{(i)}\right)+\lambda\|\mathbf{w}\|_{\rm GTV}\right]\leq{\rm gap}_{k}.$

Thus, to ensure that Algorithm 1 delivers local model parameters with sub-
optimality no larger than $\eta$, we only stop when ${\rm gap}_{k}\leq\eta$.
If the local loss functions $L_{i}\left(\cdot\right)$ and the GTV penalty
$\phi$ are convex, this condition is always fulfilled after a finite number of
iterations [29, Thm. 5.1].

##  V When Does It Work?

Section III developed Algorithm 1 as a novel FL method for distributed
training of tailored models from collections of local datasets. We obtained
Algorithm 1 by using the primal-dual method [29, Alg. 6] to solve the GTV
minimization problem (16) jointly with its dual (24). The rationale behind
Algorithm 1, and GTV minimization (16) in the first place, was that local loss
functions are clustered according to Assumption 1 and that this cluster
structure is reflected by the edges in the empirical graph $\mathcal{G}$.

To analyze the statistical properties of Algorithm 1, we must make precise the
relation between the network structure of $\mathcal{G}$ and the cluster
structure (see (10)) of local loss functions. To this end, we introduce the
notion of a well-connected cluster $\mathcal{C}\subseteq\mathcal{V}$ whose
nodes share (approximately) a common optimal choice for local model parameters
(see Assumption 1).

######  Definition 1 (Well-Connected Cluster).

Consider an empirical graph $\mathcal{G}=\left(\mathcal{V},\mathcal{E}\right)$
with edge weights $A_{i,i^{\prime}}$ for $i,i^{\prime}\in\mathcal{V}$. The
nodes $i\in\mathcal{V}$ carry local loss functions
$L_{i}\left(\mathbf{w}^{(i)}\right)$. that are partitioned into clusters
$\mathcal{P}=\\{\mathcal{C}^{(1)},\ldots,\mathcal{C}^{(k)}\\}$ according to
Assumption 1 with clustering error $\delta^{(i)}$ (see (11)). By Assumption 2,
we also require local loss functions to be smooth with Lipschitz constant
$\beta^{(i)}$ (see (12)) and their cluster-wise sum to be strongly convex with
parameter $\alpha^{(c)}$ (see (13)). A cluster
$\mathcal{C}^{(c)}\in\mathcal{P}$ (see (10)), with weighted boundary
$\big{|}\partial\mathcal{C}^{(c)}\big{|}:=\sum_{i\in\mathcal{C}^{(c)},i^{\prime}\notin\mathcal{C}^{(c)}}A_{i,i^{\prime}}$,
is well-connected if it contains a node $i_{0}\in\mathcal{C}^{(c)}$ such that

$\displaystyle\sum_{i\in\mathcal{A}}\big{[}\sum_{i^{\prime}\notin\mathcal{C}^{(c)}}A_{i,i^{\prime}}\big{]}$
$\displaystyle+\beta^{(i)}\big{|}\partial\mathcal{C}^{(c)}\big{|}/\alpha^{(c)}$
(35)
$\displaystyle+\delta^{(i)}/\lambda<\sum_{i\in\mathcal{A}}\sum_{i^{\prime}\in\mathcal{C}^{(c)}\setminus\mathcal{A}}A_{i,i^{\prime}}$
$\displaystyle\mbox{ for all subsets
}\mathcal{A}\subseteq\mathcal{C}^{(c)}\setminus\big{\\{}i_{0}\big{\\}}.$

Our main result is that GTV minimization (16) captures an underlying partition
of local datasets into nearly homogeneous data (see Assumption 1) if every
cluster is well-connected according to Definition 1. More precisely, we have
the following upper bound on the deviation between the solutions of GTV
minimization (16) and the cluster-wise optimizers (9).

######  Theorem 2.

Consider networked data with empirical graph $\mathcal{G}$ whose nodes carry
local loss functions that are clustered into
$\mathcal{P}=\\{\mathcal{C}^{(1)},\ldots,\mathcal{C}^{(k)}\\}$ according to
Assumption 1 and satisfy Assumption 2. We learn local model parameters
$\widehat{\mathbf{w}}^{(i)}$ by solving GTV minimization (16) with penalty
function $\phi$ being a norm.

If every cluster $\mathcal{C}^{(c)}\in\mathcal{P}$ is well-connected, then

  * •

the learnt parameter vectors are constant over nodes in the same cluster
$\mathcal{C}^{(c)}\in\mathcal{P}$,
$\widehat{\mathbf{w}}^{(i)}=\widehat{\mathbf{w}}^{(i^{\prime})}\mbox{ for
}i,i^{\prime}\in\mathcal{C}^{(c)}.$

  * •

the deviation between the GTV solution $\widehat{\mathbf{w}}^{(i)}$ at some
node $i\in\mathcal{C}^{(c)}$ and the cluster-wise optimizer $\overline{{\bf
w}}^{(c)}$ in (9) (for $\mathcal{C}=\mathcal{C}^{(c)}$) is bounded as

$\left\|{\widehat{\mathbf{w}}^{(i)}-\overline{{\bf w}}^{(c)}}\right\|\leq
2\big{|}\partial\mathcal{C}^{(c)}\big{|}\lambda/\alpha^{(c)}.\vspace{-2mm}$
(36)

###### Proof.

see Supplementary material IX-A. ∎

The error bound (36) involves the shape of the local loss functions via the
strong convexity parameter $\alpha^{(c)}$ for cluster $\mathcal{C}^{(c)}$.
Note that Assumption 2 requires, for each cluster
$\mathcal{C}^{(c)}\in\mathcal{P}$ of the partition (10), the cluster-wise
aggregate loss function
$\sum_{i\in\mathcal{C}^{(c)}}L_{i}\left(\mathbf{w}^{(i)}\right)$ to be
strongly convex.

##  VI Numerical Experiments

This section reports the results of some illustrative numerical experiments to
verify the performance of Algorithm 1. We provide the code to reproduce these
experiments in the supplementary material. The experiments discussed in
Section VI-A \- Section IX-E2 revolve around synthetic datasets whose
empirical graph is either a chain graph, a star graph or a realization of a
Stochastic Block Model (SBM) [47]. A numerical experiment with a handwritten
digit image dataset (“MNIST”) is discussed in Section VI-C. Section IX-E3
discusses a numerical experiment based on open weather data provided by the
Finnish Meteorological Institute (FMI).

###  VI-A Stochastic Block Model

This experiment revolves around synthetic data whose empirical graph
$\mathcal{G}^{(\rm SBM)}$ is partitioned into two equal-sized clusters
$\mathcal{P}=\\{\mathcal{C}^{(1)},\mathcal{C}^{(2)}\\}$, with
$|\mathcal{C}^{(1)}|=|\mathcal{C}^{(2)}|$. We denote the cluster assignment of
node $i\in\mathcal{V}$ by $c^{(i)}\in\\{1,2\\}$. The edges in
$\mathcal{G}^{(\rm SBM)}$ are generated via realizations of independent binary
random variables $b_{i,i^{\prime}}\in\\{0,1\\}$. These random variables are
indexed by pairs $i,i^{\prime}$ of nodes that are connected by an edge
$\\{i,i^{\prime}\\}\in\mathcal{E}$ if and only if $b_{i,i^{\prime}}=1$. Two
nodes in the same cluster are connected with probability ${\rm
Prob}\\{b_{i,i^{\prime}}=1\\}:=p_{\rm in}$ if $i,i^{\prime}$. In contrast,
${\rm Prob}\\{b_{i,i^{\prime}}=1\\}:=p_{\rm out}$ if nodes $i,i^{\prime}$
belong to different clusters. Every edge in $\mathcal{G}^{(\rm SBM)}$ has the
same weight, $A_{e}=1$ for all $e\in\mathcal{E}$.

Each node $i\in\mathcal{V}$ of the empirical graph $\mathcal{G}^{(\rm SBM)}$
holds a local dataset $\mathcal{X}^{(i)}$ of the form (4). Thus, the dataset
$\mathcal{X}^{(i)}$ consists of $m_{i}$ data points, each characterized by a
feature vector ${\bf x}^{(i,r)}\in\mathbb{R}^{d}$ and a scalar label
$y^{(i,r)}$, for $r=1,\ldots,m_{i}$. The feature vectors ${\bf
x}^{(i,r)}\sim\mathcal{N}(\mathbf{0},\mathbf{I}_{d\times d})$, are drawn
i.i.d. from a standard multivariate normal distribution. The labels of the
data points are generated by a noisy linear model

$y^{(i,r)}=\big{(}\overline{{\bf w}}^{(i)}\big{)}^{T}{\bf
x}^{(i,r)}+\sigma\varepsilon^{(i,r)}.\vspace{-2mm}$ (37)

The noise $\varepsilon^{(i,r)}\sim\mathcal{N}(0,1)$, for $i\in\mathcal{V}$ and
$r=1,\ldots,m_{i}$, are i.i.d. realizations of a standard normal distribution.
The true underlying weight vector $\overline{{\bf w}}^{(i)}$ is piece-wise
constant over the clusters in the partition
$\mathcal{P}=\\{\mathcal{C}^{(1)},\mathcal{C}^{(2)}\\}$, i.e., $\overline{{\bf
w}}^{(i)}=\overline{{\bf w}}^{(i^{\prime})}$ if
$i,i^{\prime}\in\mathcal{C}^{(c)}$.

To study the robustness of Algorithm 1 against node failures, we assume that
local datasets are only accessible in a subset
$\mathcal{M}\subseteq\mathcal{V}$. The set $\mathcal{M}$ is selected uniformly
at random among all nodes $\mathcal{V}$. We can access local datasets
$\mathcal{X}^{(i)}$ only in the subset $\mathcal{M}$ of relative size
$\rho:=|\mathcal{M}|/|\mathcal{V}|$. The inaccessibility of a local dataset
$\mathcal{X}^{(i)}$, for $i\notin\mathcal{M}$, can be modelled by using a
trivial loss function. Thus, we learn local model parameters
$\widehat{\mathbf{w}}^{(i)}$ using Algorithm 1 with local loss functions

$L_{i}\left(\mathbf{w}^{(i)}\right):=\begin{cases}\frac{1}{m_{i}}\sum_{r=1}^{m_{i}}\big{(}\big{(}{\bf
x}^{(i,r)}\big{)}^{T}\mathbf{w}^{(i)}-y^{(i,r)}\big{)}^{2}\mbox{ for
}i\in\mathcal{M}\\\ 0\mbox{ otherwise.}\end{cases}.$ (38)

The special case $\rho=1$ is obtained when all local datasets are accessible,
i.e., when $\mathcal{M}=\mathcal{V}$. For the stopping criterion in Algorithm
1 we use a fixed number $R$ of iterations, which is $R=3000$ for the results
depicted in Figure 7 and $R\\!=\\!2000$ for the results depicted in Figure 8.
The GTV regularization parameter has been set to $\lambda\\!=\\!10^{-2}$ (see
(8)). We measure the estimation error incurred by the learnt parameter vectors
$\widehat{\mathbf{w}}^{(i)}$ using the average squared estimation error (MSE),

${\rm
MSE}:=(1/|\mathcal{V}|)\sum_{i\in\mathcal{V}}\|\widehat{\mathbf{w}}^{(i)}-\overline{{\bf
w}}^{(i)}\|_{2}^{2}.\vspace{-2mm}$ (39)

####  VI-A1 High-Dimensional Linear Regression

Table I reports the results obtained by applying Algorithm 1 to
$\mathcal{G}^{(\rm SBM)}$ with $|\mathcal{C}^{(1)}|=|\mathcal{C}^{(2)}|=100$,
$p_{\rm in}=0.5$, and $p_{\rm out}=10^{-2}$. Each node carries a local dataset
of the form (4) with $m_{i}\\!=\\!10$ data points, having a feature vector
${\bf x}^{(i,r)}\in\mathbb{R}^{100}$ and labels generated according to (37)
with noise strength $\sigma=10^{-3}$. The true underlying parameter vectors in
(37) are cluster-wise constant, $\overline{{\bf w}}^{(i)}=\overline{{\bf
w}}^{(c)}$ for all $i\in\mathcal{C}^{(c)}$. The cluster-wise parameter vector
$\overline{{\bf w}}^{(c)}\in\\{0,0.5\\}^{100}$ is constructed entry-wise using
i.i.d. realizations of standard Bernoulli variables $B\in\\{0,0.5\\}$ with
${\rm Prob}(B=0)=1/2$.

We learn local model parameters using Algorithm 1 using local loss functions
(4) and assuming that each local dataset is available,
$\mathcal{M}=\mathcal{V}$. Table I reports the MSE (39) incurred by Algorithm
1 using a fixed number of $R=1000$ iterations. Table I also reports the MSE
(39) incurred by the model parameters learnt by IFCA [12] and FedAvg [48].

Method | MSE
---|---
Algorithm 1 | 1.42e-05
IFCA | 2.82
FedAvg | 2.86
TABLE I: MSE (39) of the estimation error incurred by Algorithm 1, IFCA [12]
and FedAvg [48]

###  VI-B Chain Graph

This experiment uses a dataset with $2n$ local datasets whose empirical graph
$\mathcal{G}^{(\rm chain)}$ is a chain graph. The edge set of this graph is
given by $\mathcal{E}=\\{\\{i,i+1\\}:i\in\\{1,\ldots,2n-1\\}\\}$. The nodes
$\mathcal{V}=\\{1,\ldots,2n\\}$ are partitioned into two clusters
$\mathcal{C}^{(1)}=\\{1,\ldots,n\\}$ and
$\mathcal{C}^{(2)}=\\{n+1,\ldots,2n\\}$. Every intra-cluster edge
$\\{i,i^{\prime}\\}\in\mathcal{E}$, with $i,i^{\prime}$ in the same cluster,
has weight $A_{i,i^{\prime}}=1$. The single inter-cluster edge
$e^{\prime}=\\{n,n+1\\}$ has weight $A_{n,n+1}=\varepsilon$ with some
$\varepsilon\geq 0$.

The nodes of $\mathcal{G}^{(\rm chain)}$ carry local datasets with the same
structure as in Section VI-A. In particular, each local dataset consists of
data points that are characterized by feature vectors drawn from a
multivariate normal distribution and a label that is generated from (37). The
true local parameter vectors (“ground truth”) in (37) are piece-wise constant
over clusters $\mathcal{C}^{(1)}$ and $\mathcal{C}^{(2)}$.

We use Algorithm 1 to learn the local parameter vectors in (37) using local
loss function obtained from the average squared error (38). The local datasets
are only accessible for a subset $\mathcal{M}\subseteq\mathcal{V}$ of
$\lceil\rho|\mathcal{V}|\rceil$ nodes selected uniformly at random among all
nodes $\mathcal{V}$. We apply Algorithm 1 to the local loss functions (38) and
using a fixed number $R\\!=\\!2000$ of iterations and different choices for
the GTV parameter $\lambda$ as indicated in Figure 2 and Figure 3.

Figure 2 and Figure 3 depict the MSE (39) incurred by Algorithm 1 obtained for
varying training size $\rho$, noise strength $\sigma$ (see (37)) and inter-
cluster edge weight $\varepsilon$. The curves and bars in Figures 2 and 3
represent the average and standard deviation of MSE values of $5$ simulation
runs for different choices for $\rho,\varepsilon$ and $\sigma$. Figure 2 shows
the results in the noise-less case where $\sigma=0$ in (37), and Figure 3
shows results for varying noise level $\sigma$ in (37) and fixed relative
training size $\rho=6/10$.

$0$$0.2$$0.4$$0.6$$0.8$$1$$1.2$$10^{-10}$$10^{-8}$$10^{-6}$$10^{-4}$$10^{-2}$$10^{0}$$\varepsilon$MSE$\phi({\bf
u})\\!=\\!\|{\bf
u}\|_{1},\lambda\\!=\\!0.1$$\rho\\!=\\!2/10$$\rho\\!=\\!4/10$$\rho\\!=\\!6/10$

$0$$0.2$$0.4$$0.6$$0.8$$1$$1.2$$10^{-10}$$10^{-8}$$10^{-6}$$10^{-4}$$10^{-2}$$10^{0}$$\varepsilon$MSE$\phi({\bf
u})\\!=\\!\|{\bf u}\|_{2},\lambda\\!=\\!0.1$
(nLasso)$\rho\\!=\\!2/10$$\rho\\!=\\!4/10$$\rho\\!=\\!6/10$

$0$$0.2$$0.4$$0.6$$0.8$$1$$1.2$$10^{-10}$$10^{-8}$$10^{-6}$$10^{-4}$$10^{-2}$$10^{0}$$\varepsilon$MSE$\phi({\bf
u})\\!=\\!\|{\bf u}\|_{2}^{2},\lambda\\!=\\!0.05$
(“MOCHA”)$\rho\\!=\\!0.2$$\rho\\!=\\!0.4$$\rho\\!=\\!0.6$

Figure 2:  MSE (39) of the estimation error incurred by Algorithm 1 when
learning local model parameters for a chain graph $\mathcal{G}^{(\rm chain)}$
in the noiseless case $\sigma=0$.

$0$$0.2$$0.4$$0.6$$0.8$$1$$1.2$$10^{-7}$$10^{-6}$$10^{-5}$$10^{-4}$$10^{-3}$$10^{-2}$$10^{-1}$$10^{0}$$\epsilon$MSE$\phi({\bf
u})\\!=\\!\|{\bf
u}\|_{1},\lambda\\!=\\!0.1$$\sigma\\!=\\!10^{-2}$$\sigma\\!=\\!10^{-1}$$\sigma\\!=\\!1$

$0$$0.2$$0.4$$0.6$$0.8$$1$$1.2$$10^{-7}$$10^{-6}$$10^{-5}$$10^{-4}$$10^{-3}$$10^{-2}$$10^{-1}$$10^{0}$$\varepsilon$MSE$\phi({\bf
u})\\!=\\!\|{\bf u}\|_{2},\lambda\\!=\\!0.1$
(nLasso)$\sigma\\!=\\!10^{-2}$$\sigma\\!=\\!10^{-1}$$\sigma\\!=\\!1$

$0$$0.2$$0.4$$0.6$$0.8$$1$$1.2$$10^{-7}$$10^{-6}$$10^{-5}$$10^{-4}$$10^{-3}$$10^{-2}$$10^{-1}$$10^{0}$$\varepsilon$MSE$\phi({\bf
u})\\!=\\!\|{\bf u}\|^{2}_{2},\lambda\\!=\\!0.05$
(“MOCHA”)$\sigma\\!=\\!10^{-2}$$\sigma\\!=\\!10^{-1}$$\sigma\\!=\\!1$

Figure 3:  ${\rm MSE}$ (39) incurred by Algorithm 1 when learning the local
model parameters for the chain graph $\mathcal{G}^{(\rm chain)}$ using only
local datasets at a fraction of $|\mathcal{M}|/|\mathcal{V}|=0.6$ nodes (see
(38)).

###  VI-C Handwritten Digits

This experiment revolves around a collection of $n=40$ local datasets
generated from the handwritten image dataset (“MNIST dataset”) [49]. and
represented by the empirical graph $\mathcal{G}^{(\rm
MNIST)}=\left(\mathcal{V}^{(\rm MNIST)},\mathcal{E}^{(\rm MNIST)}\right)$.
Each node $i\in\mathcal{V}^{(\rm MNIST)}=\\{1,\ldots,40\\}$ carries a local
dataset $\mathcal{X}^{(i)}$ which consists of $m_{i}=500$ data points being
images of handwritten digits. The $r$-th data point in $\mathcal{X}^{(i)}$ is
characterized by the feature vector ${\bf x}^{(i,r)}\in\mathbb{R}^{d}$ and a
label that specifies the digit that the datapoint belongs to. The entries of
the feature vector ${\bf x}^{(i,r)}\in\mathbb{R}^{d}$ are greyscale levels of
$d=28\times 28$ pixels. The nodes $\mathcal{V}^{(\rm MNIST)}$ are partitioned
into two clusters $\mathcal{C}^{(1)},\mathcal{C}^{(2)}$ (see Assumption 1),
each with $20$ nodes. The nodes in $\mathcal{C}^{(1)}$ carry local datasets
that consist of images depicting handwritten digits $0$ and $1$. The nodes
$i\in\mathcal{C}^{(2)}$ in the second cluster carry local datasets
$\mathcal{X}^{(i)}$ consisting of images that depict handwritten digits $2$
and $3$.

Besides the local dataset $\mathcal{X}^{(i)}$, the node $i\in\mathcal{V}$ is
also assigned a local model in the form of an artificial neural network (ANN),

$\displaystyle h^{(i)}\big{(}{\bf x};\mathbf{w}^{(i)}\big{)}:={\rm
SoftMax}\big{(}{\bf W}^{(i,2)}{\rm ReLU}\big{(}{\bf W}^{(i,1)}{\bf
x}\big{)}\big{)}$ (40) $\displaystyle\mbox{ with }\mathbf{w}^{(i)}={\rm
stack}\\{{\bf W}^{(i,1)},{\bf W}^{(i,2)}\\}.$

The ANN (40) consists of two densely connected layers, with the first layer
using rectified linear units ReLU [25]. The second (output) layer uses a
“soft-max” activation function [25]. The weights of connections between layers
are stored in the matrices ${\bf W}^{(i,1)}$ and ${\bf W}^{(i,2)}$,
respectively. The entries of the weight matrices are collected in the local
model parameter vector $\mathbf{w}^{(i)}$.

For learning the local model parameters $\mathbf{w}^{(i)}$ of the ANNs (40),
we split each local dataset a training and validation set,
$\mathcal{X}^{(i)}=\mathcal{X}_{\rm train}^{(i)}\cup\mathcal{X}_{\rm
val}^{(i)}$ of size $m^{(\rm train)}_{i}\\!=\\!400$ and $m^{(\rm
val)}_{i}\\!=\\!100$, respectively. The local training sets $\mathcal{X}_{\rm
train}^{(i)}$ are also used for the construction of the edges in
$\mathcal{G}^{(\rm MNIST)}$ as described next.

For each $i\in\mathcal{V}^{(\rm MNIST)}$, we apply t-SNE [50] to map the raw
feature vector ${\bf x}^{(i,r)}$ to the embedding vector ${\bf
z}^{(i,r)}\in\mathbb{R}^{2}$, for $r=1,\ldots,m^{(\rm train)}_{i}$. We then
compute a distance ${\rm dist}(i,i^{\prime})$ between any two different nodes
$i,i^{\prime}\in\mathcal{V}^{(\rm MNIST)}$ using a Kullback-Leibler divergence
estimator [51].

Given the pairwise distances ${\rm dist}(i,i^{\prime})$, for any
$i,i^{\prime}\in\mathcal{V}^{(\rm MNIST)}$, we construct the edges of
$\mathcal{E}^{(\rm MNIST)}$ and their weights as follows. Each node
$i\in\mathcal{V}^{(\rm MNIST)}$ is connected with the four other nodes
$i^{\prime}\in\mathcal{V}^{(\rm MNIST)}\setminus\\{i\\}$ of minimum distance
${\rm dist}(i,i^{\prime})$ by an edge $e=\\{i,i^{\prime}\\}$. The edge weights
are then constructed by an exponential kernel [15], $A_{i,i^{\prime}}={\rm
exp}\big{(}-{\rm dist}(i,i^{\prime})\big{)}$.

To learn the local parameters $\widehat{\mathbf{w}}^{(i)}$ of (40), we use
Algorithm 1 with local loss

$\displaystyle L_{i}\left(\mathbf{w}^{(i)}\right)$
$\displaystyle:=-\big{(}1/m^{(\rm train)}_{i}\big{)}\sum_{\left({\bf
x},y\right)\in\mathcal{X}_{\rm train}^{(i)}}\big{(}y\log{h^{(i)}\big{(}{\bf
x}\big{)}}\big{)}$ $\displaystyle-$ $\displaystyle\big{(}1/m^{(\rm
train)}_{i}\big{)}\sum_{\left({\bf x},y\right)\in\mathcal{X}_{\rm
train}^{(i)}}\big{(}(1\\!-\\!y)\log{\big{(}1-h^{(i)}\big{(}{\bf
x}\big{)}\big{)}}\big{)}.$

As stopping criterion in Algorithm 1, we use a fixed number of $R=50$
iterations. The GTV parameter has been set to $\lambda=1.0$. We measure the
quality of the model parameters learnt by Algorithm 1 via the accuracy ${\rm
ACC}^{(\rm val)}$ achieved on the validation sets $\mathcal{X}_{\rm
val}^{(i)}$. More precisely, ${\rm ACC}^{(\rm val)}$ is the fraction of
correctly classified images in the validation sets $\mathcal{X}_{\rm
val}^{(i)}$, for $i\in\mathcal{V}^{(\rm MNIST)}$.

Figure 4 shows the validation set accuracy achieved by the local model
parameters learnt by Algorithm 1 for different choices of $\lambda$. For the
extreme case $\lambda=0$, Algorithm 1 ignores the edges in $\mathcal{G}^{(\rm
MNIST)}$ and separately minimizes the local loss functions. The other extreme
case is $\lambda\rightarrow\infty$ where all local model parameter are
enforced to be identical, which is equivalent to learning a single model on
all pooled local datasets. Figure 5 compares the ${\rm ACC}^{(\rm val)}$ and
the convergence rate of the Algorithm 1 with existing methods for personalized
FL and for clustered FL. Both the ${\rm ACC}^{(\rm val)}$ and the convergence
rate of the MNIST dataset is improved by Algorithm 1.

$0$$10$$20$$30$$40$$50$$0.8$$0.85$$0.9$$0.95$$k$${\rm ACC}^{(\rm
val)}$$\phi(\mathbf{u})\\!=\\!\|\mathbf{u}\|_{2}$
(nLasso)$\lambda\\!=\\!0$$\lambda\\!=\\!1$$\lambda\\!=\\!10^{3}$ Figure 4: The
validation accuracy ${\rm ACC}^{(\rm val)}$ as a function ot the number $k$ of
iterations run by Algorithm 1. Different curves correspond to difference
choices for the GTV parameter $\lambda$.
$0$$10$$20$$30$$40$$50$$0.5$$0.6$$0.7$$0.8$$0.9$$1$$k$${\rm ACC}^{(\rm
val)}$$\lambda=1,\phi(\mathbf{u})\\!=\\!\|\mathbf{u}\|_{2}$ (nLasso)Algorithm
1IFCA [12]pFedMe [52]PyTorch optimizer Figure 5: Comparing the validation set
accuracy ${\rm ACC}^{(\rm val)}$ achieved by local model parameters learnt by
Algorithm 1 with those learnt by IFCA [12], pFedMe [52] and the PyTorch
optimizer for (16).

##  VII Conclusion

We have studied GTV minimization methods for the distributed learning of
personalized models in networked collections of local datasets. GTV
minimization is an instance of the RERM principle that is obtained using the
GTV of local model parameters as the regularization term. This approach is
built on the assumption that the statistical properties or similarities
between collections of local datasets is reflected by a known network
structure. We obtain highly-scalable NFL methods by solving GTV minimization
using distributed optimization methods that are obtained from jointly solving
GTV minimization with a dual network flow optimization problem. The resulting
message-passing algorithm exploits the known network structure of data to
adaptively pool local datasets into clusters. This pooling allows to learn
“personalized” models for local datasets that would not provide sufficient
statistical power by itself alone. Future research directions include the
joint learning of network structure and networked model weights. We also aim
at a more fine-grained convergence analysis of distributed optimization
methods for NFL that takes into account the precise cluster structure of the
data network. Another interesting avenue for future research is the extension
of our approach to allow for non-parametric local models such as decision
trees.

##  VIII Acknowledgement

We are grateful to Reza Mirzaeifard who pointed out the Moreau decomposition
property of proximity operators and provided a careful review of the
manuscript draft.

## References

  * [1] S. Cui, A. Hero, Z.-Q. Luo, and J. Moura, Eds., _Big Data over Networks_. Cambridge Univ. Press, 2016\.
  * [2] M. Wollschlaeger, T. Sauter, and J. Jasperneite, “The future of industrial communication: Automation networks in the era of the internet of things and industry 4.0,” _IEEE Industrial Electronics Magazine_ , vol. 11, no. 1, pp. 17–27, 2017.
  * [3] H. Ates, A. Yetisen, F. Güder, and C. Dincer, “Wearable devices for the detection of covid-19,” _Nature Electronics_ , vol. 4, no. 1, pp. 13–14, 2021.
  * [4] M. E. J. Newman, _Networks: An Introduction_. Oxford Univ. Press, 2010.
  * [5] A. Barabási, N. Gulbahce, and J. Loscalzo, “Network medicine: a network-based approach to human disease,” _Nature Reviews Genetics_ , vol. 12, no. 56, 2011.
  * [6] K. Grantz, H. Meredith, D. Cummings, C. Metcalf, B. Grenfell, J. Giles, S. Mehta, S. Solomon, A. Labrique, N. Kishore, C. Buckee, and A. Wesolowski, “The use of mobile phone data to inform analysis of COVID-19 pandemic epidemiology,” _Nature Communications_ , vol. 11, no. 1, p. 4961, 2020.
  * [7] B. McMahan, E. Moore, D. Ramage, S. Hampson, and B. A. y Arcas, “Communication-efficient learning of deep networks from decentralized data,” in _Proceedings of the 20th International Conference on Artificial Intelligence and Statistics_ , ser. Proceedings of Machine Learning Research, A. Singh and J. Zhu, Eds., vol. 54. Fort Lauderdale, FL, USA: PMLR, 20–22 Apr 2017, pp. 1273–1282.
  * [8] Y. Cheng, Y. Liu, T. Chen, and Q. Yang, “Federated learning for privacy-preserving ai,” _Communications of the ACM_ , vol. 63, no. 12, pp. 33–36, Dec. 2020.
  * [9] V. Smith, C.-K. Chiang, M. Sanjabi, and A. Talwalkar, “Federated multi-task learning,” in _Advances in Neural Information Processing Systems_ , vol. 30, 2017.
  * [10] J. You, J. Wu, X. Jin, and M. Chowdhury, “Ship compute or ship data? why not both?” in _18th USENIX Symposium on Networked Systems Design and Implementation (NSDI 21)_. USENIX Association, April 2021, pp. 633–651.
  * [11] F. Sattler, K. Müller, T. Wiegand, and W. Samek, “On the byzantine robustness of clustered federated learning,” in _2020 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP)_ , Barcelona, Spain, 2020, pp. 8861–8865.
  * [12] A. Ghosh, J. Chung, D. Yin, and K. Ramchandran, “An efficient framework for clustered federated learning,” in _34th Conference on Neural Information Processing Systems (NeurIPS 2020)_ , Vancouver, Canada, 2020.
  * [13] L. Jacob, J.-P. Vert, and F. Bach, “Clustered multi-task learning: A convex formulation,” in _Advances in Neural Information Processing Systems_ , D. Koller, D. Schuurmans, Y. Bengio, and L. Bottou, Eds., vol. 21. Curran Associates, Inc., 2009.
  * [14] O. Chapelle, B. Schölkopf, and A. Zien, Eds., _Semi-Supervised Learning_. Cambridge, Massachusetts: The MIT Press, 2006.
  * [15] U. von Luxburg, “A tutorial on spectral clustering,” _Statistics and Computing_ , vol. 17, no. 4, pp. 395–416, Dec. 2007.
  * [16] D. P. Bertsekas, _Network Optimization: Continuous and Discrete Models_. Athena Scientific, 1998.
  * [17] R. T. Rockafellar, _Network Flows and Monotropic Optimization_. Athena Scientific, Jul. 1998.
  * [18] A. Jung, “Learning the conditional independence structure of stationary time series: A multitask learning approach,” _IEEE Trans. Signal Processing_ , vol. 63, no. 21, Nov. 2015.
  * [19] D. Hallac, J. Leskovec, and S. Boyd, “Network lasso: Clustering and optimization in large graphs,” in _Proc. SIGKDD_ , 2015, pp. 387–396.
  * [20] A. Jung and N. Tran, “Localized linear regression in networked data,” _IEEE Sig. Proc. Lett._ , vol. 26, no. 7, Jul. 2019.
  * [21] R. Nassif, S. Vlaski, C. Richard, J. Chen, and A. Sayed, “Multitask learning over graphs: An approach for distributed, streaming machine learning,” _IEEE Signal Processing Magazine_ , vol. 37, no. 3, pp. 14–25, 2020.
  * [22] Z. Xu and K. Kersting, “Multi-task learning with task relations,” in _2011 IEEE 11th International Conference on Data Mining_ , 2011, pp. 884–893.
  * [23] L. Condat, “A primal–dual splitting method for convex optimization involving lipschitzian, proximable and linear composite terms,” _Journal of Opt. Th. and App._ , vol. 158, no. 2, pp. 460–479, Aug. 2013.
  * [24] T. Hastie, R. Tibshirani, and M. Wainwright, _Statistical Learning with Sparsity. The Lasso and its Generalizations_. CRC Press, 2015.
  * [25] I. Goodfellow, Y. Bengio, and A. Courville, _Deep Learning_. MIT Press, 2016.
  * [26] S. Boyd, N. Parikh, E. Chu, B. Peleato, and J. Eckstein, _Distributed Optimization and Statistical Learning via the Alternating Direction Method of Multipliers_. Hanover, MA: Now Publishers, 2010, vol. 3, no. 1.
  * [27] A. Nedic and A. Ozdaglar, “Distributed subgradient methods for multi-agent optimization,” _IEEE Transactions on Automatic Control_ , vol. 54, no. 1, pp. 48–61, 2009.
  * [28] K. Scaman, F. Bach, S. Bubeck, L. Massoulié, and Y. T. Lee, “Optimal algorithms for non-smooth distributed optimization in networks,” in _Advances in Neural Information Processing Systems_ , S. Bengio, H. Wallach, H. Larochelle, K. Grauman, N. Cesa-Bianchi, and R. Garnett, Eds., vol. 31. Curran Associates, Inc., 2018\.
  * [29] A. Chambolle and T. Pock, “An introduction to continuous optimization for imaging,” _Acta Numer._ , vol. 25, pp. 161–319, 2016.
  * [30] J. S. Yedidia, “Message-passing algorithms for inference and optimization,” _Journal of Statistical Physics_ , vol. 145, no. 4, pp. 860–890, 2011.
  * [31] A. Chambolle, M. Ehrhardt, P. Richtárik, and C.-B. Schönlieb, “Stochastic primal-dual hybrid gradient algorithm with arbitrary sampling and imaging applications,” _SIAM Journal on Optimization_ , vol. 28, no. 4, pp. 2783–2808, 2018.
  * [32] J. Rasch and A. Chambolle, “Inexact first-order primal–dual algorithms,” _Computational Optimization and Applications_ , vol. 76, no. 2, pp. 381–430, 2020.
  * [33] P. Bühlmann and S. van de Geer, _Statistics for High-Dimensional Data_. New York: Springer, 2011.
  * [34] D. Sun, K.-C. Toh, and Y. Yuan, “Convex clustering: Model, theoretical guarantee and efficient algorithm,” _Journal of Machine Learning Research_ , vol. 22, no. 9, pp. 1–32, 2021.
  * [35] S. Boyd and L. Vandenberghe, _Convex Optimization_. Cambridge, UK: Cambridge Univ. Press, 2004.
  * [36] A. Mahmoudi, H. Shokri-Ghadikolaei, and C. Fischione, “Cost-efficient distributed optimization in machine learning over wireless networks.” organization, 2020.
  * [37] P. J. Brockwell and R. A. Davis, _Time Series: Theory and Methods_. Springer New York, 1991.
  * [38] A. Jung, _Machine Learning: The Basics_ , 1st ed. Springer Singapore, Feb. 2022.
  * [39] T. Hastie, R. Tibshirani, and J. Friedman, _The Elements of Statistical Learning_ , ser. Springer Series in Statistics. New York, NY, USA: Springer, 2001.
  * [40] D. P. Bertsekas, A. Nedic, and A. E. Ozdaglar, _Convex Analysis and Optimization_. Athena Scientific, 2003\.
  * [41] Y. SarcheshmehPour, M. Leinonen, and A. Jung, “Federated learning from big data over networks,” in _Proc. IEEE Int. Conf. on Acoustics, Speech and Signal Processing (ICASSP), preprint: https://arxiv.org/pdf/2010.14159.pdf_ , 2021\.
  * [42] A. Jung, N. Tran, and A. Mara, “When is Network Lasso Accurate?” _Front. Appl. Math. Stat._ , vol. 3, Jan. 2018.
  * [43] N. Parikh and S. Boyd, “Proximal algorithms,” _Foundations and Trends in Optimization_ , vol. 1, no. 3, pp. 123–231, 2013.
  * [44] R. T. Rockafellar, _Convex Analysis_. Princeton, NJ: Princeton Univ. Press, 1970.
  * [45] G. Golub and C. van Loan, “An analysis of the total least squares problem,” _SIAM J. Numerical Analysis_ , vol. 17, no. 6, pp. 883–893, Dec. 1980.
  * [46] C. G. Atkeson, A. W. Moore, and S. Schaal, “Locally weighted learning,” _Artif. Intell. Rev._ , vol. 11, no. 1–5, pp. 11–73, feb 1997.
  * [47] E. Abbe, “Community detection and stochastic block models: Recent developments,” _Journal of Machine Learning Research_ , vol. 18, no. 177, pp. 1–86, 2018.
  * [48] T. Sun, D. Li, and B. Wang, “Decentralized federated averaging,” _ArXiv_ , vol. abs/2104.11375, 2021.
  * [49] Y. Lecun, L. Bottou, Y. Bengio, and P. Haffner, “Gradient-based learning applied to document recognition,” in _Proceedings of the IEEE_ , vol. 86, no. 11, Nov. 1998, pp. 2278–2324.
  * [50] L. van der Maaten and G. Hinton, “Visualizing data using t-sne,” _Jour. Mach. Learning Res._ , vol. 9, pp. 2579–2605, 2008.
  * [51] F. Perez-Cruz, “Kullback-leibler divergence estimation of continuous distributions,” in _IEEE Int. Symp. Inf. Theory_ , 2008, pp. 1666–1670.
  * [52] C. T Dinh, N. Tran, and J. Nguyen, “Personalized federated learning with moreau envelopes,” _Advances in Neural Information Processing Systems_ , vol. 33, pp. 21 394–21 405, 2020.
  * [53] D. Liu and J. Nocedal, “On the limited memory method for large scale optimization,” _Mathematical Programming B_ , vol. 5, no. 3, pp. 503–528, 1989.
  * [54] T. Pock and A. Chambolle, “Diagonal preconditioning for first order primal-dual algorithms in convex optimization,” in _IEEE ICCV_ , Barcelona, Spain, Nov. 2011.
  * [55] H. H. Bauschke and P. L. Combettes, _Convex Analysis and Monotone Operator Theory in Hilbert Spaces_. New York: Springer, 2011.
  * [56] F. M. Institute. (2022) Download observations. [Online]. Available: https://en.ilmatieteenlaitos.fi/download-observations
  * [57] C. R. Givens and R. M. Shortt, “A class of Wasserstein metrics for probability distributions.” _Michigan Mathematical Journal_ , vol. 31, no. 2, pp. 231 – 240, 1984.

Supplementary Material: Networked Federated Learning

##  IX Supplementary Material

###  IX-A Proof of Theorem 2

Our proof technique is an instance of the primal-dual witness technique which
has been championed for the analysis of Lasso methods [24, Sec. 11.4.]. A
similar instance of this technique has been applied to the study of convex
clustering which is a special case of GTV minimization [34, Theorem 5], In
particular, the proof revolves around constructing a pair of optimal networked
model parameters $\widehat{{\bf w}}\in\mathcal{W}$ and optimal flow vector
$\widehat{\mathbf{u}}\in\mathcal{U}$ that satisfies the optimality condition
(III-B). The networked model parameters $\widehat{{\bf w}}\in\mathcal{W}$
assign the local model parameters $\widehat{\mathbf{w}}^{(i)}$ to each node
$i\in\mathcal{V}$. The optimal flow vector
$\widehat{\mathbf{u}}\in\mathcal{U}$ that assigns each edge the flow vector
$\widehat{\mathbf{u}}^{(e)}$, for each edge $e\in\mathcal{E}$.

We construct the primal-dual optimal pair for the GTV minimization (16) over
the empirical graph $\mathcal{G}$ using the solutions of another GTV
minimization problem that is defined for another empirical graph
$\widetilde{\mathcal{G}}=\left(\widetilde{\mathcal{V}},\widetilde{\mathcal{E}}\right)$.
This other empirical graph $\widetilde{\mathcal{G}}$, which we will refer to
as “cluster graph”, is obtained from the original empirical graph
$\mathcal{G}$ by merging all nodes that belong to the same cluster according
to Assumption 1 (see Figure 6).

The nodes $\widetilde{\mathcal{V}}=\\{1,\ldots,k\\}$ of the cluster graph
$\widetilde{\mathcal{G}}$ represent the clusters of the partition (10) in
Assumption 1. Two different nodes $c,c^{\prime}\in\widetilde{\mathcal{V}}$ are
connected by an edges $\\{c,c^{\prime}\\}\in\widetilde{\mathcal{E}}$ if and
only if the empirical graph $\mathcal{G}$ contains two different nodes
$i\in\mathcal{C}^{(c)},i^{\prime}\in\mathcal{C}^{(c^{\prime})}$ that are
connected by an edge $\\{i,i^{\prime}\\}\in\mathcal{E}$.

Each cluster node $c$ carries the sum $f^{(c)}\left({\bf
v}\right):=\sum_{i\in\mathcal{C}^{(c)}}L_{i}\left({\bf v}\right)$ of the local
loss function assigned to the nodes $i\in\mathcal{C}^{(c)}$ of the original
graph $\mathcal{G}$ (see our clustering assumption in (9)). Two cluster nodes
$c\neq c^{\prime}$ are connected by an undirected edge $\\{c,c^{\prime}\\}$
with weight

$\widetilde{A}_{c,c^{\prime}}=\sum_{i\in\mathcal{C}^{(c)},i^{\prime}\in\mathcal{C}^{(c^{\prime})}}A_{i,i^{\prime}}.\vspace{-2mm}$
(41)

Thus, each edge in the cluster graph $\widetilde{\mathcal{G}}$ is obtained by
aggregating all edges (via summing their weights) in the original empirical
graph $\mathcal{G}$ that connect nodes in cluster $\mathcal{C}^{(c)}$ with
nodes in cluster $\mathcal{C}^{(c^{\prime})}$. Figure 6 illustrates an example
for an empirical graph whose nodes are partitioned into clusters along with
the corresponding cluster graph.

The original GTV minimization (16) over the empirical graph $\mathcal{G}$
induces the following GTV minimization problem on the cluster graph
$\widetilde{\mathcal{G}}$:

$\widetilde{\mathbf{w}}\in\underset{\mathbf{w}\in\widetilde{\mathcal{W}}}{\mathrm{arg\
min}}\sum_{c\in\widetilde{\mathcal{V}}}f^{(c)}\left(\mathbf{w}^{(c)}\right)+\lambda\sum_{e\in\widetilde{\mathcal{E}}}\widetilde{A}_{e}\|{\mathbf{w}^{(e_{+})}-\mathbf{w}^{(e_{-})}}\|.$
(42)

Here, $\widetilde{\mathcal{W}}$ denotes the space of maps
$\mathbf{w}^{(c)}:\widetilde{\mathcal{V}}\rightarrow\mathbb{R}^{d}$. Such a
map assign each node $c\in\widetilde{\mathcal{V}}$ in the cluster graph a
parameter vector $\mathbf{w}^{(c)}\in\mathbb{R}^{d}$. With slight abuse of
notation, we refer by $\mathbf{w}^{(c)}$ to both, the entire map (which is an
element of $\widetilde{\mathcal{W}}$), and the vector that is assigned to a
particular node $c\in\widetilde{\mathcal{V}}$ by this map. Beside the space
$\widetilde{\mathcal{W}}$, we associate the space $\widetilde{\mathcal{U}}$
with the cluster graph $\widetilde{\mathcal{G}}$. This space consists of all
maps $\mathbf{u}$ that assign each edge
$e=\\{c,c^{\prime}\\}\in\widetilde{\mathcal{E}}$ a vector
$\mathbf{u}^{(e)}\in\mathbb{R}^{d}$.

According to the optimality condition (III-B), any solution
$\widetilde{\mathbf{w}}$ of (42) is characterized by

$\displaystyle\sum_{e\in\widetilde{\mathcal{E}}:e_{+}=c}\widetilde{{\bf
u}}^{(e)}-\sum_{e\in\widetilde{\mathcal{E}}:e_{-}=c}\widetilde{{\bf
u}}^{(e)}=-\nabla f^{(c)}\left(\widetilde{{\bf w}}^{(c)}\right)$
$\displaystyle\quad\quad\quad\quad\mbox{ for all nodes
}c\in\widetilde{\mathcal{V}},$ (43) $\displaystyle\left\|{\widetilde{{\bf
u}}^{(e)}}\right\|_{*}\leq\lambda\widetilde{A}_{e}\mbox{ for every edge
}e\in\widetilde{\mathcal{E}},$ (44) $\displaystyle\widetilde{{\bf
w}}^{(e_{+})}$ $\displaystyle=\widetilde{{\bf w}}^{(e_{-})}\mbox{ for every
edge }e\in\widetilde{\mathcal{E}}\mbox{ with }\left\|{\widetilde{{\bf
u}}^{(e)}}\right\|_{*}<\lambda A_{e}.$ (45)

Note that (43) - (45) are necessary and sufficient for a pair $\widetilde{{\bf
w}}^{(c)}$, $\widetilde{{\bf u}}^{(e)}$ to solve (42) and its dual (which is
obtained by adapting (24) to the cluster graph $\widetilde{\mathcal{G}}$),
respectively. By combining (43) with (44), we obtain the necessary condition

$\left\|{\nabla f^{(c)}\left(\widetilde{{\bf
w}}^{(c)}\right)}\right\|_{*}\leq\lambda\sum_{c^{\prime}}\widetilde{A}_{c,c^{\prime}},\vspace{-2mm}$
(46)

which holds for any solution $\widetilde{{\bf w}}^{(c)}$ of (42).

The upper bound (46) on the (norm of the) gradient $\nabla
f^{(c)}\left(\widetilde{{\bf w}}^{(c)}\right)$ allows us to upper bound the
deviation between ${\widetilde{{\bf w}}^{(c)}}$ and the cluster-wise optimizer
(9). To this end, we use the convexity of $f^{(c)}\left(\widetilde{{\bf
w}}^{(c)}\right)$ (see Assumption 2) to obtain

$\displaystyle f^{(c)}\left(\widetilde{{\bf w}}^{(c)}\right)-$ $\displaystyle
f^{(c)}\left(\overline{{\bf w}}^{(c)}\right)\leq\big{(}\widetilde{{\bf
w}}^{(c)}-\overline{{\bf w}}^{(c)}\big{)}^{T}\nabla
f^{(c)}\left(\widetilde{{\bf w}}^{(c)}\right)$
$\displaystyle\stackrel{{\scriptstyle(a)}}{{\leq}}\left\|{\widetilde{{\bf
w}}^{(c)}-\overline{{\bf w}}^{(c)}}\right\|\left\|{\nabla
f^{(c)}\left(\widetilde{{\bf w}}^{(c)}\right)}\right\|_{*}$
$\displaystyle\stackrel{{\scriptstyle\eqref{equ_dondition_gradient_clustered_graph}}}{{\leq}}\left\|{\widetilde{{\bf
w}}^{(c)}-\overline{{\bf
w}}^{(c)}}\right\|\lambda\sum_{c^{\prime}\in\widetilde{\mathcal{V}}\setminus\\{c\\}}\widetilde{A}_{c,c^{\prime}}.$
(47)

Here, step $(a)$ uses the definition of the dual norm [35]. Next, we use the
strong convexity of $f^{(c)}\left(\cdot\right)$ (see (13)) to lower bound the
LHS of (47) as

$\displaystyle f^{(c)}\left(\widetilde{{\bf
w}}^{(c)}\right)-f^{(c)}\left(\overline{{\bf
w}}^{(c)}\right)\geq(\alpha^{(c)}/2)\left\|{\widetilde{{\bf
w}}^{(c)}-\overline{{\bf w}}^{(c)}}\right\|^{2}.$ (48)

Combining (48) with (47) yields

$\left\|{\widetilde{{\bf w}}^{(c)}-\overline{{\bf
w}}^{(c)}}\right\|\leq(2\lambda/\alpha^{(c)})\sum_{c^{\prime}\in\widetilde{\mathcal{V}}\setminus\\{c\\}}\widetilde{A}_{c,c^{\prime}}.\vspace{-2mm}$
(49)

To recap: We have introduced the cluster graph $\widetilde{\mathcal{G}}$ by
merging the nodes of the empirical graph $\mathcal{G}$ according to the
clustering assumption (see Assumption 1). This cluster graph is used to define
the “cluster-wise” GTV minimization problem (42) which is obtained from the
original GTV minimization (16) by enforcing identical local model parameters
for all nodes in the same cluster of $\mathcal{G}$. We have then derived the
upper bound (49) on the deviation between any solution $\widetilde{{\bf
w}}^{(c)}$ of (42) and the cluster-wise optimizers (9).

The main step in the proof of Theorem 2 is to show that the local model
parameters

$\displaystyle\widehat{\mathbf{w}}^{(i)}$ $\displaystyle:=\widetilde{{\bf
w}}^{(c)}$ (50) $\displaystyle\mbox{ for all nodes }i\in\mathcal{C}^{(c)}$
$\displaystyle\mbox{ in the same cluster }\mathcal{C}^{(c)}\in\mathcal{P}$
(see Assumption 1)

solve the original GTV minimization problem (16). To this end, we construct a
flow $\widehat{\mathbf{u}}^{(e)}$, for each edge $e\in\mathcal{E}$ in the
original empirical graph $\mathcal{G}$, that together with
$\widehat{\mathbf{w}}^{(i)}$ defined by (50), solves the optimality conditions
(III-B).

We construct the optimal flow $\widehat{\mathbf{u}}^{(e)}$ by considering
separately the inter-cluster edges (connecting nodes in different clusters)
and the intra-cluster edges (connecting nodes in the same cluster). The inter-
cluster (or boundary) edges are

$\partial\mathcal{P}:=\left\\{\\{i,i^{\prime}\\}\in\mathcal{E},\mbox{ with
}i\in\mathcal{C}^{(c)},i^{\prime}\in\mathcal{C}^{(c^{\prime})}\right\\}.\vspace{-2mm}$

The intra-cluster edges $\mathcal{E}\setminus\partial\mathcal{P}$ connect
nodes in the same cluster, i.e., they are given by all edges
$\\{i,i^{\prime}\\}$ with $i,i^{\prime}\in\mathcal{C}^{(c)}$ for some
$c\in\\{1,\ldots,k\\}$.

To each inter-cluster edge $e=\\{i,i^{\prime}\\}\in\partial\mathcal{P}$, we
assign the flow value

$\widehat{\mathbf{u}}^{(e)}:=\widetilde{{\bf
u}}^{(\\{c,c^{\prime}\\})}\big{(}A_{i,i^{\prime}}/\widetilde{A}_{c,c^{\prime}}\big{)}\mbox{
(see \eqref{equ_def_edge_weight_clustered})}.\vspace{-2mm}$ (51)

Note that, for any two different cluster nodes $c,c^{\prime}$ (representing
the clusters $\mathcal{C}^{(c)}$ and $\mathcal{C}^{(c^{\prime})}$ in the
original empirical graph $\mathcal{G}$),

$\displaystyle\left\|{\widehat{\mathbf{u}}^{(e)}}\right\|_{*}$
$\displaystyle=\left\|{\widetilde{{\bf
u}}^{(\\{c,c^{\prime}\\})}}\right\|_{*}\big{(}A_{i,i^{\prime}}/\widetilde{A}_{c,c^{\prime}}\big{)}\stackrel{{\scriptstyle\eqref{equ_op_condition_clusteredg_capacity_constraint}}}{{\leq}}\lambda
A_{i,i^{\prime}}$ $\displaystyle\mbox{ for any }e=\\{i,i^{\prime}\\}\mbox{
with
}i\in\mathcal{C}^{(c)},i^{\prime}\in\mathcal{C}^{(c^{\prime})}.\vspace{-2mm}$
(52)

Moreover, if the cluster-wise GTV minimization (42) delivers different
parameter vectors $\widetilde{{\bf w}}^{(c)}\neq\widetilde{{\bf
w}}^{(c^{\prime})}$ for two cluster nodes
$c,c^{\prime}\in\widetilde{\mathcal{V}}$, then (45) requires that the edge
$\\{c,c^{\prime}\\}$ must be saturated in the sense of
$\left\|{\widetilde{{\bf
u}}^{(\\{c,c^{\prime}\\})}}\right\|_{*}=\lambda\widetilde{A}_{c,c^{\prime}}$.
This implies, in turn via (51), that

$\displaystyle\left\|{\widehat{\mathbf{u}}^{(\\{i,i^{\prime}\\})}}\right\|_{*}=\lambda
A_{i,i^{\prime}}$ (53) $\displaystyle\mbox{ for any }i\in\mathcal{C}^{(c)},$
$\displaystyle i^{\prime}\in\mathcal{C}^{(c^{\prime})}\mbox{ such that
}\widetilde{{\bf w}}^{(c)}\neq\widetilde{{\bf w}}^{(c^{\prime})}.$

Note that (53) is equivalent to the third condition in (III-B) for the optimal
flow values across inter-cluster edges $e\in\partial\mathcal{P}$.

We next construct the flow values $\widehat{\mathbf{u}}^{(e)}$ for the intra-
cluster edges $e\in\mathcal{E}\setminus\partial\mathcal{P}$ such that the
optimality conditions (III-B) are satisfied. To this end, we consider a
specific cluster $\mathcal{C}^{(c)}\in\mathcal{P}$ and the corresponding
induced sub-graph
$\mathcal{G}^{(c)}=\left(\mathcal{C}^{(c)},\mathcal{E}^{(c)}\right)$ of the
empirical graph $\mathcal{G}$. Together with (IX-A) and (53), a sufficient
condition for (III-B) to hold is that, for each cluster sub-graph
$\mathcal{G}^{(c)}$,

$\displaystyle\sum_{e\in\mathcal{E}^{(c)}:i\\!=\\!e_{+}}\widehat{\mathbf{u}}^{(e)}$
$\displaystyle-\sum_{e\in\mathcal{E}^{(c)}:i\\!=\\!e_{-}}\widehat{\mathbf{u}}^{(e)}=$
$\displaystyle-\sum_{e\in\partial\mathcal{P}:i\\!=\\!e_{+}}\widehat{\mathbf{u}}^{(e)}+\sum_{e\in\partial\mathcal{P}:i\\!=\\!e_{-}}\widehat{\mathbf{u}}^{(e)}-\nabla
L_{i}\left(\widehat{\mathbf{w}}^{(i)}\right)$ $\displaystyle\quad\quad\mbox{
for each node }i\in\mathcal{C}^{(c)}\mbox{, and }$
$\displaystyle\|\widehat{\mathbf{u}}^{(e)}\|_{*}$ $\displaystyle<\lambda
A_{e}\mbox{ for each edge }e\in\mathcal{E}^{(c)}.\vspace{-2mm}$ (54)

Regarding the first (flow conservation) condition in (IX-A), we note that

$\displaystyle\left\|{\nabla
L_{i}\left(\widehat{\mathbf{w}}^{(i)}\right)}\right\|_{*}=$
$\displaystyle\quad\left\|{\nabla
L_{i}\left(\widehat{\mathbf{w}}^{(i)}\right)-\nabla L_{i}\left(\overline{{\bf
w}}^{(c)}\right)+\nabla L_{i}\left(\overline{{\bf
w}}^{(c)}\right)}\right\|_{*}$
$\displaystyle\quad\stackrel{{\scriptstyle(a)}}{{\leq}}\left\|{\nabla
L_{i}\left(\widehat{\mathbf{w}}^{(i)}\right)-\nabla L_{i}\left(\overline{{\bf
w}}^{(c)}\right)}\right\|_{*}+\left\|{\nabla L_{i}\left(\overline{{\bf
w}}^{(c)}\right)}\right\|_{*}$
$\displaystyle\quad\stackrel{{\scriptstyle\eqref{equ_upper_bound_norm_gradient}}}{{\leq}}\left\|{\nabla
L_{i}\left(\widehat{\mathbf{w}}^{(i)}\right)-\nabla L_{i}\left(\overline{{\bf
w}}^{(c)}\right)}\right\|_{*}+\delta^{(i)}$
$\displaystyle\quad\stackrel{{\scriptstyle\eqref{equ_def_cond_locallipsch}}}{{\leq}}\beta^{(i)}\|{\widehat{\mathbf{w}}^{(i)}-\overline{{\bf
w}}^{(c)}}\|+\delta^{(i)}$
$\displaystyle\quad\stackrel{{\scriptstyle\eqref{equ_def_upper_bound_clustered_custerwide}}}{{\leq}}\beta^{(i)}(2\lambda/\alpha^{(c)})|\partial\mathcal{C}^{(c)}|+\delta^{(i)}$
$\displaystyle\quad\mbox{ for each node }i\in\mathcal{C}^{(c)}.$ (55)

Here, step $(a)$ in (IX-A) used the triangle inequality for the norm
$\left\|{\cdot}\right\|_{*}$ [45]. Another application of the triangle
inequality yields, for any node $i\in\mathcal{C}^{(c)}$,

$\displaystyle\left\|{-\sum_{e\in\partial\mathcal{P}:i\\!=\\!e_{+}}\widehat{\mathbf{u}}^{(e)}+\sum_{e\in\partial\mathcal{P}:i\\!=\\!e_{-}}\widehat{\mathbf{u}}^{(e)}}\right\|_{*}$
$\displaystyle\leq\sum_{i^{\prime}\in\mathcal{N}^{(i)}\cap\mathcal{V}\setminus\mathcal{C}^{(c)}}\left\|{\widehat{\mathbf{u}}^{(\\{i,i^{\prime}\\})}}\right\|_{*}$
$\displaystyle\stackrel{{\scriptstyle\eqref{equ_def_cap_constraint_inter_cluster_edge}}}{{\leq}}\sum_{i^{\prime}\in\mathcal{N}^{(i)}\cap\mathcal{V}\setminus\mathcal{C}^{(c)}}\lambda
A_{i,i^{\prime}}$
$\displaystyle\leq\sum_{i^{\prime}\notin\mathcal{C}^{(c)}}\lambda
A_{i,i^{\prime}}.$ (56)

Let us rewrite the requirement (IX-A), for some cluster $\mathcal{C}^{(c)}$,
as

$\displaystyle\sum_{e\in\mathcal{E}_{c}:i\\!=\\!e_{+}}\widehat{\mathbf{u}}^{(e)}-\sum_{e\in\mathcal{E}_{c}:i\\!=\\!e_{-}}\widehat{\mathbf{u}}^{(e)}$
$\displaystyle=\mathbf{s}^{(i)}\mbox{ for each node }i\in\mathcal{C}^{(c)},$
$\displaystyle\mbox{ and }\|\widehat{\mathbf{u}}^{(e)}\|_{*}<\lambda$
$\displaystyle A_{e}\mbox{ for each edge }e\in\mathcal{E}^{(c)}.$ (57)

with vector-valued demands
$\mathbf{s}^{(i)}=-\sum_{e\in\partial\mathcal{P}:i\\!=\\!e_{+}}\widehat{\mathbf{u}}^{(e)}+\sum_{e\in\partial\mathcal{P}:i\\!=\\!e_{-}}\widehat{\mathbf{u}}^{(e)}-\nabla
L_{i}\left(\widehat{\mathbf{w}}^{(i)}\right)$. According to (IX-A) and (IX-A),
we have

$\left\|{\mathbf{s}^{(i)}}\right\|_{*}\leq\sum_{i^{\prime}\notin\mathcal{C}^{(c)}}\lambda
A_{i,i^{\prime}}+\beta^{(i)}(2\lambda/\alpha^{(c)})|\partial\mathcal{C}^{(c)}|+\delta^{(i)}.$

Since Theorem 2 requires each cluster $\mathcal{C}^{(c)}$ to be well-connected
(see (35)), we can apply Lemma 3 to verify the existence of a flow
$\widehat{\mathbf{u}}^{(e)}$ that satisfies (IX-A).

$\mathcal{C}^{(1)}$$\mathcal{C}^{(2)}$$\mathcal{C}^{(3)}$$i^{\prime}$$i$$A_{i,i^{\prime}}$$c\\!=\\!1$$c^{\prime}\\!=\\!3$$\widetilde{A}_{c,c^{\prime}}$
Figure 6: Left: Empirical graph $\mathcal{G}$ of a networked dataset whose
nodes are partitioned into three clusters
$\mathcal{P}=\\{\mathcal{C}^{(1)},\ldots,\mathcal{C}^{(k)}\\}$. Right: Cluster
graph $\widetilde{\mathcal{G}}$ whose nodes represent the clusters in
$\mathcal{P}$.

######  Lemma 3 (Main Helper).

Consider an empirical graph
$\mathcal{G}^{\prime}=\big{(}\mathcal{V}^{\prime},\mathcal{E}^{\prime}\big{)}$
with edge weights $A^{\prime}_{e}$. We orient the edges
$e=\\{i,i^{\prime}\\}\in\mathcal{E}^{\prime}$ by defining its head and tail as
$e_{+}=\min\\{i,i^{\prime}\\}$ and $e_{-}=\max\\{i,i^{\prime}\\}$,
respectively. Each node $i\in\mathcal{V}^{\prime}=\\{1,\ldots,n^{\prime}\\}$
carries a (“demand”) vector $\mathbf{s}^{(i)}\in\mathbb{R}^{d}$ such that
$\sum_{i\in\mathcal{V}^{\prime}}\mathbf{s}^{(i)}=\mathbf{0}$ and
$\left\|{\mathbf{s}^{(i)}}\right\|\leq\delta^{(i)}$ with prescribed upper
bounds $\delta^{(i)}\in\mathbb{R}_{+}$.

If there is a node $i_{0}\in\mathcal{V}^{\prime}$ such that

$\sum_{i\in\mathcal{A}}\delta^{(i)}<\sum_{i\in\mathcal{A}}\sum_{i^{\prime}\in\mathcal{V}^{\prime}\setminus\mathcal{A}}A^{\prime}_{i,i^{\prime}}\mbox{
for all subsets
}\mathcal{A}\subseteq\mathcal{V}^{\prime}\setminus\big{\\{}i_{0}\big{\\}},$
(58)

we can construct a flow $\mathbf{u}^{(e)}\in\mathbb{R}^{d}$, for each edge
$e$, such that

$\displaystyle\sum_{e\in\mathcal{E}^{\prime}:i\\!=\\!e_{+}}\mathbf{u}^{(e)}-\sum_{e\in\mathcal{E}^{\prime}:i\\!=\\!e_{-}}\mathbf{u}^{(e)}$
$\displaystyle=\mathbf{s}^{(i)}\mbox{ for all nodes
}i\in\mathcal{V}^{\prime},$
$\displaystyle\left\|{\mathbf{u}^{(e)}}\right\|<A^{\prime}_{e}$
$\displaystyle\mbox{ for every edge }e\in\mathcal{E}^{\prime}.$ (59)

###### Proof.

We first apply the _Feasible Distribution Theorem_ [17, p. 71] to (58) which
yields the existence of a scalar-valued flow $u^{(e)}$ such that

$\displaystyle\sum_{e\in\mathcal{E}^{\prime}:i_{0}\\!=\\!e_{+}}u^{(e)}-\sum_{e\in\mathcal{E}^{\prime}:i_{0}\\!=\\!e_{-}}u^{(e)}$
$\displaystyle=-\sum_{i\in\mathcal{V}^{\prime}\setminus\\{i_{0}\\}}\left\|{\mathbf{s}^{(i)}}\right\|,$
$\displaystyle\sum_{e\in\mathcal{E}^{\prime}:i\\!=\\!e_{+}}u^{(e)}-\sum_{e\in\mathcal{E}^{\prime}:i\\!=\\!e_{-}}u^{(e)}$
$\displaystyle=\left\|{\mathbf{s}^{(i)}}\right\|$ $\displaystyle\mbox{ for all
nodes }i\in\mathcal{V}^{\prime}$
$\displaystyle\setminus\left\\{i_{0}\right\\}\mbox{, and }$
$\displaystyle|u^{(e)}|<A^{\prime}_{e}\mbox{ for every}$ $\displaystyle\mbox{
edge }e\in\mathcal{E}^{\prime}.$ (60)

Next, we decompose the flow $u^{(e)}$ using the _Conformal Decomposition
Theorem_ [16, Prop 1.1] into elementary (scalar-valued) flows
$u^{(e)}_{r}\in\\{0,1\\}$ with the index $r$ representing a simple path,

$u^{(e)}=\sum_{r=1}^{n^{\prime}+|\mathcal{E}^{\prime}|}a_{r}u^{(e)}_{r}.$ (61)

Each elementary flow $u^{(e)}_{r}$ in (61) is associated with separate path
$\mathcal{P}^{(r)}$ that starts from a node in
$\mathcal{V}\setminus\\{i_{0}\\}$ and terminates at the node $i_{0}$.
Moreover, each flow $u^{(e)}_{r}$ in (61) conforms to the flow $u^{(e)}$ in
the sense of [16, p.8]

$\displaystyle u^{(e)}$ $\displaystyle>0\mbox{ for any edge
}e\in\mathcal{E}^{\prime}\mbox{ with }u^{(e)}_{r}>0\mbox{, and }$
$\displaystyle u^{(e)}$ $\displaystyle<0\mbox{ for any edge
}e\in\mathcal{E}^{\prime}\mbox{ with }u^{(e)}_{r}<0.$ (62)

Let the set $\mathcal{S}^{(i)}$ collect the indices
$r\in\\{1,\ldots,n+|\mathcal{E}|\\}$ of all paths $\mathcal{P}^{(r)}$ starting
at node $i\in\mathcal{V}^{\prime}\setminus\\{i_{0}\\}$. The final step of the
proof is to multiply the elementary flows $u^{(e)}_{r}$, for each path
$r\in\mathcal{S}^{(i)}$, with the vector-valued demand $\mathbf{s}^{(i)}$,
resulting in the vector-valued flow

$\mathbf{u}^{(e)}=\sum_{i\in\mathcal{V}\setminus\\{i_{0}\\}}\mathbf{s}^{(i)}\sum_{r\in\mathcal{S}^{(i)}}u^{(e)}_{r}.$
(63)

∎

###  IX-B Primal Updates for Specific Local Models

Networked Linear Regression. We can use Algorithm 1 to learn a separate linear
predictor $\hat{y}=h({\bf x};{\bf w}^{(i)})=\big{(}{\bf
w}^{(i)}\big{)}^{T}{\bf x}$, for each $i\in\mathcal{V}$. For a local dataset
$\mathcal{X}^{(i)}$, we measure the quality of the linear predictor $h({\bf
x};\mathbf{w}^{(i)})$ via the squared error loss

$L_{i}\left(\mathcal{X}^{(i)},{\bf
v}\right)\\!:=\\!(1/m_{i})\sum_{r=1}^{m_{i}}\big{(}y_{r}^{(i)}\\!-\\!({\bf
x}^{(i,r)})^{T}{\bf v}\big{)}^{2}$ (64)

Inserting (64) into (32) yields the primal update operator

$\mathcal{PU}^{(i)}\left\\{{\bf
v}\right\\}=\left[\mathbf{I}\\!+\\!2\tau_{i}\widetilde{\bf
Q}^{(i)}\right]^{-1}\left[{\bf v}+2\tau_{i}\tilde{{\bf y}}^{(i)}\right].$ (65)

Here, we used $\widetilde{\mathbf{Q}}^{(i)}\\!:=\\!(1/m_{i})\big{(}{\bf
X}^{(i)}\big{)}^{T}{\bf X}^{(i)}$ with ${\bf X}^{(i)}:=\big{(}{\bf
x}^{(i,1)},\ldots,{\bf x}^{(i,m_{i})}\big{)}^{T}$ and the normalized label
vector $\tilde{{\bf y}}^{(i)}\\!:=\\!(1/m_{i})\big{(}{\bf
X}^{(i)}\big{)}^{T}{\bf y}^{(i)}$ with ${\bf
y}^{(i)}:=\big{(}y^{(i)}_{1},\ldots,y_{m_{i}}^{(i)}\big{)}^{T}$.

Networked Logistic Regression. We can use Algorithm 1 to learn linear
classifiers $h^{(i)}\big{(}{\bf x};\mathbf{w}^{(i)}\big{)}={\bf
x}^{T}\mathbf{w}^{(i)}$ for $i\in\mathcal{V}$. To measure the quality of
$h({\bf x};{\bf v})$ at node $i\in\mathcal{V}$, we use the logistic loss given
by (for label values $y\in\\{-1,1\\}$),

$L_{i}\left(\mathcal{X}^{(i)},{\bf v}\right)=(1/m_{i})\sum_{r=1}^{m_{i}}{\rm
log}\left(1\\!+\\!\exp\big{(}-y^{(i)}_{r}{\bf v}^{T}{\bf
x}^{(i,r)}\big{)}\right).$

The corresponding primal update operator is

$\displaystyle\mathcal{PU}^{(i)}\left\\{{\bf
v}\right\\}=\operatorname*{argmin}_{{\bf
z}\in\mathbb{R}^{d}}(1/m_{i})\sum_{r=1}^{m_{i}}$ $\displaystyle{\rm
log}\left(1\\!+\\!\exp\big{(}-y^{(i)}_{r}{\bf z}^{T}{\bf
x}^{(i,r)}\big{)}\right)$ $\displaystyle\\!+\\!(1/2\tau_{i})\|{\bf
z}\\!-\\!{\bf v}\|^{2}.$ (66)

While lacking a closed-form solution, (IX-B) can be solved numerically with
efficient optimization methods such as L-BFGS [53].

Networked Regularized Linear Regression. Consider a local dataset
$\mathcal{X}^{(i)}$ of the form (4), consisting of $m_{i}$ data points, each
characterized by a feature vector ${\bf x}^{(i,r)}\in\mathbb{R}^{d}$ and label
$y^{(i,r)}\in\mathbb{R}$. Many application domains generate high-dimensional
data with a number of datapoints much smaller than number of features,
$m_{i}\ll n$ [33, 24]. Regularized linear regression methods allow to handle
this high-dimensional regime by augmenting the squared error loss (64) with a
regularization term [38, Ch. 7]

$L_{i}\left({\bf w}\right):=\left(y^{(r)}-{\bf w}^{T}{\bf
x}^{(r)}\right)^{2}+\eta\mathcal{R}({\bf w}).$ (67)

Adding $\eta\mathcal{R}({\bf w})$ to the squared error loss aims at avoiding
overfitting and favours parameter vectors ${\bf w}$ yielding a hypothesis that
generalized well to data points not contained in $\mathcal{X}^{(i)}$.

Two well-known special cases are ridge regression, using $\mathcal{R}({\bf
w})=\left\|{{\bf w}}\right\|_{2}^{2}$. and the Lasso, which is uses
$\mathcal{R}({\bf w})=\left\|{{\bf w}}\right\|_{1}$ [24]. Plugging the
regularized loss (67) into (32) yields the primal update operator

$\displaystyle\mathcal{PU}^{(i)}\left\\{{\bf
v}\right\\}=\operatorname*{argmin}_{{\bf z}\in\mathbb{R}^{d}}$
$\displaystyle(1/m_{i})\sum_{r=1}^{m_{i}}\big{(}{\bf z}^{T}{\bf
x}^{(i,r)}\\!-\\!y_{r}^{(i)}\big{)}^{2}\\!$
$\displaystyle+\\!\eta\mathcal{R}({\bf z})\\!+\\!(1/2\tau_{i})\|{\bf
z}\\!-\\!{\bf v}\|^{2}.$ (68)

###  IX-C Dual Updates for Specific Penalty Functions

NLasso Penalty.

$\displaystyle\phi^{*}({\bf v})$ $\displaystyle{:=\sup_{{\bf
z}\in\mathbb{R}^{d}}{\bf v}^{T}{\bf z}-\phi({\bf z})}$
$\displaystyle=\sup_{{\bf z}\in\mathbb{R}^{d}}{\bf v}^{T}{\bf z}-\|{\bf
z}\|_{2}$ $\displaystyle=\begin{cases}0\quad\text{if\ }\|{\bf v}\|_{2}\leq
1\\\ \infty\quad\text{otherwise}\end{cases}$

For inserting IX-C into IV if $\|{\bf v}\|_{2}\leq\lambda A_{e}$ so
$\phi^{*}({\bf z}/(\lambda A_{e}))=0$ and the dual update operator would be

$\displaystyle\mathcal{DU}^{(e)}\big{\\{}{\bf v}\big{\\}}$
$\displaystyle{:=\operatorname*{argmin}_{{\bf z}\in\mathbb{R}^{d}}\lambda
A_{e}\phi^{*}\big{(}{\bf z}/(\lambda A_{e})\big{)}\\!+\\!(1/2\sigma_{e})\|{\bf
v}\\!-\\!{\bf z}\|^{2}}$
$\displaystyle\stackrel{{\scriptstyle\eqref{nlasso_penatly_formula_star}}}{{:=}}\operatorname*{argmin}_{{\bf
z}\in\mathbb{R}^{d}}(1/2\sigma_{e})\|{\bf v}\\!-\\!{\bf z}\|^{2}$
$\displaystyle:={\bf v}.$ (70)

And if $\|{\bf v}\|_{2}>\lambda A_{e}$ then $\|(\lambda A_{e}{\bf z}/\|{\bf
z}\|_{2})\|_{2}\leq\lambda A_{e}$, so $\mathcal{DU}^{(e)}\big{\\{}{\bf
v}\big{\\}}:=\mathcal{T}^{(\lambda A_{e})}\big{\\{}{\bf v}\big{\\}}$.
MOCHA Penalty.

$\displaystyle\phi^{*}({\bf v})$ $\displaystyle:=\sup_{{\bf
z}\in\mathbb{R}^{d}}{\bf v}^{T}{\bf z}-\phi({\bf z})$
$\displaystyle=\sup_{{\bf z}\in\mathbb{R}^{d}}{\bf v}^{T}{\bf
z}-\frac{1}{2}\|{\bf z}\|_{2}^{2}$ $\displaystyle=\frac{1}{2}\|{\bf
v}\|_{2}^{2}$ (71)

Inserting IX-C into IV

$\displaystyle\mathcal{DU}^{(e)}\big{\\{}{\bf v}\big{\\}}$
$\displaystyle{:=\operatorname*{argmin}_{{\bf z}\in\mathbb{R}^{d}}\lambda
A_{e}\phi^{*}\big{(}{\bf z}/(\lambda A_{e})\big{)}\\!+\\!(1/2\sigma_{e})\|{\bf
v}\\!-\\!{\bf z}\|^{2}}$
$\displaystyle\stackrel{{\scriptstyle\eqref{mocha_penatly_formula_star}}}{{:=}}\operatorname*{argmin}_{{\bf
z}\in\mathbb{R}^{d}}\frac{\lambda A_{e}}{2}\|\big{(}{\bf z}/(\lambda
A_{e})\big{)}\|_{2}^{2}\\!+\\!(1/2\sigma_{e})\|{\bf v}\\!-\\!{\bf z}\|^{2}$
$\displaystyle:=\operatorname*{argmin}_{{\bf
z}\in\mathbb{R}^{d}}\frac{1}{2\lambda A_{e}}\|{\bf
z}\|_{2}^{2}\\!+\\!(1/2\sigma_{e})\|{\bf v}\\!-\\!{\bf z}\|^{2}$
$\displaystyle:=\frac{{\bf v}}{1+(\sigma_{e}/(\lambda A_{e}))}$ (72)

1 Norm Penalty.

$\displaystyle\phi^{*}({\bf v})$ $\displaystyle{:=\sup_{{\bf
z}\in\mathbb{R}^{d}}{\bf v}^{T}{\bf z}-\phi({\bf z})}$
$\displaystyle=\sup_{{\bf z}\in\mathbb{R}^{d}}{\bf v}^{T}{\bf z}-\|{\bf
z}\|_{1}$ $\displaystyle=\begin{cases}\infty\quad\text{if\
}\|v_{j}\|_{1}>1\text{ for some $j\in\\{1,d\\}$ }\\\ 0\quad\text{otherwise}\\\
\end{cases}$

For inserting IX-C into IV if $\|v_{j}/(\lambda A_{e})\|_{1}\leq 1$ for all
$j\in\\{1,d\\}$

$\displaystyle\mathcal{DU}^{(e)}\big{\\{}{\bf v}\big{\\}}$
$\displaystyle{:=\operatorname*{argmin}_{{\bf z}\in\mathbb{R}^{d}}\lambda
A_{e}\phi^{*}\big{(}{\bf z}/(\lambda A_{e})\big{)}\\!+\\!(1/2\sigma_{e})\|{\bf
v}\\!-\\!{\bf z}\|^{2}}$
$\displaystyle\stackrel{{\scriptstyle\eqref{nomr1_penatly_formula_star}}}{{:=}}\operatorname*{argmin}_{{\bf
z}\in\mathbb{R}^{d}}(1/2\sigma_{e})\|{\bf v}\\!-\\!{\bf z}\|^{2}$
$\displaystyle:={\bf v}$ (74)

Otherwise if $\|z_{j}\|_{1}>\lambda A_{e}$ for some $j\in\\{1,d\\}$ then
$\|(\lambda A_{e}z_{j}/\|z_{j}\|_{1})\|_{2}=\|(\lambda
A_{e}z_{j}/\|z_{j}\|_{2})\|_{2}\leq\lambda A_{e}$ , so
$\mathcal{DU}^{(e)}\big{\\{}{\bf v}\big{\\}}:=(\mathcal{T}^{(\lambda
A_{e})}(v_{1}),\cdots,\mathcal{T}^{(\lambda A_{e})}(v_{n}))^{T}$.

###  IX-D Derivation of Algorithm 1

We can rewrite the optimality condition (23) as [29]

${\bf M}^{-1}\begin{pmatrix}\partial f&\mathbf{D}^{T}\\\ -\mathbf{D}&\partial
g^{*}\end{pmatrix}\begin{pmatrix}\widehat{{\bf w}}\\\
\widehat{\mathbf{u}}\end{pmatrix}\ni\mathbf{0}\mbox{ with }{\bf
M}:=\begin{pmatrix}\mathbf{T}^{-1}&-\mathbf{D}^{T}\\\
-\mathbf{D}&{\bf\Sigma}^{-1}\end{pmatrix}.$ (76)

The condition (76) involves the invertible matrix $\mathbf{M}$ which depends
on the structure of the empirical graph via the (block-) incidence matrix
$\mathbf{D}$ (6). Condition (76) also uses two positive definite block-
diagonal matrices

$\displaystyle\big{(}{\bf\Sigma}\big{)}_{e,e}:=\sigma_{e}\mathbf{I},\mbox{ for
each edge }e\\!\in\\!\mathcal{E}\mbox{, }$
$\displaystyle\big{(}\mathbf{T}\big{)}_{i,i}:=\tau_{i}\mathbf{I}\mbox{ for
each node }i\\!\in\\!\mathcal{V},$ (77)

with $\sigma_{e}\\!:=\\!1/2$ for $e\\!\in\\!\mathcal{E}$ and
$\tau_{i}\\!:=\\!1/|\mathcal{N}^{(i)}|$ for $i\\!\in\\!\mathcal{V}$.

In principle, we can use any invertible matrix $\mathbf{M}$ in (76) such that
it is a necessary and sufficient condition for a solution of (16). The
specific choice (76) serves two main purposes which are fulfilled by its off-
diagonal and main-diagonal blocks, respectively. The off-diagonal blocks in
(76) decouples (16) into two separate optimizations of (regularized variants
of) $f({\bf w})$ and $g(\mathbf{u})$, respectively. These separate
optimization steps are coupled via applications of the block-incidence matrix
$\mathbf{D}$ (6) and its transpose $\mathbf{D}^{T}$. The main diagonal blocks
of $\mathbf{M}$ serve as a preconditioning of (16) for irregular network
structures with highly varying node degrees [54].

The optimality conditions (76) can be solved by an application of the proximal
point algorithm [55, Thm. 23.41]. This algorithm consists of the iterations

$\begin{pmatrix}\widehat{{\bf w}}_{k\\!+\\!1}\\\
\widehat{\mathbf{u}}_{k\\!+\\!1}\end{pmatrix}:=\left(\mathbf{I}\\!+\\!\mathbf{M}^{-1}\begin{pmatrix}\partial
f&\mathbf{D}^{T}\\\ -\mathbf{D}&\partial
g^{*}\end{pmatrix}\right)^{-1}\begin{pmatrix}\widehat{{\bf w}}_{k}\\\
\widehat{\mathbf{u}}_{k}\end{pmatrix}.$ (78)

The iteration (78) can be rewritten as [29]

$\displaystyle\widehat{{\bf w}}_{k+1}$
$\displaystyle\\!:=\\!\big{(}\mathbf{I}\\!+\\!\mathbf{T}\partial
f\big{)}^{-1}(\widehat{{\bf
w}}_{k}-\mathbf{T}\mathbf{D}^{T}\widehat{\mathbf{u}}_{k})$ (79)
$\displaystyle\widehat{\mathbf{u}}_{k+1}$
$\displaystyle\\!:=\\!\big{(}\mathbf{I}\\!+\\!{\bf\Sigma}\partial
g^{*}\big{)}^{-1}(\widehat{\mathbf{u}}_{k}\\!+\\!{\bf\Sigma}\mathbf{D}(2\widehat{{\bf
w}}_{k\\!+\\!1}\\!-\\!\widehat{{\bf w}}_{k})).$ (80)

For convex $f({\bf w})$ and $g(\mathbf{u})$, the sequences $\widehat{{\bf
w}}_{k}$ and $\widehat{\mathbf{u}}_{k}$ obtained from (79), (80) converge to
solutions of (16) and (19), respectively. This convergence is guaranteed for
any initialization $\widehat{{\bf w}}_{0}$ and $\widehat{\mathbf{u}}_{0}$
[54].

To compute the updates (79) and (80), we can use the identity [54]

$\displaystyle\big{(}\mathbf{I}\\!+\\!\mathbf{S}\partial h\big{)}^{-1}({\bf
a})\\!:=\\!\operatorname*{argmin}_{{\bf b}\in\mathbb{R}^{d}}h({\bf
b})\\!+\\!(1/2)\|{\bf a}\\!-\\!{\bf b}\|^{2}_{\mathbf{S}^{-1}},$ (81)

which is valid for any convex function $h$ and positive definite matrix
$\mathbf{S}\in\mathbb{R}^{d\times d}$. By inserting (6) into (79), and using
(81), we can decompose the primal update into node-wise updates

$\mathbf{w}^{(i)}_{k+1}\\!:=\\!\mathcal{PU}^{(i)}\left\\{\mathbf{w}^{(i)}_{k}-\tau_{i}\sum_{e\in\mathcal{E}}D_{e,i}\mathbf{u}^{(e)}_{k}\right\\}\mbox{,
for each node }i\in\mathcal{V}.$ (82)

Here, we used the primal update operator (32). Similarly, the dual update (80)
decomposes into independent updates for each edge $e\in\mathcal{E}$,

$\displaystyle\mathbf{u}^{(e)}_{k+1}\\!:=\\!$
$\displaystyle\mathcal{DU}^{(e)}\big{\\{}\mathbf{u}^{(e)}_{k}\\!+\\!\sigma_{e}\big{(}2\big{(}\mathbf{w}^{(e_{+})}_{k+1}\\!-\\!\mathbf{w}^{(e_{-})}_{k+1}\big{)}\\!-\\!\big{(}\mathbf{w}^{(e_{+})}_{k}\\!-\\!\mathbf{w}^{(e_{-})}_{k}\big{)}\big{)}\big{\\}},$
$\displaystyle\mbox{ for each }e\in\mathcal{E},$ (83)

with the dual update operator (IV).

###  IX-E More Numerical Experiments

####  IX-E1 SBM Low-Dimensional Linear Regression

Figure 7 and 8 depict results for data with empirical graph $\mathcal{G}^{(\rm
SBM)}$ consisting of two clusters,
$|\mathcal{C}^{(1)}|=|\mathcal{C}^{(2)}|=100$. The intra-cluster edge
probability is fixed to $p_{\rm in}\\!=\\!0.5$, but different choices for the
inter-cluster edge probability $p_{\rm out}$ are used. The local datasets
$\mathcal{X}^{(i)}$ consist of $m_{i}\\!=\\!5$ data points with feature
vectors ${\bf x}^{(i,r)}\in\mathbb{R}^{2}$ and labels generated via (37) where
$\overline{{\bf w}}^{(1)}=\big{(}2,2\big{)}^{T}$ and $\overline{{\bf
w}}^{(2)}=\big{(}-2,2\big{)}^{T}$. For each choice of $\rho$ (fraction of
accessible local datasets), $\sigma$ (label noise) and $p_{\rm out}$ (inter-
cluster edge probability), we draw $5$ realizations of $\mathcal{G}^{(\rm
SBM)}$ and compute the average and standard deviation of the resulting MSE
values (39) obtained for each realization. Figure 7 depicts the average and
standard deviation for the MSE values in the noise-less case where $\sigma=0$
in (37). Figure 8 depicts average and standard deviation of the MSE obtained
when having access to a fraction of $\rho=0.6$ local datasets.

$0$$0.1$$0.2$$0.3$$0.4$$0.5$$0.6$$10^{-16}$$10^{-12}$$10^{-8}$$10^{-4}$$10^{0}$$p_{\rm
out}$MSE$\phi(\mathbf{u})\\!=\\!\|\mathbf{u}\|_{1},\lambda\\!=\\!10^{-2}$$\rho\\!=\\!2/10$$\rho\\!=\\!4/10$$\rho\\!=\\!6/10$

$0$$0.1$$0.2$$0.3$$0.4$$0.5$$0.6$$10^{-16}$$10^{-12}$$10^{-8}$$10^{-4}$$10^{0}$$p_{\rm
out}$MSE$\phi(\mathbf{u})\\!=\\!\|\mathbf{u}\|_{2},\lambda\\!=\\!10^{-2}$
(nLasso)$\rho\\!=\\!2/10$$\rho\\!=\\!4/10$$\rho\\!=\\!6/10$

$0$$0.1$$0.2$$0.3$$0.4$$0.5$$0.6$$10^{-16}$$10^{-12}$$10^{-8}$$10^{-4}$$10^{0}$$p_{\rm
out}$MSE$\phi(\mathbf{u})\\!=\\!\|\mathbf{u}\|_{2}^{2},\lambda\\!=\\!10^{-2}$
(“MOCHA”)$\rho\\!=\\!0.2$$\rho\\!=\\!0.4$$\rho\\!=\\!0.6$

Figure 7:  MSE (39) incurred by Algorithm 1 used to learn the local model
parameters (37) of $\mathcal{G}^{(\rm SBM)}$. These results are obtained for
local datasets obtained from a noise-free linear model ($\sigma=0$ in (37))
and which are available only for the nodes in the subset $\mathcal{M}$(see
(38)) of relative size $\rho=|\mathcal{M}|/|\mathcal{V}|$.

$0$$0.1$$0.2$$0.3$$0.4$$0.5$$0.6$$10^{-4}$$10^{-3}$$10^{-2}$$10^{-1}$$10^{0}$$p_{\rm
out}$MSE$\phi(\mathbf{u})\\!=\\!\|\mathbf{u}\|_{1},\lambda\\!=\\!10^{-2}$$\sigma\\!=\\!10^{-2}$$\sigma\\!=\\!10^{-1}$$\sigma\\!=\\!1$

$0$$0.1$$0.2$$0.3$$0.4$$0.5$$0.6$$10^{-4}$$10^{-3}$$10^{-2}$$10^{-1}$$10^{0}$$p_{\rm