# A Large-Scale Simulation Method for Neuromorphic Circuits

Amir Shahhosseini, Thomas Chaffey and Rodolphe Sepulchre A. Shahhosseini is
with KU Leuven, Department of Electrical Engineering (STADIUS), KasteelPark
Arenberg, 10, B-3001 Leuven, Belgium<EMAIL_ADDRESS>T.
Chaffey is with the University of Cambridge, Department of Engineering,
Trumpington Street, CB2 1PZ<EMAIL_ADDRESS>R. Sepulchre is with both the KU
Leuven, Department of Electrical Engineering (STADIUS), KasteelPark Arenberg,
10, B-3001 Leuven, Belgium, and the University of Cambridge, Department of
Engineering, Trumpington Street, CB2 1PZ<EMAIL_ADDRESS>research leading to these results has received funding from the European
Research Council under the Advanced ERC Grant Agreement SpikyControl
n.101054323. The work of T. Chaffey was supported by Pembroke College,
Cambridge.

###### Abstract

Splitting algorithms are well-established in convex optimization and are
designed to solve large-scale problems. Using such algorithms to simulate the
behavior of nonlinear circuit networks provides scalable methods for the
simulation and design of neuromorphic systems. For circuits made of linear
capacitors and inductors with nonlinear resistive elements, we propose a
splitting that breaks the network into its LTI lossless component and its
static resistive component. This splitting has both physical and algorithmic
advantages and allows for separate calculations in the time domain and in the
frequency domain. To demonstrate the scalability of this approach, a network
made from one hundred neurons modeled by the well-known FitzHugh-Nagumo
circuit with all to all diffusive coupling is simulated.

##  I Introduction

Neuromorphic engineering promises a new technological paradigm for computing
machines. Neuromorphic technologies are designed to be more energy efficient
than their digital counterparts [1] to handle ill-conditioned inputs
significantly better [2]. In recent years, spiking neural networks [3],
neuromorphic computing devices [4], and neuromorphic sensors [5] have drawn
the attention of scientists in various disciplines. This interest creates new
demand in the simulation, analysis, and design of large-scale electrical
circuits that are at the heart of any neuromorphic system.

The common methods of simulating spiking neurons rely on the numerical
integration of nonlinear differential equations [6]. This approach enabled the
study and simulation of simple neurons and contributed to the development of
computational neuroscience. Nevertheless, numerical integration methods are
known to struggle with neural networks with spiking behaviors due to the
stiffness of the governing equations. Additionally, such methods do not scale
up adequately with an increase in the size of the network and become
prohibitive. It is also impractical to analyze the robustness of the neural
networks against the variability and uncertainty of the components. This
motivates and necessitates exploring alternative methods for the efficient
simulation and analysis of large-scale electrical circuits that represent
neuromorphic systems.

This paper takes inspiration from passive linear circuit theory in developing
an alternative. There, studying the network is computationally tractable, and
passivity provides insight into the physics of the problem [7]. Motivated to
extend these properties to networks with nonlinear resistors, Minty put
forward the idea of monotonicity that corresponds to incremental passivity
[8]. In the meantime, Rockafellar’s works on convex analysis have placed
monotonicity at the core of large-scale convex optimization algorithms [9].
Thus, like passivity in the case of linear networks, monotonicity is the
concept that reconciles the physics of nonlinear circuits with algorithmic
tractability. Thus, monotonicity is the key concept that enables the
development of new efficient methods for the simulation of large-scale
nonlinear networks.

To this end, the recent papers [10, 11] explore monotone circuits, formulate
their dynamics as a monotone zero inclusion problem, and use splitting
algorithms to obtain their solution [12]. Splitting algorithms are central to
these studies since they provide computational tractability to large-scale
optimization methods, by allowing computational steps to be performed
separately for each circuit component. The results of these studies confirm
that methods of large-scale convex optimization can be used instead of
numerical integration to efficiently simulate monotone circuits and networks.
Nevertheless, models of spiking neurons contain a mixture of monotone and
anti-monotone elements that represent a mixed feedback system [13]. The recent
work [14] explores a method that is inspired by convex-concave programming
[15] to solve mixed monotone systems using this approach, but no attention is
given to the splitting itself.

Similar to how monotonicity signifies the physics of the element, adequate
splitting should exploit the circuit topology. Thus, in addition to providing
physical intuition into the structure of the network, a proper splitting also
makes simulation and computations significantly faster. This paper exploits a
splitting of the circuit into its lossless storage elements (capacitors and
inductors) and its dissipative nonlinear resistors. With this splitting, the
zero inclusion problem is broken into two subproblems, one being a linear
time-invariant (LTI) lossless operator (representing inductors, capacitors,
and interconnections with no dissipation) and the other a mixed-monotone
operator (the difference of two monotone operators representing the nonlinear
resistors). This is similar to the element extraction approach adopted in
[11], but the reactive elements are grouped with the interconnection matrix,
rather than with the dissipative elements. Computational steps for the
lossless LTI operator are performed in the frequency domain, while
computational steps for the nonlinear resistive operator are performed in the
time domain. We call this approach _time-frequency splitting_.

The remainder of this paper is organized as follows. Section 2 formulates the
problem, and describes the energy-based splitting and its connection with the
energy-based methods of the literature. Section 3 presents an algorithm that
can solve this zero inclusion problem. It will also explain how using both the
time and frequency domains can allow the exploitation of the LTI structure of
the lossless part and make the computations efficient. Section 4 demonstrates
the presented algorithm by simulating the FitzHugh-Nagumo model [16] of a
spiking neuron and, extends this idea to a network of 100 heterogeneous
neurons with diffusive coupling.

##  II Energy-based Splitting of Nonlinear Circuits

In this section, the energy-based splitting of RLC networks with static
nonlinear resistors is introduced. Furthermore, the relation between this
representation and existing energy-based representations of nonlinear circuits
is described.

###  II-A Preliminaries

We begin by introducing some necessary preliminaries from the theory of
monotone operators. $\mathcal{L}^{2}_{\mathbb{T}}$ is the Hilbert space of
square-integrable signals over the time axis $\mathbb{T}$, equipped with the
inner product

$\left\langle{u},{y}\right\rangle=\int_{0}^{\mathbb{T}}u^{\top}(t)y(t)dt<\infty$
(1)

We use the shorthand notation $\mathcal{L}^{2}$ to denote
$\mathcal{L}^{2}_{[0,\infty)}$, and $\mathcal{L}^{2}_{T}$ to denote
$\mathcal{L}^{2}_{[0,T]}$. The latter may be associated with $T$-periodic
signals, restricted to a single period. A trajectory $v(t)$ is $T$-periodic if
$v(t+T)=v(t)$ for any $t$. The Hilbert space of $T$-periodic signals is of
interest since we mostly deal with periodic behavior that represents rhythmic
behavior. $l^{2}_{T}$ is the discrete-time counterpart of
$\mathcal{L}^{2}_{T}$.

Definition 1. An operator $A$ on the space $\mathcal{L}^{2}$ is a set-valued
mapping $A:\mathcal{L}^{2}\rightrightarrows\mathcal{L}^{2}$. The graph of the
operator $A$ is defined as

$\text{Gra}A=\\{(u,y)|y\in Au\\}\subseteq\mathcal{L}^{2}\times\mathcal{L}^{2}$
(2)

Definition 2. An operator $A:\mathcal{L}^{2}\rightrightarrows\mathcal{L}^{2}$
is monotone if

$\left\langle{u_{1}-u_{2}},{A(u_{1})-A(u_{2})}\right\rangle\geq 0$ (3)

for all $u_{1},u_{2}\in\mathrm{dom}\,{A}$. The operator $A$ is said to be
_maximal monotone_ if its graph is not properly contained in the graph of any
other monotone operator.

The subgradient of a closed, convex, and proper function is maximal monotone.
For the converse to be true, a stronger property is required: maximal cyclic
monotonicity.

Definition 3. An operator $A:l_{2}\rightrightarrows l_{2}$ is n-cyclic
monotone if, for every $(u_{1},u_{2},...,u_{n+1})\in{l_{2}}^{n+1}$ and every
$(y_{1},y_{2},...,y_{n})\in{l_{2}}^{n}$ where $y_{i}=A(u_{i})$

$\sum_{i=1}^{n}\left\langle y_{i+1}-y_{i},u_{i}\right\rangle\leq 0$ (4)

If $A$ is $n$-cyclically monotone for every integer $n\geq 2$, then $A$ is
cyclically monotone. An operator A is maximal cyclic monotone if it is maximal
monotone and cyclic monotone.

Theorem 1. (Rockafellar’s theorem [9]) An operator
$A:\mathcal{H}\rightrightarrows\mathcal{H}$ is maximal cyclic monotone if and
only if it is the subgradient of a closed, convex, and proper function from
$\mathcal{H}$ to $(-\infty,\infty]$.

Definition 4. Given an operator
$A:\mathcal{L}^{2}\rightrightarrows\mathcal{L}^{2}$, its resolvent, denoted
$J_{\alpha A}:\mathcal{L}^{2}\rightrightarrows\mathcal{L}^{2}$ is defined as

$J_{\alpha A}=(I+\alpha A)^{-1}$ (5)

where $\alpha$ is a positive real number.

###  II-B Neuromorphic Circuits and Energy-based Splitting

####  II-B1 From Neural Networks to Electrical Circuits

Neuromorphic engineering is inspired by mimicking the functionality of the
nervous system. Neurons are the building blocks of these nervous systems, and
their working mechanisms have been the focus of studies within neuroscience
for decades [17]. They function and transfer information through action
potentials (spikes) that are nothing more than a sudden discharge of ions. The
ion channels that govern the passing of these ions are best modeled as
nonlinear memristive elements. The neuron’s membrane can be modeled using a
capacitor.

This fundamental alignment of the physical language of the neurons with
electrical circuits motivates modeling them using electrical elements. This
compatibility gives rise to the class of conductance-based models of neurons
[18]. This analogy between neural networks and electrical circuits motivates
studying electrical circuits.

It is possible to utilize network analysis within circuit theory to express
the dynamics of the resulting electrical circuits. However, the representation
of the dynamics is not unique. The representation of this paper relies on the
energy dissipation properties of the elements and breaks the network down into
two subcircuits. The first subcircuit contains all the resistive elements and
the second subcircuit contains the rest (inductors, capacitors, and
interconnections that are all LTI lossless elements and have no dissipation).
This is concurrent with the resistance extraction method in circuit theory
[19]. Figure 1 illustrates the splitting of the general network $\mathscr{N}$
into a skew lossless subnetwork $\mathscr{N}_{L}$ and a resistive subnetwork
$\mathscr{N}_{R}$.

Figure 1: Decomposition of the circuit network $\mathscr{N}$ into a lossless
part and a resistive part

####  II-B2 Energy-based Splitting of Nonlinear Electrical Circuits

The main idea is to express the behavior of the electrical circuit as a zero
inclusion problem with an operator governing the dynamics. Mathematically
speaking,

$A\begin{pmatrix}i(t)\\\ v(t)\end{pmatrix}=0,$ (6)

where $A$ comes from the KVL and KCL and the physics of the elements. In Eq.
(6), $i(t)$ and $v(t)$ are trajectories in the signal space
$\mathcal{L}^{2}[0,\infty)$. Still, even for simple circuits, $A$ can become
complicated and thus, it is rational to break it into smaller operators.
Splitting the operator into smaller operators helps in two ways. First, if
done properly, the smaller operators will bear physical meaning (like the
splitting of this paper). Second, splitting the operator makes the problem
computationally easier. This is because the resolvent operator of the original
operator is usually difficult to compute while the smaller operators are
simpler to deal with [12]. This is the main motivation for using splitting
algorithms in large-scale convex optimization algorithms [20]. Consequently,
it is rational to exploit the structure of the smaller operators to make the
computations even faster. Thus, Eq. (6) can be rewritten as

$S\begin{pmatrix}i(t)\\\ v(t)\end{pmatrix}+R\begin{pmatrix}i(t)\\\
v(t)\end{pmatrix}=0,$ (7)

where $R$ represents the resistive portion of the circuit and $S$ represents
the lossless portion. We call this an “energy-based splitting” as the operator
is split into those components which dissipate energy, and those which do not.

###  II-C Mixed Monotone Description of Nonlinear Resistors

The resistive elements in conductance-based models of neurons have regions of
negative resistance, and thus cannot be modeled as monotone operators. They
can, however, be modeled as the difference of monotone operators. Thus, Eq.
(7) can be reformulated as

$S\begin{pmatrix}i(t)\\\ v(t)\end{pmatrix}+M_{1}\begin{pmatrix}i(t)\\\
v(t)\end{pmatrix}-M_{2}\begin{pmatrix}i(t)\\\ v(t)\end{pmatrix}=0,$ (8)

where $M_{1}$ and $M_{2}$ are cyclic monotone operators. Cyclic monotonicity
of $M_{1}$ and $M_{2}$ implies that they are the subgradients of convex
functionals, via Rockafellar’s theorem (Theorem 1). These convex functionals
represent the energy supplied or dissipated by the resistors. The reason for
breaking the operator $R$ further into its (cyclic) monotone components is
mainly computational. In fact, most algorithms are only guaranteed to converge
when the operators are monotone, and mixed monotone operators disrupt their
functionality and convergence [12]. Beyond computational advantages, this
further splitting also bears physical meaning. $M_{1}$ corresponds to the part
of the nonlinear resistor that dissipates energy while $M_{2}$ is the part
that supplies energy. This interplay of active and dissipative elements is
central to the spiking behavior of the neuron and its rhythmic nature.

###  II-D Connection with Energy-based Representations

The splitting of this paper relies on the energy characteristics of the
circuit components. In this section, we explore the connection of this
representation with other energy-based representations of electrical circuits.
For circuits with linear elements, dynamics can be obtained by using the
Hamiltonian method. For the case of conservative systems, it is only necessary
to obtain the Hamiltonian/Co-Hamiltonian. For linear dissipative elements,
Wells introduced a “power function” (corresponding to Rayleigh’s dissipation
function in mechanics) that could be incorporated into the Hamiltonian (or
Lagrangian) equations and obtain the dynamics of linear RLC circuits [21].

The extension from linear electrical circuits to nonlinear complicated this
energy framework. The nonlinear resistive elements could no longer be
expressed in terms of simple dissipation functions. However, this issue was
resolved when Millar proposed the idea of content and co-content which were
generalizations of the “power functions” [22, 23].

In fact, the dynamics of an RLC network with nonlinear resistors can be
written as

$\left\\{\begin{matrix}\frac{d}{dt}\nabla_{i}{H}^{*}(i,v)+\nabla_{i}(v^{T}Ni)=-\nabla_{i}\mathcal{R}(i,v)\\\
\frac{d}{dt}\nabla_{v}{H}^{*}(i,v)-\nabla_{v}(v^{T}Ni)=-\nabla_{v}\mathcal{R}(i,v),\end{matrix}\right.$
(9)

where ${H}^{*}(i,v)$ is the co-Hamiltonian and is defined as

${H}^{*}(i,v)=\mathcal{T}^{*}(i)+\mathcal{V}^{*}(v).$ (10)

Here, $\mathcal{T}^{*}(i)$ and $\mathcal{V}^{*}(v)$ are the inductor’s
magnetic co-energy and capacitor’s electric co-energy [24]. The term
$\mathcal{R}(i,v)$ is defined as

$\mathcal{R}(i,v)=\mathcal{D}(i)+\mathcal{D}^{*}(v)$ (11)

where $\mathcal{D}(i)$ and $\mathcal{D}^{*}(v)$ are the corresponding content
and co-content of the resistive elements and describe the dissipated power.
The relations of Eq. (9) are identical to the zero inclusion problem of Eq.
(7). Here, the LHS only represents lossless LTI elements and the
interconnection of the elements. Moreover, the RHS also describes the
resistive part of the network.

##  III Time-Frequency Approach in Mixed Monotone Circuits

Zero inclusion problems of the form (7) are commonly solved by reformulating
them to fixed-point iterations (FPIs). These methods involve iteratively
evaluating the resolvent of the original operator. However, computing the
resolvent of complicated operators is a difficult task. The issue is also
compounded if the operator contains both dynamic elements and static nonlinear
elements. To this end, operators are split into simpler components, and the
FPIs are reformulated to iterate the resolvents of the components. This
approach offers significant computational advantages and simplifies solving
involved zero inclusion problems. The interested reader is referred to [12]
for relevant literature and the importance of monotonicity in this framework.

###  III-A Difference of Monotone Splitting Algorithms

The monotonicity of the operator in the zero inclusion problem ensures that if
the problem has a solution, it can be found systematically using
computationally tractable algorithms [12]. However, the problem is more
complex when the operator is not monotone since it can have no solutions, one
solution, or many solutions. Figure 2 demonstrates this matter for a simple
single-variable operator. Additionally, given different initializations for
the FPI algorithm, it is possible to converge to different solutions.

Figure 2: Red crosses indicate the solution to the zero inclusion problem for
a single-variable operator; (a) operator is monotone (b) operator is not
monotone

The Douglas-Rachford splitting algorithm is a powerful splitting algorithm
that finds the zero of the sum of two monotone operators. This method was
recently extended to the class of difference of monotone operators [25] that
have the form

$S\begin{pmatrix}x\end{pmatrix}+B\begin{pmatrix}x\end{pmatrix}-C\begin{pmatrix}x\end{pmatrix}=0,$
(12)

where $S$, $B$ and $C$ are monotone operators and $x$ is a trajectory in the
signal space. As can be seen, this is identical to the format of Eq. (8) and
thus, this problem can be solved using this method. The pseudo-code for the
difference of monotone Douglas-Rachford (DMDR) algorithm is provided below.
Here, $J_{\alpha A}$ represents the resolvent operator of operator $A$
(Definition 4), and $\alpha$ is the step size. $z$ defines the auxiliary
variable that is used within the structure of the algorithm and the
superscripts denote the iteration count of the algorithm.

Algorithm 1 DMDR

1:for $j=1,2,\ldots,max-iteration$ do

2: $x^{j+1}=J_{\alpha S}\left(z^{j}\right)$

3: $z^{j+1}=z^{j}-x^{j+1}+J_{\alpha B}\left(2x^{j+1}-z^{j}+\alpha
C\left(x^{j+1}\right)\right)$

4: $j=j+1$

5:end for

The approach of this paper is not limited to the DMDR algorithm and can be
used with any splitting algorithm that can solve the zero inclusion problem of
the form of Eq. (12).

###  III-B Time-Frequency Approach

In previous works [10, 11, 14], the dynamical components were discretized in
the time domain, using a backward Euler discretization of the differentiation
operator as

$D=\frac{1}{h}\left[\begin{array}[]{ccccc}+1&0&\cdots&0&-1\\\ -1&+1&&&0\\\
\vdots&&\ddots&&\vdots\\\ 0&&-1&+1&0\\\ 0&0&\cdots&-1&+1\end{array}\right]$
(13)

where $h$ is the time increment of the trajectory and the -1 at the end of the
first row enforces periodicity of the operator.

LTI systems (that are all contained in the operator $S$) are then represented
as rational functions of this matrix. Computing the resolvent of such matrices
in the time domain, which involves inversion (as indicated in Definition 4),
creates a dense matrix by inverting a structured and very sparse matrix. This
does not exploit the structure of the LTI elements and makes the computations
inefficient.

One of the innovations of the method we present here is that such resolvents
are now computed in the frequency domain. LTI systems are diagonal in the
frequency domain, and thus, the corresponding computations can be reduced from
matrix multiplication to Hadamard products. To carry out this step in the
frequency domain, let us first write the first step (line 2 of pseudo-code) as

$\displaystyle(\mathbb{I}+\alpha S)\begin{pmatrix}i^{j+1}(t)\\\
v^{j+1}(t)\end{pmatrix}=\qquad\qquad\text{ \>}$ (14)
$\displaystyle(\mathbb{I}+\alpha\begin{bmatrix}P(D)&N^{*}\\\
-N&Q(D)\end{bmatrix})\begin{pmatrix}i^{j+1}(t)\\\
v^{j+1}(t)\end{pmatrix}=\begin{pmatrix}z_{1}^{j+1}(t)\\\
z_{2}^{j+1}(t)\end{pmatrix}$

where the signal $x$ in the first step of the DMDR algorithm is broken into
$v(t)$ and $i(t)$ to provide a form that is compatible with Eq. (8). The
lossless operator $S$ is also expanded to show its structure. $P$ and $Q$ are
related to the inductors and capacitors and are functions of the derivative
operator. $N$ incorporates the interconnection and is a matrix with elements
that can only be +1, -1, or 0.

Upon taking the FFT of the signals and transforming the LTI operators to the
frequency domain, the operators $P$ and $Q$ become diagonal. Mathematically

$\underset{\text{Frequency
domain}}{\underbrace{(\mathbb{I}+\alpha\begin{bmatrix}\hat{P}(j\omega)&N^{*}\\\
-N&\hat{Q}(j\omega)\end{bmatrix})}}\begin{pmatrix}I^{j+1}(j\omega)\\\
V^{j+1}(j\omega)\end{pmatrix}=\begin{pmatrix}Z_{1}^{j+1}(j\omega)\\\
Z_{2}^{j+1}(j\omega)\end{pmatrix}$ (15)

where $\hat{P}$ and $\hat{Q}$ are now diagonal matrices and thus, new
iterations of the signals $I(j\omega)$ and $V(j\omega)$ are obtained in the
frequency domain. By taking the iFFT of these signals, the new iteration of
the signal $x$, which includes both $i(t)$ and $v(t)$, is obtained in the time
domain. For the sake of clarity, the entire process of first step is
summarized as

$\displaystyle
iFFT\begin{pmatrix}(\mathbb{I}+\alpha\begin{bmatrix}\hat{P}(j\omega)&N^{*}\\\
-N&\hat{Q}(j\omega)\end{bmatrix})^{-1}FFT\begin{pmatrix}z_{1}^{j+1}(t)\\\
z_{2}^{j+1}(t)\end{pmatrix}\end{pmatrix},$ (16)

but here, the inversion does not deal with a dense matrix but a block diagonal
one.

The computations related to the second step of the DMDR algorithm (line 3 of
the pseudo-code) should be carried out in the time domain since static
nonlinear elements are in fact diagonal in the time domain. To compute the
resolvent, efficient proximal algorithms such as the guarded Newton method can
be used [20].

By switching between the time and frequency domain, the LTI structure of the
lossless components is exploited and the computational cost of simulating the
network is drastically decreased. In previous works, and upon computing the
resolvent offline, the first step involved a matrix-vector multiplication that
has a complexity of $\mathcal{O}(n^{2})$ but this is now reduced to a vector
by vector multiplication, and the computations of this step are dominated by
the computational cost of FFT and iFFT that is $\mathcal{O}(nlog(n))$.

##  IV Results: From a Single Neuron to a Heterogeneous Network

To show the capability of the proposed method, two examples are provided. The
first example is the FitzHugh-Nagumo model of the spiking neuron. This example
is solved in detail to provide a tutorial of the approach. The second example
is a network of 100 neurons with all to all connections. This example serves
as a proof of concept for the scalability of this approach.

###  IV-A FitzHugh-Nagumo Model

The FitzHugh-Nagumo model captures the rhythmic behavior of a spiking neuron.
The circuit representation of the model is shown in Fig. (3).

Figure 3: The FitzHugh-Nagumo model circuit that is made three branches; a
linear capacitor, a linear inductor in series with a linear resistor and
finally, a static nonlinear resistor

Using KVL and KCL, the governing equations of motion can be obtained as

$\displaystyle C\frac{d}{dt}{v}=I_{g}(v)-i$ (17) $\displaystyle
L\frac{d}{dt}{i}=-Ri+v$

where $I_{g}(v)=v-\frac{v^{3}}{3}$. This can be rewritten as

$\begin{bmatrix}CD&+\mathbb{I}\\\
-\mathbb{I}&LD\end{bmatrix}\begin{bmatrix}v(t)\\\
i(t)\end{bmatrix}-\begin{bmatrix}I_{g}(v(t)))\\\ -Ri(t)\end{bmatrix}=0$ (18)

where $D$, again, defines the differentiation operator ($s$ in the Laplace
domain and $j\omega$ in the frequency domain) and $\mathbb{I}$ is the identity
operator. In the splitting of Eq. (18), the first operator is a lossless
operator that incorporates the capacitor, inductor, and interconnection of
elements. The second operator represents the resistive elements of the
circuit. This representation is exactly in the form of Eq. (7). However, the
representation can be further broken down as

$\begin{bmatrix}CD&+\mathbb{I}\\\
-\mathbb{I}&LD\end{bmatrix}\begin{bmatrix}v(t)\\\
i(t)\end{bmatrix}+\begin{bmatrix}\frac{(\cdot)^{3}}{3}&0\\\
0&R\end{bmatrix}\begin{bmatrix}v(t)\\\
i(t)\end{bmatrix}-\begin{bmatrix}\mathbb{I}&0\\\
0&0\end{bmatrix}\begin{bmatrix}v(t)\\\ i(t)\end{bmatrix}=0$ (19)

to be in the format of Eq. (8). This breaks the resistive part into the
difference of two monotone operators. With this representation, it is
straightforward to apply DMDR and obtain the solution for this circuit. Figure
4 illustrates the steady-state response of the FitzHugh-Nagumo circuit in the
space of discrete $T$-periodic square-integrable signals (${l}^{2}_{T}$) with
parameters $L=20$, $C=R=1$. To proceed with the computations, the signals
$v(t)$ and $i(t)$ were both discretized to vectors with 556 elements. The step
size ($\alpha$) for the DMDR algorithm was set to 0.1 and the sampling
frequency of the FFT was chosen to be 10 Hz. The simulation takes around 28
milliseconds to converge. The numerical integration method (two-step Adams-
Bashforth method with a step size of 0.01 seconds) computes the steady-state
in 5 milliseconds.

Figure 4: Simulation of the FitzHugh-Nagumo model; Time evolution that
corresponds to its steady-state behavior using the splitting method compared
with numerical integration

###  IV-B A Spiking Neural Network

The next example studies a network of 100 neurons modeled via the FitzHugh-
Nagumo circuit. The network has all to all connectivity meaning that each
neuron is connected to every other neuron. These synaptic connections are
modeled using diffusive coupling which means that the neurons are connected to
each other via a linear resistor. Figure 5 illustrates a small portion of this
network.

Figure 5: A small cut of the network of neurons; each node represents a neuron
that is modeled using the FitzHugh-Nagumo model and the synaptic connections
are diffusive.

Using the resistance extraction method, the dynamics of the network can be
obtained and by reformulating the dynamics of the network to the format of Eq.
(8), we have

$S=\left[\begin{array}[]{ccc:ccc}C_{1}D&&\emptyset&&&\\\
&\ddots&&&+\mathbb{I}&\\\ \emptyset&&C_{100}D&&&\\\
\hdashline&&&L_{1}D&&\emptyset\\\ &-\mathbb{I}&&&\ddots&\\\
&&&\emptyset&&L_{100}D\end{array}\right]$ (20)

and

$M_{1}=\left[\begin{array}[]{ccc:ccc}\ddots&&\emptyset&&&\\\
&\frac{(\cdot)^{3}}{3}+\sum_{j=1,j\neq
k}^{100}\frac{1}{R_{c_{k,j}}}&&&\emptyset&\\\ \emptyset&&\ddots&&&\\\
\hdashline&&&R_{1}&&\emptyset\\\ &\emptyset&&&\ddots&\\\
&&&\emptyset&&R_{100}\end{array}\right]$ (21)

where here $k$ indicates the row number and $j$ indicates the column number.
Also

$M_{2}=\left[\begin{array}[]{cccc:ccc}\mathbb{I}&\frac{1}{R_{c_{1,2}}}\mathbb{I}&\ldots&\frac{1}{R_{c_{1,100}}}\mathbb{I}&&&\\\
\frac{1}{R_{c_{2,1}}}\mathbb{I}&\mathbb{I}&&\vdots&&&\\\
\vdots&&\ddots&&&\emptyset&\\\
\frac{1}{R_{c_{100,1}}}\mathbb{I}&\ldots&&\mathbb{I}&&&\\\ \hdashline&&&&&&\\\
&\emptyset&&&&\emptyset&\\\ &&&&&&\end{array}\right]$ (22)

where $\emptyset$ is the zero operator of adequate dimension, and $L_{j}$,
$C_{j}$ and , $R_{j}$ are the inductance, capacitance and resistance of the
$j^{\text{th}}$ neuron. $R_{c_{i,j}}$ represents the strength of the diffusive
coupling between neurons $i$ and $j$ and thus $R_{c_{i,j}}=R_{c_{j,i}}$. The
nominal values of these parameters are represented in Table 1. Nevertheless,
to incorporate heterogeneity into the simulation, a random deviation, up to 20
% of the nominal values, is added to all the parameters.

TABLE I: Nominal Values of neural network parameters | $C$ | $L$ | $R$ | $R_{c}$
---|---|---|---|---
Nominal Value | 1 | 20 | 1 | 5

Simulating this network, the algorithm converges to the response in around 2.8
seconds. The numerical integration method computes the steady-state in 2.45
seconds. Compared with the example of a single neuron, the better scalability
of the proposed method is observed. All the simulations are performed in an HP
tower Z2 computer. Figure 6 demonstrates the response for a select few neurons
and synchrony of the network’s response can be observed.

Figure 6: Simulation of the spiking neural network; Steady-state for 4
randomly selected neurons demonstrates synchrony.

##  V Conclusions

The simulation of neuromorphic circuits using methods of large-scale convex
optimization is explored. To this end, an energy-based splitting is proposed
that breaks the governing operator of the nonlinear network into a linear
lossless part and a nonlinear resistive part. The nonlinear resistive part is
further broken into a difference of two cyclic monotone operators. Fixed point
iteration methods that deal with difference of monotone systems are introduced
and used to solve these networks. It is also shown that by switching between
the time domain and the frequency domain, the LTI structure of the lossless
operator can be exploited and the computations can become efficient. A large-
scale network of spiking neurons modeled by the Fitzhugh-Nagumo circuit and
diffusive coupling is simulated using this method to demonstrate its
computational tractability and scalability.

##  VI Code Availability

All the code corresponding to the methods and simulations of this paper is
publicly available at this link. The codes are written in two versions. The
first version follows the exact structure of this paper and is intended for
pedagogical purposes. The second version is optimized to provide the best
computational performance.

## References

  * [1] C. Mead, “Neuromorphic electronic systems,” _Proceedings of the IEEE_ , vol. 78, no. 10, pp. 1629–1636, 1990.
  * [2] S. Furber, “Large-scale neuromorphic computing systems,” _Journal of neural engineering_ , vol. 13, no. 5, p. 051001, 2016.
  * [3] J. L. Lobo, J. Del Ser, A. Bifet, and N. Kasabov, “Spiking neural networks and online learning: An overview and perspectives,” _Neural Networks_ , vol. 121, pp. 88–100, 2020.
  * [4] A. R. Young, M. E. Dean, J. S. Plank, and G. S. Rose, “A review of spiking neuromorphic hardware communication systems,” _IEEE Access_ , vol. 7, pp. 135 606–135 620, 2019.
  * [5] A. Vanarse, A. Osseiran, and A. Rassau, “A review of current neuromorphic approaches for vision, auditory, and olfactory sensors,” _Frontiers in neuroscience_ , vol. 10, p. 171438, 2016.
  * [6] T. J. Sejnowski, C. Koch, and P. S. Churchland, “Computational neuroscience,” _Science_ , vol. 241, no. 4871, pp. 1299–1306, 1988.
  * [7] R. E. Kalman, “Lyapunov functions for the problem of lur’e in automatic control,” _Proceedings of the national academy of sciences_ , vol. 49, no. 2, pp. 201–205, 1963.
  * [8] G. J. Minty, “Monotone networks,” _Proceedings of the Royal Society of London. Series A. Mathematical and Physical Sciences_ , vol. 257, no. 1289, pp. 194–212, 1960.
  * [9] H. H. Bauschke, P. L. Combettes, H. H. Bauschke, and P. L. Combettes, _Convex analysis and monotone operator theory in Hilbert spaces_. Springer, 2017.
  * [10] T. Chaffey and R. Sepulchre, “Monotone one-port circuits,” _IEEE Transactions on Automatic Control_ , 2023.
  * [11] T. Chaffey, S. Banert, P. Giselsson, and R. Pates, “Circuit analysis using monotone+ skew splitting,” _European Journal of Control_ , vol. 74, p. 100854, 2023.
  * [12] E. K. Ryu and W. Yin, _Large-scale convex optimization: algorithms & analyses via monotone operators_. Cambridge University Press, 2022.
  * [13] R. Sepulchre, “Spiking control systems,” _Proceedings of the IEEE_ , vol. 110, no. 5, pp. 577–589, 2022.
  * [14] A. Das, T. Chaffey, and R. Sepulchre, “Oscillations in mixed-feedback systems,” _Systems & Control Letters_, vol. 166, p. 105289, 2022.
  * [15] X. Shen, S. Diamond, Y. Gu, and S. Boyd, “Disciplined convex-concave programming,” in _2016 IEEE 55th conference on decision and control (CDC)_. IEEE, 2016, pp. 1009–1014.
  * [16] E. M. Izhikevich and R. FitzHugh, “Fitzhugh-nagumo model,” _Scholarpedia_ , vol. 1, no. 9, p. 1349, 2006.
  * [17] L. Squire, D. Berg, F. E. Bloom, S. Du Lac, A. Ghosh, and N. C. Spitzer, _Fundamental neuroscience_. Academic press, 2012.
  * [18] E. M. Izhikevich, _Dynamical systems in neuroscience_. MIT press, 2007.
  * [19] C. A. Desoer and E. S. Kuh, “Basic circuit theory,” _(No Title)_ , 1969.
  * [20] N. Parikh, S. Boyd _et al._ , “Proximal algorithms,” _Foundations and trends® in Optimization_ , vol. 1, no. 3, pp. 127–239, 2014.
  * [21] D. Wells, “A“power function”for the determination of lagrangian generalized forces,” _Journal of Applied Physics_ , vol. 16, no. 9, pp. 535–538, 1945.
  * [22] W. Millar, “Cxvi. some general theorems for non-linear systems possessing resistance,” _The London, Edinburgh, and Dublin Philosophical Magazine and Journal of Science_ , vol. 42, no. 333, pp. 1150–1160, 1951.
  * [23] D. Jeltsema, “Modeling and control of nonlinear networks: a power-based perspective,” 2005.
  * [24] D. Jeltsema and J. M. Scherpen, “Multidomain modeling of nonlinear networks and systems,” _IEEE Control Systems Magazine_ , vol. 29, no. 4, pp. 28–59, 2009.
  * [25] C.-S. Chuang, H. He, and Z. Zhang, “A unified douglas–rachford algorithm for generalized dc programming,” _Journal of Global Optimization_ , pp. 1–19, 2022.


# Federated Learning with Fair Averaging

Zheng Wang1 Xiaoliang Fan1,***Corresponding Author Jianzhong Qi2
Chenglu Wen1 Cheng Wang 1&Rongshan Yu 1 1Fujian Key Laboratory of Sensing and
Computing for Smart Cities, School of Informatics, Xiamen University, Xiamen,
China
2School of Computing and Information Systems, University of Melbourne,
Melbourne, Australia<EMAIL_ADDRESS>
<EMAIL_ADDRESS><EMAIL_ADDRESS>{clwen, cwang,
<EMAIL_ADDRESS>
###### Abstract

Fairness has emerged as a critical problem in federated learning (FL). In this
work, we identify a cause of unfairness in FL – _conflicting_ gradients with
large differences in the magnitudes. To address this issue, we propose the
federated fair averaging (FedFV) algorithm to mitigate potential conflicts
among clients before averaging their gradients. We first use the cosine
similarity to detect gradient conflicts, and then iteratively eliminate such
conflicts by modifying both the direction and the magnitude of the gradients.
We further show the theoretical foundation of FedFV to mitigate the issue
conflicting gradients and converge to Pareto stationary solutions. Extensive
experiments on a suite of federated datasets confirm that FedFV compares
favorably against state-of-the-art methods in terms of fairness, accuracy and
efficiency. The source code is available at https://github.com/WwZzz/easyFL.

##  1 Introduction

Federated learning (FL) has emerged as a new machine learning paradigm that
aims to utilize clients’ data to collaboratively train a global model while
preserving data privacy McMahan et al. (2017). The global model is expected to
perform better than any locally trained model, since it has much more data
available for the training. However, it is difficult for the global model to
treat each client fairly Kairouz et al. (2019); Li et al. (2019); Mohri et al.
(2019). For example, a global model for face recognition may work well for
younger users (clients) but may suffer when being used by more senior users,
as the younger generation may use mobile devices more frequently and
contribute more training data.

Techniques have been proposed to address the fairness issue in FL. AFL Mohri
et al. (2019) aims at good-intent fairness, which is to protect the worst-case
performance on any client. This technique only works on a small network of
dozens of clients because it directly treats each client as a domain, which
may suffer in generalizability. Li _et al._ Li et al. (2019) seek to balance
the overall performance and fairness using a fair resource allocation method.
Hu _et al._ Hu et al. (2020) find a common degrade direction for all clients
without sacrificing anyone’s benefit and demonstrate robustness to inflating
loss attacks. Li _et al._ Li et al. (2020a) further explore a trade-off
between a more general robustness and fairness, and they personalize each
client’s model differently.

Figure 1: When there are conflicting gradients ($g_{a}\cdot g_{b}\textless
0,g_{a}\cdot g_{c}\textless 0,g_{b}\cdot g_{c}\textless 0$) and large
differences in the gradient magnitudes,
$||g_{a}||\textgreater||g_{b}||\textgreater||g_{c}||$, the average of the
original gradients $\bar{\theta}$ will be dominated by $g_{a}$, which will be
far away from the global optimal $\theta^{*}$, resulting in unfairness to the
dominated clients $b$ and $c$. We argue that if projecting the gradients to
mitigate the conflicts before averaging, the update direction
$\bar{g}^{\prime}$ will be corrected, which is closer to the optimal and
fairer.

Different from these studies, we observe that _conflicting_ gradients with
large differences in the magnitudes might fatally bring unfairness in FL.
Since gradients are calculated locally by the clients, an update direction of
some clients may hinder the model performance on other clients. The gradient
conflicts seem not detrimental on their own since it has been shown that a
simple average of gradients is valid to optimize the global objective function
Li et al. (2020c). However, when the conflicting gradients have large
differences in their magnitudes, model accuracy can suffer significant
reductions for some of the clients, as shown in Figure 1.

We highlight three characteristics of FL that make it critical in ensuring
fairness for the clients.

  * •

Heterogeneous and imbalanced data. Different clients have their own data
characteristics and hence potentially different data distributions. When the
data is non-IID distributed, the divergence among the gradients of difference
clients gets even larger Zhao et al. (2018). These contribute conflicting
gradients in FL. Further, different clients may also have different dataset
sizes. If the same local batch size is used across the clients, a client with
more data will take more steps to train the local model, thus leading to large
differences in gradient magnitudes. As a result, the selected users in a
communication round can largely conflict with each other, leading to an unfair
update to some of them – we call such conflicts _internal conflicts_.

  * •

Party selection. In each communication round, the server only randomly selects
a subset of clients to train the global model with their local datasets. There
is no guarantee that the data distribution in each round is representative of
the real population distribution. The data distribution can be highly
imbalanced when the server repeatedly chooses a particular type of clients,
e.g., due to sampling strategies based on computing power, latency, battery
etc., the average of gradients may become dominated by the repeatedly selected
clients. Therefore, apart from the internal conflicts, we also identify
_external conflicts_ between those selected and those ignored by the server.

  * •

Client dropping out. Another source of _external conflicts_ comes from network
interruption. When a client drops out due to network interruption (or other
hardware issues), its magnitude is considered as 0. An update in the global
model may also be unfair for the client.

To address the issues above, we propose Federated Fair aVeraging (FedFV)
algorithm to mitigate the conflicts among clients before averaging their
gradients. We first use the cosine similarity to detect gradient conflicts and
then iteratively eliminate such conflicts by modifying both the direction and
magnitude of the gradients. We further show how well FedFV can mitigate the
conflicts and how it converges to either a pareto stationary solution or the
optimal on convex problems. Extensive experiments on a suite of federated
datasets confirm that FedFV compares favorably against state-of-the-art
methods in terms of fairness, accuracy and efficiency.

The contributions of this work are summarized as follows:

  * •

We identify two-fold gradients conflicts (i.e., internal conflict and external
conflict) that are major causes of unfairness in FL.

  * •

We propose FedFV that complements existing fair FL systems and prove its
ability to mitigate two gradient conflicts and converge to Pareto stationary
solutions.

  * •

We perform extensive experiments on a suite of federated datasets to validate
the competitiveness of FedFV against state-of-the-art methods in terms of
fairness, accuracy and efficiency.

##  2 Preliminaries

FL aims to find a shared parameter $\theta$ that minimizes the weighted
average loss of all clients McMahan et al. (2017):

$\displaystyle\min_{\theta}F(\theta)=$
$\displaystyle\sum_{k=1}^{K}p_{k}F_{k}(\theta)$ (1)

where $F_{k}(\theta)$ denotes the local objective of the $k$th client with
weight $p_{k}$, $p_{k}\geq 0$ and $\sum_{k=1}^{K}p_{k}=1$. The local objective
is usually defined by an empirical risk on a local dataset, i.e.,
$F_{k}(\theta)=\frac{1}{n_{k}}\sum_{\xi_{i}\in D_{k}}l(\theta,\xi_{i})$, where
$D_{k}$ denotes the local dataset of the $k$th client, and $n_{k}$ is the size
of $D_{k}$. To optimize this objective, McMahan _et al._ McMahan et al. (2017)
propose _FedAvg_ , an iterative algorithm where the server randomly samples a
subset $S_{t}$ of $m$ clients, $0\textless m\leq K$, to train the global model
with their own datasets and aggregates local updates
$G_{t}=\\{g_{1}^{t},g_{2}^{t},...,g_{m}^{t}\\}$ into an average
$g^{t}=\sum_{k=1}^{m}p_{k}g_{k}^{t}$ with weight
$p_{k}=\frac{n_{k}}{\sum_{j=1}^{m}n_{j}}$ in the $t$th iteration. FedAvg has
been shown to be effective in minimizing the objective with low communication
costs while preserving privacy. However, it may lead to an unfair result where
accuracy distribution is imbalanced among different clients Li et al. (2019);
Mohri et al. (2019).

In this paper, we aim to minimize Equation (1) while achieving a fair accuracy
distribution among the clients.

##  3 Proposed Approach

We first present our approach, FedFV, for easing the negative impact of
conflicting gradients with largely different magnitudes on FL fairness. We
then give a theoretical analysis on how well FedFV achieves fairness and its
convergence. We start by defining gradient conflicts following Yu et al. Yu et
al. (2020).

######  Definition 1.

Client $i$’s gradient $g_{i}\in R^{d}$ conflicts with client $j$’s gradient
$g_{j}\in R^{d}$ iff $g_{i}\cdot g_{j}\textless 0$.

###  3.1 FedFV

FedFV aims to mitigate the conflicts among clients before averaging their
gradients. In the $t$th communication round of FedFV, after receiving the
selected clients’ updates $G_{t}=\\{g_{1}^{t},g_{2}^{t},...,g_{m}^{t}\\}$ and
training losses $L_{t}=\\{l_{1}^{t},l_{2}^{t},...,l_{m}^{t}\\}$, the server
updates the gradient history $GH$ to keep a trace of the latest gradient of
each client. Then, the server sorts gradients in $G_{t}$ in ascending order of
their corresponding clients’ losses to obtain $PO_{t}$, which provides the
order of each gradient to be used as a projection target. Finally, FedFV
mitigates the internal conflicts and the external conflicts sequentially.
Algorithm 1 summarizes the steps of FedFV.

Algorithm 1 FedFV

Input:$T,m,\alpha,\tau,\eta,\theta^{0},p_{k},k=1,...,K,$

1: Initialize $\theta_{0}$ and gradient history $GH=[]$.

2: for $t=0,1,...,T-1$ do

3: Server samples a subset $S_{t}$ of $m$ clients with the prob. $p_{k}$ and
sends the model $\theta^{t}$ to them.

4: Server receives the updates $g_{k}^{t}$ and the training loss $l_{k}^{t}$
from each client $k\in S_{t}$, where
$g_{k}^{t}\leftarrow\theta^{t}-\theta_{k}^{t}$ and $\theta_{k}^{t}$ is updated
by client $k$. Then server updates the trace
$GH=[g_{1}^{t_{1}},g_{2}^{t_{2}},...,g_{K}^{t_{K}}]$ of the latest updates of
all clients.

5: Server sorts the clients’ updates into a projecting order list
$PO_{t}=[g_{i_{1}}^{t},...,g_{i_{m}}^{t}]$ based on their losses, where
$l_{i_{j}}^{t}\leq l_{i_{j+1}}^{t}$.

6: $g^{t}\leftarrow$MitigateInternalConflict$(PO_{t},G_{t},\alpha)$

7: if $t\geq\tau$ then

8: $g^{t}\leftarrow$MitigateExternalConflict$(g^{t},GH,\tau)$

9: end if

10: $g^{t}=g^{t}/||g^{t}||*||\frac{1}{m}\sum_{k}^{m}g_{k}^{t}||$

11: Server updates the model $\theta^{t+1}\leftarrow\theta^{t}-g^{t}$

12: end for

###  3.2 Mitigating Internal Conflicts

We first handle the fairness among the selected clients in each communication
round.

######  Definition 2.

In the $t$th communication round, there are internal conflicts if there is at
least a pair of client gradients $<g_{i}^{t},g_{j}^{t}>$ such that $g_{i}^{t}$
conflicts with $g_{j}^{t}$, where $g_{i}^{t},g_{j}^{t}\in
G_{t}=\\{g_{1}^{t},g_{2}^{t},...,g_{m}^{t}\\}$ and $G^{t}$ is the selected
clients’ updates.

Internal conflicts account for the unfairness among selected clients. Consider
learning a binary classifier. if clients with data of one class are in the
majority in a communication round and there exist conflicts between gradients
of these two classes, the global model will be updated to favor clients of the
majority class, sacrificing the accuracy on clients of the other class.
Further, even when the proportion of clients of the two classes is balanced,
the magnitudes of gradients may still largely differ due to different dataset
sizes on different clients.

To address the internal conflicts, FedFV iteratively projects a client’s
gradient onto the normal plane of another client with a conflicting gradient.
Instead of projecting in a random order, we design a loss-based order to
decide which gradient should be the projecting target of other gradients based
on Theorem 1. We show that the later a client gradient is used as the
projection target, the fewer conflicts it will have with the final average
gradient computed by FedFV. Algorithm 2 summarizes our steps to mitigate
internal conflicts.

######  Theorem 1.

Suppose there is a set of gradients $G=\\{g_{1},g_{2},...,g_{m}\\}$ where
$g_{i}$ always conflicts with $g_{j}^{(t_{j})}$ before projecting
$g_{j}^{(t_{j})}$ to $g_{i}$’s normal plane, and $g_{i}^{(t_{i})}$ is obtained
by projecting $g_{i}$ to the normal planes of other gradients in $G$ for
$t_{i}$ times. Assuming that
$|\cos{<g_{i}^{(t_{i})},g_{j}^{(t_{j})}>}|\leq\epsilon,0\textless\epsilon\leq
1$, for each $g_{i}\in G$, as long as we iteratively project $g_{i}$ onto
$g_{k}$’s normal plane (skipping $g_{i}$ it self) in the ascending order of
$k$ where $k=1,2,...,m$, the larger $k$ is, the smaller the upper bound of
conflicts between the average
$\bar{g}^{\prime}=\frac{1}{m}\sum_{i=1}^{m}g_{i}^{(m)}$ and $g_{k}$ is.

###### Proof.

See Appendix A.1. ∎

Since projecting gradients favors for the gradient that serves as the
projecting target later, we put the clients with larger training losses at the
end of the projecting order list $PO_{t}$ in the $t$th communication round to
improve the model performance on the less well trained data. We allow $\alpha
m$ clients to keep their original gradients to further enhance fairness. We
will detail $\alpha$ in Section 3.4.

Algorithm 2 MitigateInternalConflict

Input: $PO_{t},G_{t},\alpha$

1: Server selects a set $S_{t}^{\alpha}$ of the top $\alpha$ clients in
$PO_{t}$.

2: Server sets $g_{k}^{PC}\leftarrow g_{k}^{t}$.

3: for each client $k\in S_{t}^{\alpha}$ in parallel do

4: for each $g_{i_{j}}^{t}\in PO_{t},j=1,...,m$ do

5: if $g_{k}^{PC}\cdot g_{i_{j}}^{t}\textless 0$ and $k\neq i_{j}$ then

6: Server computes $g_{k}^{PC}\leftarrow g_{k}^{PC}-\frac{(g_{i_{j}}^{t})\cdot
g_{k}^{PC}}{||g_{i_{j}}^{t}||^{2}}g_{i_{j}}^{t}$.

7: end if

8: end for

9: end for

10: $g^{t}\leftarrow\frac{1}{m}\sum_{k=1}^{m}g_{k}^{PC}$

11: return $g^{t}$

###  3.3 Mitigate External Conflicts

Due to party selection and client dropping out, there is a sample bias during
each communication round in FL Kairouz et al. (2019). A client that is not
selected in the $t$th round can suffer a risk of being forgotten by the model
whenever the combined update $g^{t}$ conflicts with its imaginary gradient
$g_{img}^{t}$. We can consider such clients to have a weight of zero in each
round. However, we cannot directly detect conflicts for such clients because
we have no access to their real gradients. To address this issue, we estimate
their real gradients according to their recent gradients, and we call such
gradient conflicts external conflicts:

######  Definition 3.

In the $t$th communication round, there are external conflicts if there is a
client $h\notin S_{t}$ whose latest gradient $g_{h}^{t-k}$ conflicts with the
combined update $g^{t}$, where $0\textless
k\textless\tau,0\textless\tau\textless t$.

An external conflict denotes the conflict between the assumed gradient of a
client that has not been selected and the combined update $g^{t}$. We add
extra steps to prevent the model from forgetting data of clients beyond
current selection. We also iteratively project the update $g^{t}$ onto the
normal plane of the average of conflicting gradients for each previous round
in a time-based order. The closer the round is, the later we make it the
target to mitigate more conflicts between the update $g^{t}$ and more recent
conflicting gradients according to Theorem 1. We scale the length of the
update to $||\frac{1}{m}\sum_{k}^{m}g_{k}^{t}||$ in the end because the length
of all gradients can be enlarged by the projection. Algorithm 3 summarizes our
steps to mitigate external conflicts.

Algorithm 3 MitigateExternalConflict

Input: $g^{t},GH,\tau$

1: for round $t-i,i=\tau,\tau-1,...,1$ do

2: $g_{con}\leftarrow 0$

3: for each client $k=1,2,...,K$ do

4: if $t_{k}=t-i$ then

5: if $g^{t}\cdot g_{k}^{t_{k}}\textless 0$ then

6: $g_{con}\leftarrow g_{con}+g_{k}^{t_{k}}$

7: end if

8: end if

9: end for

10: if $g^{t}\cdot g_{con}\textless 0$ then

11: Server computes $g^{t}\leftarrow g^{t}-\frac{g^{t}\cdot
g_{con}}{||g_{con}||^{2}}g_{con}$.

12: end if

13: end for

14: return $g^{t}$

###  3.4 Analysis

We analyze the capability of FedFV to mitigate gradient conflicts and its
convergence. We show how FedFV can find a pareto stationary solution on convex
problems.

######  Theorem 2.

Suppose that there is a set of gradients $G=\\{g_{1},g_{2},...,g_{m}\\}$ where
$g_{i}$ always conflicts with $g_{j}^{(t_{j})}$ before projecting
$g_{j}^{(t_{j})}$ to $g_{i}$’s normal plane and $g_{i}^{(t_{i})}$ is obtained
by projecting $g_{i}$ to the normal planes of different gradients in $G$ for
$t_{i}$ times. If
$\epsilon_{1}\leq|\cos{<g_{i}^{(t_{i})},g_{j}^{(t_{j})}>}|\leq\epsilon_{2},0\textless\epsilon_{1}\leq\epsilon_{2}\leq
1$, then as long as we iteratively project $g_{i}$ onto $g_{k}$’s normal plane
(skipping $g_{i}$ it self) in the ascending order of $k$ where $k=1,2,...,m$,
the maximum value of $|g_{k}\cdot\bar{g}^{\prime}|$ is bounded by
$\frac{m-1}{m}(\max_{i}{||g_{i}||})^{2}f(m,k,\epsilon_{1},\epsilon_{2})$,
where
$f(m,k,\epsilon_{1},\epsilon_{2})=\frac{\epsilon_{2}^{2}(1-\epsilon_{1}^{2})^{\frac{1}{2}}(1-(1-\epsilon_{1}^{2})^{\frac{m-k}{2}})}{1-(1-\epsilon_{1}^{2})^{\frac{1}{2}}}$.

###### Proof.

See Appendix A.2. ∎

According to Theorem 2, it is possible for FedFV to bound the maximum conflict
for any gradient by choosing $k$ to let
$f(m,k,\epsilon_{1},\epsilon_{2})\textless\epsilon_{2}$ with any given
$m,m\geq 2$, since the possible maximum value of conflicts between $g_{k}$ and
the original average gradient is
$\max_{k}{|g_{k}\cdot\bar{g}|}\leq\frac{m-1}{m}\epsilon_{2}(\max_{k}{||g_{k}||})^{2}$.

In practice, we mitigate all the conflicts that have been detected to limit
the upper bound of gradient conflicts of clients. We allow $\alpha m$ clients
with large training losses to keep their original gradients to further enhance
fairness. When $\alpha=1$, all clients keep their original gradients so that
FedFV covers FedAvg, and when $\alpha=0$, all clients are enforced to mitigate
conflicts with others. Parameter $\alpha$ thus controls the degree of
mitigating conflicts and can be applied to seek for a balance.

Next, we show that FedFV can find a pareto stationary point according Theorems
3 and 4.

######  Theorem 3.

Assume that there are only two types of users whose objective functions are
$F_{1}(\theta)$ and $F_{2}(\theta)$, and each objective function is
differentiable, $L$-smooth and convex. For the average objective
$F(\theta)=\frac{1}{2}\sum F_{i}(\theta)$, FedFV with stepsize
$\eta\leq\frac{1}{L}$ will converge to either 1) a pareto stationary point, 2)
or the optimal $\theta^{*}$.

###### Proof.

See Appendix A.3. ∎

######  Theorem 4.

Assume that there are $m$ objective functions $F_{1}(\theta)$,
$F_{2}(\theta)$,…, $F_{m}(\theta)$, and each objective function is
differentiable, $L$-smooth and convex. If
$\cos{<\bar{g},\bar{g}^{\prime}>}\geq\frac{1}{2}$ and
$||\bar{g}||\geq||\bar{g}^{\prime}||$ where $\bar{g}$ is the true gradient for
the average objective $F(\theta)=\frac{1}{m}\sum F_{i}(\theta)$, then FedFV
with step size $\eta\leq\frac{1}{L}$ will converge to either 1) a pareto
stationary point, 2) or the optimal $\theta^{*}$.

###### Proof.

See Appendix A.3. ∎

##  4 Related Work

###  4.1 Fairness in FL

Federated Learning (FL) is first proposed by McMahan et al McMahan et al.
(2017) to collaboratively train a global model distributedly while preserving
data privacyKairouz et al. (2019); Li et al. (2020b); McMahan et al. (2017). A
number of studies focus on collaborative fairness where the server allocates
different models to clients according to their contribution Lyu et al. (2020,
2019); Xu and Lyu (2020). A few other studies address the fairness of a
uniform accuracy distribution across devices Cho et al. (2020); Hu et al.
(2020); Li et al. (2019, 2020a); Mohri et al. (2019). Mohri _et al._ Mohri et
al. (2019) and Li _et al._ Li et al. (2019) propose different federated
objectives AFL and q-FFL to further improve fairness. Hu _et al._ Hu et al.
(2020) observe the competing relation between fairness and robustness to
inflating loss attacks. Abay _et al._ Abay et al. (2020) analyze the potential
causes of bias in FL which leads to unfairness, and they also point out the
negative impact of sample bias due to the party selection. Cho _et al._ Cho et
al. (2020) also show that client selection has an impact on fairness, and they
propose a bandit-based communication-efficient client selection strategy to
overcome biasness. Huang _et al._ Huang et al. (2020) reweight clients
according to their accuracy and numbers of times of being selected to achieve
fairness. They use double momentum to accelerate the convergence. Different
from these, we identify the conflicts among clients to be a potential cause of
unfairness in FL. We mitigate such conflicts by computing a fairer average of
gradients to achieve fair model performance across devices.

###  4.2 Gradient Projection

Gradient projection has been well studied in continual learning to mitigate
the adversary impact of gradient updates to previously learned tasks Chaudhry
et al. (2018); Farajtabar et al. (2019); Guo et al. (2019); Lopez-Paz and
Ranzato (2017). Lopez-Paz _et al._ Lopez-Paz and Ranzato (2017) project
gradient by solving a quadratic programming problem. Chaudhry _et al._
Chaudhry et al. (2018) project gradient onto the normal plane of the average
gradient of previous tasks. Farajtabar _et al._ Farajtabar et al. (2019)
project the current task gradients onto the orthonormal set of previous task
gradients. Yu _et al._ Yu et al. (2020) focus on adversary influence between
task gradients when simultaneously learning multiple tasks. They iteratively
project each task gradient onto the normal plane of conflicting gradients,
which motivates our solution in this paper. To the best of our knowledge, we
are the first to take the adversary gradient interference into consideration
in FL. Specifically, our proposed FedFV method can build a connection between
fairness and conflicting gradients with large differences in the magnitudes,
and we prove its ability to mitigate two gradient conflicts (i.e., internal
and external conflicts) and converge to Pareto stationary solutions.

##  5 Experiments

###  5.1 Experimental Setup

####  5.1.1 Datasets and Models

We evaluate FedFV on three public datasets: CIFAR-10 Krizhevsky (2012),
Fashion MNIST Xiao et al. (2017) and MNIST LeCun et al. (1998). We follow
McMahan et al. (2017) to create non-I.I.D. datasets. For CIFAR-10, we sort all
data records based on their classes, and then split them into 200 shards. We
use 100 clients, and each client randomly picks 2 shards without replacement
so that each has the same data size. We use 200 clients for MNIST and
preprocess MNIST in the same way as CIFAR-10. The local dataset is split into
training and testing data with percentages of 80% and 20%. For Fashion MNIST,
we simply follow the setting in Li et al. (2019). We use a feedforward neural
network with 2 hidden layers on CIFAR-10 and Fashion MNIST. We use a CNN with
2 convolution layers on MNIST.

####  5.1.2 Baselines

We compare with the classical method FedAvg McMahan et al. (2017) and FL
systems that address fairness in FL, including AFL Mohri et al. (2019),
q-FedAvg Li et al. (2019), FedFa Huang et al. (2020) and FedMGDA+ Hu et al.
(2020). We compare with q-fedavg, FedFa and fedmgda+ on all three datasets and
compare with AFL only on Fashion MNIST, because AFL is only suitable for small
networks with dozens of clients.

####  5.1.3 Hyper-parameters

For all experiments, we fix the local epoch $E=1$ and use batchsize
$B_{CIFAR-10|MNIST}\in\\{full,64\\},B_{FashionMNIST}\in\\{full,400\\}$ to run
Stochastic Gradient Descent (SGD) on local datasets with stepsize
$\eta\in\\{0.01,0.1\\}$. We verify the different methods with hyper-parameters
as listed in Table 1. We take the best performance of each method for the
comparison.

Method |  Parameters
---|---
AFL |  $\eta_{\lambda}\in\\{0.01,0.1,0.5\\}$
qFedAvg |  $q\in\\{0.1,0.2,1,2,5,15\\}$
FedMGDA+ |  $\epsilon\in\\{0.01,0.05,0.1,0.5,1\\}$
FedFA |  $(\alpha,\beta)\in\\{(0.5,0.5)\\},(\gamma_{s},\gamma_{c})\in\\{(0.5,0.9)\\}$
FedFV |  $\alpha\in\\{0,0.1,0.2,\frac{1}{3},0.5,\frac{2}{3}\\},\tau\in\\{0,1,3,10\\}$
Table 1: Method specific hyper-parameters.

####  5.1.4 Implementation

All our experiments are implemented on a 64g-MEM Ubuntu 16.04.6 server with 40
Intel(R) Xeon(R) CPU E5-2630 v4 @ 2.20GHz and 4 NVidia(R) 2080Ti GPUs. All
code is implemented in PyTorch version 1.3.1. Please see
https://github.com/WwZzz/easyFL for full details.

###  5.2 Experimental Results

####  5.2.1 Fairness

We first verify the advantage of FedFV in fairness. In Table 2, we list the
the mean, variance, the worst 5% and the best 5% of test accuracy on 100
clients created by splitting CIFAR-10 where $10\%$ of clients are sampled in
each communication round. While reaching a higher mean of test accuracy across
clients, FedFV also yields the lowest variance among all experiments except
for q-FedAvg$|_{q=5.0}$. Note that q-FedAvg sacrifices the best performance of
clients from $69.28\%$ to $63.60\%$ and causes a significant reduction of the
mean of accuracy from $46.85\%$ to $45.25\%$. When reaching a similar variance
to q-FedAvg$|_{q=5.0}$ ($9.59$ versus $9.72$), FedFV improves the mean of
accuracy to $50.42\%$. Further, the worst 5% performance of all the
experiments on FedFV is higher than that of any others, which indicates
FedFV’s effectiveness on protecting the worst performance.

Method | Ave. | Var. | Worst 5% | Best 5%
---|---|---|---|---
FedAvg | 46.85$\pm$0.65 | 12.57$\pm$1.50 | 19.84$\pm$6.55 | 69.28$\pm$1.17
qFedAvg$|_{q=0.1}$ | 47.02$\pm$0.89 | 13.16$\pm$1.84 | 18.72$\pm$6.94 | 70.16$\pm$2.06
qFedAvg$|_{q=0.2}$ | 46.91$\pm$0.90 | 13.09$\pm$1.84 | 18.88$\pm$7.00 | 70.16$\pm$2.10
qFedAvg$|_{q=1.0}$ | 46.79$\pm$0.73 | 11.72$\pm$1.00 | 22.80$\pm$3.39 | 68.00$\pm$1.60
qFedAvg$|_{q=2.0}$ | 46.36$\pm$0.38 | 10.85$\pm$0.76 | 24.64$\pm$2.17 | 66.80$\pm$2.02
qFedAvg$|_{q=5.0}$ | 45.25$\pm$0.42 | 9.59$\pm$0.36 | 26.56$\pm$1.03 | 63.60$\pm$1.13
FedFa | 46.43$\pm$0.56 | 12.79$\pm$1.54 | 19.28$\pm$6.78 | 69.36$\pm$1.40
FedMGDA+$|_{\epsilon=0.01}$ | 45.65$\pm$0.21 | 10.94$\pm$0.87 | 25.12$\pm$2.34 | 67.44$\pm$1.20
FedMGDA+$|_{\epsilon=0.05}$ | 45.58$\pm$0.21 | 10.98$\pm$0.81 | 25.12$\pm$1.87 | 67.76$\pm$2.27
FedMGDA+$|_{\epsilon=0.1}$ | 45.52$\pm$0.17 | 11.32$\pm$0.86 | 24.32$\pm$2.24 | 68.48$\pm$2.68
FedMGDA+$|_{\epsilon=0.5}$ | 45.34$\pm$0.21 | 11.63$\pm$0.69 | 24.00$\pm$1.93 | 68.64$\pm$3.11
FedMGDA+$|_{\epsilon=1.0}$ | 45.34$\pm$0.22 | 11.64$\pm$0.66 | 24.00$\pm$1.93 | 68.64$\pm$3.11
FedFV$|_{\alpha=0.1,\tau=0}$ | 48.57$\pm$0.76 | 10.97$\pm$1.02 | 28.32$\pm$2.01 | 69.76$\pm$2.45
FedFV$|_{\alpha=0.2,\tau=0}$ | 48.54$\pm$0.64 | 10.56$\pm$0.96 | 28.88$\pm$1.92 | 69.20$\pm$2.01
FedFV$|_{\alpha=0.5,\tau=0}$ | 48.14$\pm$0.38 | 10.60$\pm$0.76 | 28.72$\pm$1.92 | 68.80$\pm$1.77
FedFV$|_{\alpha=0.1,\tau=1}$ | 49.34$\pm$0.56 | 10.74$\pm$0.89 | 28.56$\pm$2.54 | 70.24$\pm$1.49
FedFV$|_{\alpha=0.1,\tau=3}$ | 50.00$\pm$0.74 | 10.85$\pm$1.14 | 28.24$\pm$3.03 | 70.96$\pm$1.00
FedFV$|_{\alpha=0.1,\tau=10}$ | 50.42$\pm$0.55 | 9.70$\pm$0.96 | 32.24$\pm$2.10 | 69.68$\pm$2.84
Table 2: The average, the variance, the worst and the best of the test
accuracy of all clients on CIFAR-10. All experiments are running over 2000
rounds with full batch size, learning rate $\eta=0.1$ and local epochs $E=1$.
The reported results are averaged over 5 runs with different random seeds.

Table 3 lists the model test accuracy on each client’s data and the average
and variance of the accuracy on Fashion-MNIST dataset. Different from settings
in Table 2, we use a small network with only 3 clients, and the server selects
all the clients to participate the training process in each communication
round to eliminate the sample bias, which results in the absence of external
conflicts. Therefore, we only use FedFV with
$\alpha\in\\{0,\frac{1}{3},\frac{2}{3}\\},\tau=0$ in this case. In our
observation, even though FedFV with $\alpha=0$ can still find the solution
with the highest mean of test accuracy $80.48\%$, it’s a unfair solution where
the variance is up to $13.76$. Instead, FedFV with $\alpha=\frac{2}{3}$ also
finds a solution with high accuracy $80.28\%$, and simultaneously keep the
variance in a low level. We notice that AFL$|_{\eta_{\lambda}=0.5}$ has the
lowest variance. However, AFL only works for small works where the number of
clients is limited. Despite of this, FedFV$|_{\alpha=\frac{2}{3}}$ outperforms
AFL$|_{\eta_{\lambda}=0.5}$ on all the clients (i.e. shirt, pullover and
T-shirt), as well as the mean of accuracy. Further, when compared to other
methods except AFL, our method has the lowest variance down to 1.77, which
validates the ability of FedFV to find a fair and accurate solution without
the presence of sample bias.

From Table 2 and Table 3, we observe that FedFV is able to ease the negative
impact of conflicting gradients with largely different magnitudes, which
results in a fairer and accurate model with lower variance.

Method | $shirt$ | $pullover$ | $T-shirt$ | Ave. | Var.
---|---|---|---|---|---
FedAvg |  64.26$\pm$1.07 |  87.03$\pm$1.40 |  89.97$\pm$1.30 |  80.42$\pm$0.76 |  11.50$\pm$0.95
AFL$|_{\eta_{\lambda}=0.01}$ |  71.09$\pm$0.91 |  81.89$\pm$1.30 |  84.29$\pm$1.23 |  79.09$\pm$0.73 |  5.76$\pm$0.75
AFL$|_{\eta_{\lambda}=0.1}$ |  76.34$\pm$0.63 |  79.00$\pm$1.04 |  79.06$\pm$1.21 |  78.13$\pm$0.72 |  1.27$\pm$0.66
AFL$|_{\eta_{\lambda}=0.5}$ |  76.57$\pm$0.58 |  78.77$\pm$0.99 |  79.09$\pm$1.15 |  78.14$\pm$0.71 | 1.12$\pm$0.61
qFedAvg$|_{q=5}$ |  71.29$\pm$0.85 |  81.46$\pm$1.08 |  82.86$\pm$0.97 |  78.53$\pm$0.64 |  5.16$\pm$0.69
qFedAvg$|_{q=15}$ |  77.09$\pm$0.88 |  75.40$\pm$2.07 |  60.69$\pm$1.79 |  71.06$\pm$0.81 |  7.46$\pm$0.88
FedFa |  63.80$\pm$1.10 |  86.86$\pm$1.28 |  89.97$\pm$1.30 |  80.21$\pm$0.73 |  11.68$\pm$0.93
FedMGDA+$|_{\epsilon=0.05}$ |  44.63$\pm$1.63 |  57.77$\pm$5.72 | 99.09$\pm$0.23 |  67.16$\pm$1.65 |  23.39$\pm$0.57
FedMGDA+$|_{\epsilon=0.1}$ |  72.26$\pm$3.12 |  79.71$\pm$4.47 |  86.03$\pm$4.27 |  79.33$\pm$1.11 |  6.45$\pm$2.21
FedMGDA+$|_{\epsilon=1.0}$ |  72.46$\pm$3.32 |  79.74$\pm$4.48 |  85.66$\pm$4.88 |  79.29$\pm$1.14 |  6.42$\pm$2.23
FedFV$|_{\alpha=0.0,\tau=0}$ |  61.06$\pm$1.21 | 89.31$\pm$1.38 |  91.06$\pm$1.29 | 80.48$\pm$0.68 |  13.76$\pm$1.02
FedFV$|_{\alpha=\frac{1}{3},\tau=0}$ |  77.09$\pm$1.35 |  80.69$\pm$1.53 |  81.06$\pm$0.53 |  79.61$\pm$0.85 |  1.91$\pm$0.57
FedFV$|_{\alpha=\frac{2}{3},\tau=0}$ | 77.91$\pm$0.80 |  81.46$\pm$1.99 |  81.46$\pm$1.51 |  80.28$\pm$1.09 |  1.77$\pm$0.87
Table 3: Test accuracy on the different clothing classes of Fashion MNIST
dataset. All experiments are running over 200 rounds with full batch size,
learning rate $\eta=0.1$ and local epochs $E=1$. The reported results are
averaged over 5 runs with different random seeds.

####  5.2.2 Accuracy and Efficiency

| FedFV | FedFVRandom | FedFVReverse
---|---|---|---
$Var._{CIFAR10}$ | 13.19$\pm$1.06 | 14.12$\pm$0.55 | 16.28$\pm$0.72
$Var._{FMNIST}$ | 13.76$\pm$1.02 | 20.14$\pm$0.61 | 22.05$\pm$0.68
Table 4: The effects of projecting order to fairness. The results are averaged
over 5 runs with different random seeds.

We further show that FedFV outperforms existing works addressing fairness in
terms of accuracy and efficiency on CIFAR-10 and MNIST. All the methods are
tuned to their best performance. As shown in Figure 2, when we only mitigate
the internal conflicts, which is done by FedFV($\tau=0$), we already converge
faster while keeping a lower variance than the others. In addition, when we
mitigate both the internal and external conflicts by FedFV($\tau\textgreater
0$), there is a substantial improvement in the accuracy and efficiency again.
FedFV outperforms state-of-the-art methods by up to 7.2% on CIFAR-10 and 78%
on MNIST, which verifies the ability of FedFV to reach a higher accuracy with
less communication rounds while still keeping a low variance. Thus, we
demonstrate that FedFV’s advantages in saving communication cost and finding a
better generalization.

(a) CIFAR-10

(b) MNIST

Figure 2: The mean (left) and the variance (right) of test accuracy on all
clients on (a) CIFAR-10, and (b) MNIST. The results are averaged over 5 runs
with different random seeds.

####  5.2.3 Effects of Projecting Order

To verify the effectiveness of our projecting order, we compare FedFV with
another two cases: 1) projecting gradients in a random order of the projection
targets; and 2) projecting gradients in an order of the projection targets
that is reverse to FedFV, as shown in Table 4.

For CIFAR-10, 20% of clients are sampled in each communication round. For
Fashion MNIST, all clients are selected in each communication round. We set
$\alpha=0,\tau=0$ for all groups to confirm the effectiveness of the
projecting order. If we project gradients to the targets in the loss-based
order of FedFV, the variance is the lowest. Projecting gradients to the
targets in a random order is also fairer than in an order reverse to FedFV,
which indicates that the loss-based order used by FedFV helps improve
fairness.

##  6 Conclusions and Future Work

We identified the issue of conflicting gradients with large differences in the
magnitudes, which brings unfairness in FL. To address this issue, we propose
the Federated Fair aVeraging (FedFV) algorithm to mitigate the potential
conflicts among clients before averaging their gradients. In addition, we show
how well FedFV can mitigate the conflicts and how it converges to either a
pareto stationary solution or the optimal on convex problems. Extensive
experiments on a suite of federated datasets confirm that FedFV compares
favorably against state-of-the-art methods in terms of fairness, accuracy and
efficiency. In the future, we plan to build a complete theoretical analysis of
FedFV on mitigating complex external conflicts.

##  Appendix A Proof

###  A.1 Proof of Theorem 1

###### Proof.

For each gradient $g_{i}\in G$, we project $g_{i}$ onto $g_{k}$’s normal plane
in an increasing order with $k$. Thus, we have update rules

$\displaystyle\begin{split}\left\\{\begin{array}[]{ll}g_{i}^{(0)}=g_{i},&k=0\\\
g_{i}^{(k)}=g_{i}^{(k-1)}-\frac{g_{i}^{(k-1)}\cdot
g_{k}}{||g_{k}||^{2}}g_{k},&k=1,2,...,m,k\neq i\\\
g_{i}^{(k)}=g_{i}^{(k-1)},&k=i\end{array}\right.\end{split}$

Therefore, $g_{i}^{(k)}$ never conflicts with $g_{k}$, since all potential
conflicts have been eliminated by the updating rules. We then focus on how
much the final gradient $\bar{g}^{\prime}$ conflicts with each gradient in
$G$. Firstly, the last in the projecting order is $g_{m}$, and the final
gradient $\bar{g}^{\prime}$ will not conflict with it. Then, we focus on the
last but second gradient $g_{m-1}$. Following the update rules, we have

$\begin{split}g_{i}^{(m)}=g_{i}^{(m-1)}-\frac{g_{i}^{(m-1)}\cdot
g_{m}}{||g_{m}||^{2}}g_{m}\\\ \end{split}$ (2)

Let $\phi_{i,j}^{(k)}$ denote the angle between $g_{i}^{(k)}$ and $g_{j}$, and
$\phi_{i,j}=\phi_{i,j}^{(0)}$, then the average of projected gradients is:

$\begin{split}&\bar{g}^{\prime}=\frac{1}{m}\sum_{i}^{m}g_{i}^{(m)}\\\
&=\frac{1}{m}(\sum_{i\neq m}^{m}(g_{i}^{(m-1)}-\frac{g_{i}^{(m-1)}\cdot
g_{m}}{||g_{m}||^{2}}g_{m})+g_{m}^{(m-1)})\\\
&=\frac{1}{m}\sum_{i}^{m}g_{i}^{(m-1)}-\frac{1}{m}\sum_{i\neq
m}^{m}||g_{i}^{(m-1)}||\cos{\phi_{i,m}^{(m-1)}}\frac{g_{m}}{||g_{m}||}\end{split}$
(3)

Since $g_{m-1}\cdot\sum_{i}^{m}g_{i}^{(m-1)}\geq 0$,

$\begin{split}&g_{m-1}\cdot\bar{g}^{\prime}\geq
g_{m-1}\cdot\frac{-1}{m}\sum_{i\neq
m}^{m}||g_{i}^{(m-1)}||\cos{\phi_{i,m}^{(m-1)}}\frac{g_{m}}{||g_{m}||}\\\
&=-\frac{1}{m}\sum_{i\neq
m}^{m}||g_{i}^{(m-1)}||||g_{m-1}||\cos{\phi_{i,m}^{(m-1)}}\cos{\phi_{m-1,m}}\\\
&\geq-\frac{\epsilon^{2}}{m}||g_{m-1}||\sum_{i\neq
m}^{m}||g_{i}^{(m-1)}||\end{split}$ (4)

Similarly, we can compute the conflicts between any gradient $g_{k}\in G$ and
$\bar{g}^{\prime}$ by removing the $g_{i}^{(k)}$ from $\bar{g}^{\prime}$

$\begin{split}g_{k}\cdot\bar{g}^{\prime}&\geq
g_{k}\cdot\sum_{j=k}^{m-1}(\frac{-1}{m}\sum_{i\neq
j+1}^{m}||g_{i}^{(j)}||\cos{\phi_{i,j+1}^{(j)}}\frac{g_{j+1}}{||g_{j+1}||})\\\
&=-\frac{||g_{k}||}{m}\sum_{j=k}^{m-1}\sum_{i\neq
j+1}^{m}||g_{i}^{(j)}||\cos{\phi_{i,j+1}^{(j)}}\cos{\phi_{k,j+1}}\\\
&\geq-\frac{\epsilon^{2}}{m}||g_{k}||\sum_{j=k}^{m-1}\sum_{i\neq
j+1}^{m}||g_{i}^{(j)}||\\\
|e_{k}\cdot\bar{g}^{\prime}|&=|\frac{g_{k}}{||g_{k}||}\cdot\bar{g}^{\prime}|\leq\frac{\epsilon^{2}}{m}\sum_{j=k}^{m-1}\sum_{i\neq
j+1}^{m}||g_{i}^{(j)}||\end{split}$ (5)

Therefore, the later $g_{k}$ serves as the projecting target of others, the
smaller the upper bound of conflicts between $\bar{g}^{\prime}$ and $g_{k}$
is, since $\sum_{j=k}^{m-1}\sum_{i\neq
j+1}^{m}||g_{i}^{(j)}||\leq\sum_{j=k-1}^{m-1}\sum_{i\neq
j+1}^{m}||g_{i}^{(j)}||$. ∎

###  A.2 Proof of Theorem 2

###### Proof.

According to (5), the upper bound of the conflict between any client gradient
$g_{k}$ and the final average gradient $\bar{g}^{\prime}$ can be expressed as

$\begin{split}|g_{k}\cdot\bar{g}^{\prime}|\leq\frac{\epsilon_{2}^{2}}{m}||g_{k}||\sum_{j=k}^{m-1}\sum_{i\neq
j+1}^{m}||g_{i}^{(j)}||\end{split}$ (6)

With the update rules of FedFV, we can infer that

$\begin{split}||g_{k}^{(i)}||^{2}&=||g_{k}^{(i-1)}-||g_{k}^{(i-1)}||\cos{\phi_{k,i}^{(i-1)}}\frac{g_{i}}{||g_{i}||}||^{2}\\\
&=||g_{k}^{(i-1)}||^{2}-2||g_{k}^{(i-1)}||^{2}\cos^{2}{\phi_{k,i}^{(i-1)}}+\\\
&||g_{k}^{(i-1)}||^{2}\cos^{2}{\phi_{k,i}^{(i-1)}}\\\
&=(1-\cos^{2}{\phi_{k,i}^{(i-1)}})||g_{k}^{(i-1)}||^{2}\\\
&\leq(1-\epsilon_{1}^{2})||g_{k}^{(i-1)}||^{2}\end{split}$ (7)

Therefore, the maximum value of gradient conflict is bounded by

$\begin{split}&|g_{k}\cdot\bar{g}^{\prime}|\leq\frac{\epsilon_{2}^{2}}{m}||g_{k}||\sum_{j=k}^{m-1}\sum_{i\neq
j+1}^{m}||g_{i}^{(j)}||\\\
&\leq\frac{\epsilon_{2}^{2}}{m}(\max_{i}{||g_{i}||})\sum_{j=k}^{m-1}\sum_{i\neq
j+1}^{m}||g_{i}^{(0)}||(1-\epsilon_{1}^{2})^{\frac{j}{2}}\\\
&\leq\frac{m-1}{m}\epsilon_{2}^{2}(\max_{i}{||g_{i}||})^{2}\sum_{j=k}^{m-1}(1-\epsilon_{1}^{2})^{\frac{j}{2}}\\\
&=\frac{m-1}{m}(\max_{i}{||g_{i}||})^{2}\frac{\epsilon_{2}^{2}(1-\epsilon_{1}^{2})^{\frac{1}{2}}(1-(1-\epsilon_{1}^{2})^{\frac{m-k}{2}})}{1-(1-\epsilon_{1}^{2})^{\frac{1}{2}}}\\\
&=\frac{m-1}{m}(\max_{i}{||g_{i}||})^{2}f(m,k,\epsilon_{1},\epsilon_{2})\end{split}$
(8)

∎

###  A.3 Proof of the Convergence

###### Proof of Theorem 3.

Let $g_{1}$ and $g_{2}$ be the two types of users’ updates of parameters in
the $t$th communication round. When $F(\theta)$ is $L$-smooth, we have

$F(\theta^{t+1})\leq F(\theta^{t})+\nabla
F(\theta)^{T}(\theta^{t+1}-\theta^{t})+\frac{1}{2}L||\theta^{t+1}-\theta^{t}||_{2}^{2}$
(9)

There are two cases: conflicting gradients exist or otherwise. If there is no
conflict between $g_{1}$ and $g_{2}$, which indicates $g_{1}\cdot g_{2}\geq
0$, FedFV updates as FedAvg does, simply computing an average of the two
gradients and adding it to the parameters to obtain the new parameters
$\theta^{t+1}=\theta^{t}-\eta\bar{g}$, where $g=\frac{1}{2}(g_{1}+g_{2})$.
Therefore, when using step size $\eta\leq\frac{1}{L}$, it will strictly
decrease the objective function $F(\theta)$ Li et al. (2020c). However, if
$g_{1}\cdot g_{2}\textless 0$, FedFV will compute the update as:

$\theta^{t+1}=\theta^{t}-\eta\bar{g}^{\prime}=\theta^{t}-\eta(g_{1}+g_{2}-\frac{g_{1}\cdot
g_{2}}{||g_{1}||^{2}}g_{1}-\frac{g_{1}\cdot g_{2}}{||g_{2}||^{2}}g_{2})$ (10)

Combining with (9), we have:

$\begin{split}&F(\theta^{t+1})\leq
F(\theta^{t})-\eta(g_{1}+g_{2})^{T}(g_{1}+g_{2}-\frac{g_{1}\cdot
g_{2}}{||g_{1}||^{2}}g_{1}-\\\ &\frac{g_{1}\cdot
g_{2}}{||g_{2}||^{2}}g_{2})+\frac{L}{2}\eta^{2}||g_{1}+g_{2}-\frac{g_{1}\cdot
g_{2}}{||g_{1}||^{2}}g_{1}-\frac{g_{1}\cdot
g_{2}}{||g_{2}||^{2}}g_{2}||_{2}^{2}\\\
&=F(\theta^{t})-\eta||g_{1}+g_{1}||^{2}+\eta(g_{1}+g_{2})^{T}(\frac{g_{1}\cdot
g_{2}}{||g_{1}||^{2}}g_{1}\\\ &+\frac{g_{1}\cdot
g_{2}}{||g_{2}||^{2}}g_{2})+\frac{L}{2}\eta^{2}||g_{1}+g_{2}-\frac{g_{1}\cdot
g_{2}}{||g_{1}||^{2}}g_{1}-\frac{g_{1}\cdot
g_{2}}{||g_{2}||^{2}}g_{2}||_{2}^{2}\end{split}$ (11)

Since $g_{1}\cdot g_{2}=||g_{1}||||g_{2}||\cos{\phi_{12}}$, where $\phi_{12}$
denotes the angle between $g_{1}$ and $g_{2}$, after rearranging the items in
this inequality, we have

$\begin{split}F(\theta^{t+1})\leq
F(\theta^{t})-(\eta-\frac{L}{2}\eta^{2})(1-\cos^{2}\phi_{12})(||g_{1}||^{2}+\\\
||g_{2}||^{2})-Lt^{2}(1-\cos^{2}\phi_{12})||g_{1}||||g_{2}||\cos\phi_{12}\end{split}$
(12)

To decrease the objective function, the inequality below should be satisfied

$\begin{split}-(\eta-\frac{L}{2}\eta^{2})(1-\cos^{2}\phi_{12})(||g_{1}||^{2}+||g_{2}||^{2})-L\eta^{2}(1-\\\
\cos^{2}\phi_{12})||g_{1}||||g_{2}||\cos\phi_{12}\leq 0\end{split}$ (13)
$\begin{split}(\eta-\frac{L}{2}\eta^{2})(||g_{1}||^{2}+||g_{2}||^{2})+\frac{L}{2}\eta^{2}2||g_{1}||||g_{2}||\cos\phi_{12}\geq
0\end{split}$ (14)
$\begin{split}(\eta-L\eta^{2})(||g_{1}||^{2}+||g_{2}||^{2})+\frac{L}{2}\eta^{2}(||g_{1}||^{2}+||g_{2}||^{2}+\\\
2||g_{1}||||g_{2}||\cos\phi_{12})\geq 0\end{split}$ (15)
$\begin{split}(\eta-L\eta^{2})(||g_{1}||^{2}+||g_{2}||^{2})+\frac{L}{2}\eta^{2}||g_{1}+g_{2}||^{2}\geq
0\end{split}$ (16)

With $\eta\leq\frac{1}{L}$, we have $\eta-L\eta^{2}\geq 0$, which promises the
decrease of the objective function. On the other hand, we can infer that

$\begin{split}-\eta(1-\frac{L}{2}\eta)=\eta(\frac{L}{2}\eta-1)\leq\eta(\frac{1}{2}-1)=-\frac{\eta}{2}\end{split}$
(17)

Combine (17) with (13), we have:

$\begin{split}&F(\theta^{t+1})\leq
F(\theta^{t})-\frac{\eta}{2}(1-\cos^{2}\phi_{12})(||g_{1}||^{2}+||g_{2}||^{2})\\\
&-\frac{\eta}{2}(1-\cos^{2}\phi_{12})2||g_{1}||||g_{2}||\cos\phi_{12}\\\
&=F(\theta_{t})-\frac{\eta}{2}(1-\cos^{2}\phi_{12})(||g_{1}||^{2}+||g_{2}||^{2}+2||g_{1}||||g_{2}||\\\
&\cos\phi_{12})\\\
&=F(\theta_{t})-\frac{\eta}{2}(1-\cos^{2}\phi_{12})||g_{1}+g_{2}||^{2}\\\
&=F(\theta_{t})-\frac{\eta}{2}(1-\cos^{2}\phi_{12})||g||^{2}\end{split}$ (18)

Therefore, the objective function will always degrade unless 1)$||g||=0$,
which indicates it will reach the optimal $\theta^{*}$,
2)$\cos{\phi_{12}}=-1$, then we can suppose $g_{1}=-\beta g_{2}$ where
$\beta\textgreater 0$, and we will get

$\begin{split}\frac{1}{1+\beta}g_{1}+\frac{\beta}{1+\beta}g_{2}=0\end{split}$
(19)

So $\theta^{t}$ is a pareto stationary point as defined below:

######  Definition 4.

For smooth criteria $F_{k}(\theta)(1\leq k\leq K)$, $\theta^{0}$ is called a
Pareto-stationary iff there exists some convex combination of the gradietns
$\nabla F_{k}(\theta^{0})$ that equals zero Désidéri (2012).

∎

###### Proof of Theorem 4.

With the assumption that $F(\theta)$ is $L$-smooth, we can get the inequality:

$F(\theta^{t+1})\leq F(\theta^{t})+\nabla
F(\theta)^{T}(\theta^{t+1}-\theta^{t})+\frac{1}{2}L||\theta^{t+1}-\theta^{t}||_{2}^{2}.$
(20)

Similar to the proof for Theorem 2, if there exist conflicts, then

$\begin{split}F(\theta^{t+1})&\leq
F(\theta^{t})-\eta\bar{g}\cdot\bar{g}^{\prime}+\frac{1}{2}L\eta^{2}||\bar{g}^{\prime}||^{2}\\\
&\leq
F(\theta^{t})-\frac{\eta}{2}||\bar{g}||||\bar{g}^{\prime}||+\frac{L}{2}\eta^{2}||\bar{g}^{\prime}||^{2}\\\
&\leq
F(\theta^{t})-\frac{\eta}{2}||\bar{g}||||\bar{g}^{\prime}||+\frac{L}{2}\eta^{2}||\bar{g}||||\bar{g}^{\prime}||\\\
&\leq
F(\theta^{t})+(-\frac{\eta}{2}+\frac{L}{2}\eta^{2})||\bar{g}||||\bar{g}^{\prime}||\\\
\end{split}$ (21)

Since $(-\frac{\eta}{2}+\frac{L}{2}\eta^{2})\leq 0$ with
$\eta\leq\frac{1}{L}$, the average objective will always degrade if we
repeatedly apply the update rules of FedFV unless 1)$||\bar{g}||=0$, which
indicates it will finally reach the optimal $\theta^{*}$, 2)
$||\bar{g}^{\prime}||=0$, which indicates that all pairs of conflicting
gradients have a cosine similarity of $-1$, leading to an existence of convex
combination $p$ which satisfies $\sum p_{i}g_{i}=0$, since we can easily
choose a pair of conflicting gradients $g_{i}$ and $g_{j}$ as the proof of
theorem 3 does and then set the weights of the rest zero. ∎

## References

  * Abay et al. (2020) Annie Abay, Yi Zhou, Nathalie Baracaldo, Shashank Rajamoni, Ebube Chuba, and Heiko Ludwig.  Mitigating bias in federated learning.  arXiv e-prints, page arXiv:2012.02447, December 2020.
  * Chaudhry et al. (2018) Arslan Chaudhry, Marc’Aurelio Ranzato, Marcus Rohrbach, and Mohamed Elhoseiny.  Efficient lifelong learning with A-GEM.  CoRR, abs/1812.00420, 2018.
  * Cho et al. (2020) Yae Jee Cho, Samarth Gupta, Gauri Joshi, and Osman Yağan.  Bandit-based communication-efficient client selection strategies for federated learning.  arXiv e-prints, page arXiv:2012.08009, 2020.
  * Désidéri (2012) Jean-Antoine Désidéri.  Multiple-gradient descent algorithm (mgda) for multiobjective optimization.  Comptes Rendus Mathematique, 350(5):313 – 318, 2012.
  * Farajtabar et al. (2019) Mehrdad Farajtabar, Navid Azizan, Alex Mott, and Ang Li.  Orthogonal gradient descent for continual learning.  CoRR, abs/1910.07104, 2019.
  * Guo et al. (2019) Yunhui Guo, Mingrui Liu, Tianbao Yang, and Tajana Rosing.  Learning with long-term remembering: Following the lead of mixed stochastic gradient.  CoRR, abs/1909.11763, 2019.
  * Hu et al. (2020) Zeou Hu, Kiarash Shaloudegi, Guojun Zhang, and Yaoliang Yu.  Fedmgda+: Federated learning meets multi-objective optimization.  arXiv e-prints, page arXiv:2006.11489, 2020.
  * Huang et al. (2020) Wei Huang, Tianrui Li, Dexian Wang, Shengdong Du, and Junbo Zhang.  Fairness and accuracy in federated learning.  arXiv e-prints, page arXiv:2012.10069, 2020.
  * Kairouz et al. (2019) Peter Kairouz, H. Brendan McMahan, Brendan Avent, Aurélien Bellet, Mehdi Bennis, and et al.  Advances and open problems in federated learning.  CoRR, abs/1912.04977, 2019.
  * Krizhevsky (2012) Alex Krizhevsky.  Learning multiple layers of features from tiny images.  University of Toronto, 05 2012.
  * LeCun et al. (1998) Yann LeCun, Léon Bottou, Yoshua Bengio, and Patrick Haffner.  Gradient-based learning applied to document recognition.  Proceedings of the IEEE, 86(11):2278–2324, 1998.
  * Li et al. (2019) Tian Li, Maziar Sanjabi, and Virginia Smith.  Fair resource allocation in federated learning.  CoRR, abs/1905.10497, 2019.
  * Li et al. (2020a) Tian Li, Shengyuan Hu, Ahmad Beirami, and Virginia Smith.  Federated multi-task learning for competing constraints.  arXiv e-prints, page arXiv:2012.04221, 2020.
  * Li et al. (2020b) Tian Li, Anit Kumar Sahu, Manzil Zaheer, Maziar Sanjabi, Ameet Talwalkar, and Virginia Smith.  Federated optimization in heterogeneous networks.  arXiv e-prints, page arXiv:1812.06127, 2020.
  * Li et al. (2020c) Xiang Li, Kaixuan Huang, Wenhao Yang, Shusen Wang, and Zhihua Zhang.  On the convergence of fedavg on non-iid data.  arXiv e-prints, page arXiv:1907.02189, 2020.
  * Lopez-Paz and Ranzato (2017) David Lopez-Paz and Marc’Aurelio Ranzato.  Gradient episodic memory for continuum learning.  CoRR, abs/1706.08840, 2017.
  * Lyu et al. (2019) Lingjuan Lyu, Jiangshan Yu, Karthik Nandakumar, Yitong Li, Xingjun Ma, and Jiong Jin.  Towards fair and decentralized privacy-preserving deep learning with blockchain.  CoRR, abs/1906.01167, 2019.
  * Lyu et al. (2020) Lingjuan Lyu, Xinyi Xu, and Qian Wang.  Collaborative fairness in federated learning.  arXiv e-prints, page arXiv:2008.12161, 2020.
  * McMahan et al. (2017) Brendan McMahan, Eider Moore, Daniel Ramage, Seth Hampson, and Blaise Aguera y Arcas.  Communication-Efficient Learning of Deep Networks from Decentralized Data.  In Aarti Singh and Jerry Zhu, editors, Proceedings of the 20th International Conference on Artificial Intelligence and Statistics, volume 54 of Proceedings of Machine Learning Research, pages 1273–1282, Fort Lauderdale, FL, USA, 20–22 Apr 2017. PMLR.
  * Mohri et al. (2019) Mehryar Mohri, Gary Sivek, and Ananda Theertha Suresh.  Agnostic federated learning.  CoRR, abs/1902.00146, 2019.
  * Xiao et al. (2017) Han Xiao, Kashif Rasul, and Roland Vollgraf.  Fashion-mnist: a novel image dataset for benchmarking machine learning algorithms.  CoRR, abs/1708.07747, 2017.
  * Xu and Lyu (2020) Xinyi Xu and Lingjuan Lyu.  Towards building a robust and fair federated learning system.  arXiv e-prints, page arXiv:2011.10464, 2020.
  * Yu et al. (2020) Tianhe Yu, Saurabh Kumar, Abhishek Gupta, Sergey Levine, Karol Hausman, and Chelsea Finn.  Gradient Surgery for Multi-Task Learning.  arXiv e-prints, page arXiv:2001.06782, January 2020.
  * Zhao et al. (2018) Yue Zhao, Meng Li, Liangzhen Lai, Naveen Suda, Damon Civin, and Vikas Chandra.  Federated Learning with Non-IID Data.  arXiv e-prints, page arXiv:1806.00582, June 2018.