# Fairness-aware Federated Minimax Optimization with Convergence Guarantee
††thanks: This work was supported by a grant from the NSFC/RGC Joint Research
Scheme sponsored by the Research Grants Council of the Hong Kong Special
Administrative Region, China and National Natural Science Foundation of China
(Project No. N_HKUST656/22).

Gerry Windiarto Mohamad Dunda  Electronic and Computer Engineering
The Hong Kong University of Science and Technology
Clear Water Bay, Hong Kong
<EMAIL_ADDRESS>Shenghui Song  Electronic and Computer Engineering
The Hong Kong University of Science and Technology
Clear Water Bay, Hong Kong
<EMAIL_ADDRESS>
###### Abstract

Federated learning (FL) has garnered considerable attention due to its
privacy-preserving feature. Nonetheless, the lack of freedom in managing user
data can lead to group fairness issues, where models are biased towards
sensitive factors such as race or gender. To tackle this issue, this paper
proposes a novel algorithm, fair federated averaging with augmented Lagrangian
method (FFALM), designed explicitly to address group fairness issues in FL.
Specifically, we impose a fairness constraint on the training objective and
solve the minimax reformulation of the constrained optimization problem. Then,
we derive the theoretical upper bound for the convergence rate of FFALM. The
effectiveness of FFALM in improving fairness is shown empirically on CelebA
and UTKFace datasets in the presence of severe statistical heterogeneity.

###### Index Terms:

federated learning, group fairness, convergence rate, augmented Lagrangian

##  I Introduction

Federated learning (FL) [1] is a distributed machine learning approach that
enables model training on potentially sensitive data from different entities
without the necessity for data sharing. This technique is promising in diverse
domains such as computer vision (CV) as it can facilitate training of models
on a large-scale, diverse set of data while preserving data privacy. However,
a direct implementation of existing federated algorithms may violate group
fairness [2], which refers to the equitable treatment of different groups in a
population. Group fairness is required by law such as in Europe [3], enforcing
that the decision making by predictive models does not exhibit bias towards
any particular sensitive group, such as race or gender. For example, an AI
model used in a hiring process may have been trained on historical data that
reflects biased hiring patterns, leading to discriminatory outcomes for
underrepresented groups in the workforce. There are more examples [4] that
further motivate raising awareness in training fair deep learning models.

The sources of group unfairness or bias mainly come from dataset, which may
reflect measurement or historical bias from the annotators, and the training
algorithm, which may learn unwanted biased features from such dataset. The
aforementioned sources also induce statistical heterogeneity. Under severe
heterogeneity, the trained model may have poor task performance and biased
model. Finding FL algorithms that are fair and robust to statistical
heterogeneity or non-identical and independently distributed (non-iid) data is
an arduous task, and currently it is an open problem [5].

###  I-A Contributions

Since we have little control on the clients data in FL, efforts to mitigate
bias from the federated algorithm are important. Therefore, the focus of this
paper is to improve group fairness in binary classification tasks, involving
binary sensitive attributes, which is common in fairness literature [6]. In
particular, we propose a new federated algorithm that is effective in reducing
bias while maintaining similar task performance as existing federated
algorithms. The main contributions are summarized below.

  * •

We propose a fairness-aware algorithm, fair federated averaging with augmented
Lagrangian method (FFALM). Firstly, we formulate a constrained minimization
problem on the global loss function satisfying a fairness metric. Inspired by
augmented Lagrangian method [7], we solve the problem by leveraging the local
training as a sub-solver to find the optimal model parameters given dual
iterates. Then, the dual iterates are locally updated given the updated model
parameter, and they are aggregated using weighted average.

  * •

We propose a theoretical upper bound for the convergence rate of FFALM over a
nonconvex-strongly-concave objective function, which is
$\mathcal{O}(\frac{1}{T^{2/3}})$.

  * •

We empirically assess the proposed method on publicly accessible CV datasets
(CelebA and UTKFace) containing sensitive attributes (gender and skin color
respectively) based on two common fairness metrics: demographic parity
difference (DPD) and equal opportunity difference (EOD). Experimental results
show that FFALM improves DPD by $10\%$ and EOD by $14\%$ in the attractiveness
prediction task, and DPD by $6\%$ and EOD by $5\%$ in the youth prediction
task compared to FedAvg under severe data heterogeneity.

##  II Related Work

There have been some engaging results in tackling the fairness issues in deep-
learning models. In the following, we categorize some prior related works
based on how the training is conducted, either centralized or federated.

###  II-A Ensuring fairness in centralized learning

In centralized learning, it is not uncommon to modify the training framework
to achieve a suitable degree of group fairness. The authors of [8]
decorrelated the input images and the protected group attributes by using
adversarial training. Knowledge transfer techniques and multi-classifiers can
also be adopted as a debiasing method [9]. Augmenting each image sample with
its perturbed version generated from generative models can potentially reduce
biases as well [10]. The aforementioned works require additional components to
the model, thus increasing the computation cost. This might not be suitable
for FL. A possible alternative is to alter the loss function to take into
account group fairness. The authors of [11] introduced a loss function
obtained from the upper bound of the Lagrangian based on a constrained
optimization formulation from a finite hypothesis space perspective.

###  II-B Ensuring fairness in FL

Some prior works considered group fairness in FL. Due to system constraints,
most innovations came from the objective formulation to include fairness or
more information exchange between the server and the clients. The example for
the latter is FairFed [12], where the client coefficients are adaptively
adjusted during the aggregation phase based on the deviation of each client’s
fairness metric from the global average. Along the line of the objective
formulation, FCFL [13] proposed a two-stage optimization to solve a multi-
objective optimization with fairness constraints, which demands more
communication rounds. FPFL [14] utilized differential multipliers optimization
method to solve main objective by taking into account fairness as a
constraint, which is similar to this work. FedFB [15] adjusted the weight of
the local loss function for each sensitive group during the aggregation phase.
The objective formulation of FPFL, however, is not smooth, which may hinder
the convergence of gradient-based learning. Moreover, the theoretical
convergence guarantee is missing for the aforementioned works.

##  III Preliminaries

In this section, we introduce some mathematical notations and group fairness
notions. After that, we briefly describe the framework of minimax FL.

###  III-A Notations

Throughout this paper, we primarily focus on supervised binary classification
tasks with binary sensitive attributes. The dataset is denoted as
$\mathcal{D}=X\times Y\times S$ with size $|\mathcal{D}|$ constituting of an
input image $X$, a label $Y=\\{0,1\\}$, and a sensitive attribute
$S=\\{0,1\\}$. We slightly abuse the notation of $\mathcal{D}$, $X$, $Y$, and
$S$ to represent both the set and the distribution. The datasets can also be
partitioned based on sensitive attributes, $\mathcal{D}^{s_{0}}=X\times
Y\times S_{0}$ and $\mathcal{D}^{s_{1}}=X\times Y\times S_{1}$.

Some mathematical notations are stated as follows. $[N]$ denotes
$\\{1,2,...,N\\}$, $\|.\|$ represents the $\ell_{2}$-norm, and $\mathbbm{1}$
denotes the indicator function. We use $\mathcal{W}\subseteq\mathbb{R}^{d}$
and $\Lambda$ to represent the parameter spaces of the model $w$ and an
additional learnable training parameter $\lambda$, respectively. Denote
$f_{w}:X\rightarrow Y$ as a deep-learning model parameterized by $w$, taking
$X$ as an input, and outputting the predicted label $Y$, and denote
$q_{w}:X\rightarrow\mathbb{R}^{2}$ as the logits of the model, where the first
element corresponds to $Y=0$ and the second element corresponds to $Y=1$.

###  III-B Group Fairness Metrics

To evaluate the group fairness performance of a machine learning model $f$,
there are various notions based on how likely the model predicts a favorable
outcome ($\hat{Y}=f(X)=1$) for each group. Demographic parity (DP) [16] is
commonly used for assessing the fairness of the model. $f$ satisfies DP if the
model prediction of favorable label is independent of $S$, i.e.,

$\displaystyle\mathbb{E}_{X|S=0}[f(X)=1]=\mathbb{E}_{X|S=1}[f(X)=1].$ (1)

Another way to define the notion of group fairness is accuracy parity (AP)
[17]. To satisfy this notion, $f$ conforms to the following equality

$\displaystyle\mathbb{E}_{\mathcal{D}^{s_{0}}}[\mathbbm{1}_{f(X)\neq
Y}]=\mathbb{E}_{\mathcal{D}^{s_{1}}}[\mathbbm{1}_{f(X)\neq Y}].$ (2)

In some use cases where the preference of users belonging to a sensitive group
is considered, it is amenable to adopt equal opportunity (EO) [18] of positive
outcomes for each sensitive attribute as a fairness notion, mathematically
written as

$\displaystyle\mathbb{E}_{X|S=0,Y=1}[f(X)=1]=\mathbb{E}_{X|S=1,Y=1}[f(X)=1].$
(3)

In practice, it is difficult to achieve perfect fairness imposed by
aforementioned fairness notions. To measure how close $f$ is to satisfy DP, we
employ demographic parity difference metric $(\Delta_{DP})$ on favorable
label, which is defined as

$\displaystyle\Delta_{DP}=|\mathbb{E}_{X|S=0}[\mathbbm{1}_{f(X)=1}]-\mathbb{E}_{X|S=1}[\mathbbm{1}_{f(X)=1}]|.$
(4)

Similarly, the closeness measure to satisfy EO condition, equal opportunity
difference ($\Delta_{EO}$) is defined as

$\displaystyle\Delta_{EO}=|\mathbb{E}_{X|Y=1,S=0}[\mathbbm{1}_{f(X)=1}]-\mathbb{E}_{X|Y=1,S=1}[\mathbbm{1}_{f(X)=1}]|.$
(5)

These two closeness metrics are commonly used for assessing group fairness in
machine learning models. Since only samples rather than the true data
distribution are available, the metrics are estimated using the samples.

Figure 1: The minimax FL framework. In each communication round, there are
four steps in FL training: 1. Broadcasting phase 2. Local training phase 3.
Client-to-server communication phase 4. Aggregation phase.

###  III-C Minimax FL Framework

For a FL system with $N$ clients and one server, the goal is to train a global
deep learning model $f_{w}$ on each client dataset $\mathcal{D}_{i}$
$(i\in[N])$ without sharing their datasets. In addition, an additional
learning parameter (e.g. Lagrangian dual) $\lambda\in\Lambda$ that aids the
training can be exchanged between the server and clients, and processed on the
clients and the server. Its generic procedure consists of four phases.
Firstly, the clients receive the global model from the server (broadcasting
phase) and train the model on their own dataset (local training phase). After
that, the clients send their model update and the dual parameter update to the
server (client-to-server communication phase), and the server aggregates the
received updates from each participating client to get an updated model
(aggregation phase). This process is repeated until convergence or a specified
communication round, as illustrated in Fig. 1.

During local training phase, each client aims to minimize their local risk
function $F_{i}(w,\lambda)$ to update their model. At the same time, it
maximizes the local risk function to update $\lambda$, in contrast to the
conventional FL framework. For any regularization-based federated algorithms,
the explicit formulation for the true local risk function of the $i-$th client
represented by a loss function $l(f_{w}(x),y)=l(x,y;w)$ and a regularization
function $g$ is given by

$\displaystyle
F_{i}(w,\lambda):=\mathbb{E}_{(x_{j},y_{j})\sim\mathcal{D}_{i}}[l(x_{j},y_{j};w)]+g(\mathcal{D}_{i};\lambda,w),$
(6)

and the corresponding empirical risk function is given by

$\displaystyle
F_{i,S}(w,\lambda):=\frac{1}{|\mathcal{D}_{i}|}\sum_{(x_{j},y_{j})\sim\mathcal{D}_{i}}l(x_{j},y_{j};w)+g(\mathcal{D}_{i};\lambda,w).$
(7)

The ultimate goal is to solve the following minimax objective of the true
global risk function $F(w,\lambda)$

$\displaystyle\min_{w}\max_{\lambda}F(w,\lambda):=\sum_{i=1}^{N}p_{i}F_{i}(w,\lambda),$
(8)

where $p_{i}$ is the client coefficient with $\sum_{i=1}^{N}p_{i}=1$ and
$p_{i}\in[0,1]$. In FedAvg, the coefficient is set to the proportion of the
samples from each client. Since the clients only have access to samples rather
than the data distribution, the objective is replaced with the global
empirical risk function defined as

$\displaystyle F_{S}(w,\lambda):=\sum_{i=1}^{N}p_{i}F_{i,S}(w,\lambda).$ (9)

##  IV FFALM

We first introduce the problem formulation for FL with group fairness
constraints. Subsequently, we describe the proposed algorithm to achieve the
objective. Lastly, we offer the convergence rate of the proposed algorithm.

###  IV-A Solving Group Fairness Issue

####  IV-A1 Problem formulation

The objective of this work is to ensure group fairness on the FL-trained
binary classification model. We tackle the problem by enforcing fairness
during the local training. Specifically, the local training aims to minimize
the local risk function while satisfying a notion of fairness. The strategy is
to reformulate the local risk function as a sum of the main objective, which
is related to the task performance, and the term related to fairness
constraint, weighted with a learnable $\lambda$. In this way, we can use
minimax FL framework as described in the previous section.

One of the essential requirements of minimax FL framework is to have a smooth
local risk function. Since the indicator functions appearing in the fairness
notions in Section III-B are not differentiable functions with respect to the
model parameters, we need to replace all indicator functions with their
corresponding surrogate continuous functions. In the case of AP, the choice of
such continuous functions is readily available, which is cross entropy loss
$CE(y,q_{w}(x))=-\log\sigma(q_{w}(x)_{y})$, where
$\sigma(x)=\frac{1}{1+e^{-x}}$ is a sigmoid function. This is similar to how
0-1 loss ($\mathbb{E}_{\mathcal{D}}[\mathbbm{1}_{\hat{Y}\neq Y}]$) in the
basic gradient-based learning can be replaced with cross-entropy loss [19].
The only difference in the formulation of AP is that it is conditioned on each
sensitive attribute.

To this end, we write the objective of the local training of the $i$-th client
as a constrained optimization

$\displaystyle\min_{w}L_{S}(w,\mathcal{D}_{i})\textrm{ s.t.
}\mu(w,\mathcal{D}_{i}^{s_{0}})=\mu(w,\mathcal{D}_{i}^{s_{1}}),$ (10)

where
$L_{S}(w,\mathcal{D}):=\frac{1}{|\mathcal{D}|}\sum_{(x_{j},y_{j})\sim\mathcal{D}}l(x_{j},y_{j};w)$
and
$\mu(w,\mathcal{D}):=\frac{1}{|\mathcal{D}|}\sum_{(x,y)\in\mathcal{D}}CE(y,q_{w}(x))$.
It is worth mentioning that we also estimate (2) from the samples, as shown in
the definition of $\mu(w,\mathcal{D})$.

####  IV-A2 Local training phase

This problem can be approximately solved by following similar techniques from
the augmented Lagrangian approach [7] by treating the constraint as a soft
constraint. Specifically, it seeks a saddle solution of an augmented
Lagrangian function $\mathcal{L}_{S}$ parameterized by a suitable choice of
penalty coefficient ($\beta$)

$\displaystyle\mathcal{L}_{S}(w_{t},\lambda_{t-1},\mathcal{D}_{i})$
$\displaystyle:=L_{S}(w_{t},\mathcal{D}_{i})+\frac{\beta}{2}\Delta_{\mu}(w_{t},\mathcal{D}_{i})^{2}$
(11) $\displaystyle+\lambda_{t-1}\Delta_{\mu}(w_{t},\mathcal{D}_{i}),$

where
$\Delta_{\mu}(w_{t},\mathcal{D}_{i}):=\mu(w_{t},\mathcal{D}_{i}^{s_{0}})-\mu(w_{t},\mathcal{D}_{i}^{s_{1}})$
by introducing a sub-optimizer $O$ to seek $w_{t}$ such that
$\|\nabla_{w}\mathcal{L}_{S}(w_{t},\beta,\mathcal{D}_{i})\|$ is sufficiently
small. Afterwards, $\lambda_{t-1}$ is updated to close the infeasibility gap,
and the process is repeated. If the algorithm converges to the solution
$(w^{*},\lambda^{*})$ of (11) that satisfies second-order sufficient
conditions [7], $w^{*}$ is the global solution to (10).

Translating this view into FL, we can assign the sub-optimizer $O$ to the
local training and the iteration index $t$ to the communication round. Hence,
we can formulate the local training of $i$-th client as a two-stage process

$\displaystyle w_{i,t}$
$\displaystyle=\min_{w}F_{i,S}(w,\lambda_{t-1})=\min_{w}\mathcal{L}_{S}(w,\lambda_{t-1},\mathcal{D}_{i})$
(12) $\displaystyle\lambda_{i,t}$
$\displaystyle=\lambda_{t-1}+\eta_{\lambda,t}\Delta_{\mu}(w_{t},\mathcal{D}_{i}).$
(13)

Note that in the original augmented Lagrangian method, $\eta_{\lambda,t}$ is
set to $\beta$. As shown later in the experiment results section, this
proposed two-stage optimization gives more competitive results in terms of
fairness performance.

Stochastic gradient descent is used to solve (12), similar to FedAvg.
Specifically, the $i$-th client computes the stochastic gradient of
$\mathcal{L}^{(t,k)}_{i}=\mathcal{L}_{S}(w_{i}^{(t,k)},\lambda_{i,t},\mathcal{B}_{i})$
at communication round $t$ and local iteration $k$ from its batch samples
$\mathcal{B}_{i}$ sampled from its local distribution $\mathcal{D}_{i}$ as

$\displaystyle\nabla_{w}\mathcal{L}_{i}^{(t,k)}$
$\displaystyle=\nabla_{w}\biggl{(}L_{S}(w_{i}^{(t,k)},\mathcal{B}_{i})+\lambda\Delta_{\mu}(w_{i}^{(t,k)},\mathcal{B}_{i})$
(14)
$\displaystyle+\frac{\beta}{2}\Delta_{\mu}(w_{i}^{(t,k)},\mathcal{B}_{i})^{2}\biggr{)}.$

####  IV-A3 Aggregation phase

The server receives model updates, $w_{i,t}$, as well as the dual updates,
$\lambda_{i,t}$, from the clients. Following FedAvg, the received dual update
from each client is aggregated by weighted average with the same client
coefficient ($p_{i}$) as model aggregation

$\displaystyle
w_{t}=\sum_{i=1}^{N}p_{i}w_{i,t}\quad\text{and}\quad\lambda_{t}=\sum_{i=1}^{N}p_{i}\lambda_{i,t}.$
(15)

The proposed algorithm is summarized in Algorithm 1.

Algorithm 1 FFALM Algorithm

1:Inputs: $N$, $\\{\mathcal{D}_{i}\\}_{i=1}^{N}$, $\beta$, $\eta_{w,t}$, the
number of local iteration $E$, and the maximum communication round $T$.

2:Randomly initialize the global model $w_{0}$ and set $\lambda_{0}=0$ on the
server side

3:for $t=1$ to $T$ do

4: Broadcast $w_{t-1}$ and $\lambda_{t-1}$ to all clients

5: for each $i\in[N]$ do

6: $w_{i}^{(t,0)}\leftarrow w_{t-1}$

7: for $k=1$ to $E$ do

8: Randomly sample the batch $\mathcal{B}_{i}$ from $\mathcal{D}_{i}$

9: Compute $\nabla_{w}\mathcal{L}_{i}^{(t,k-1)}$ from (14)

10: $w_{i}^{(t,k)}\leftarrow
w_{i}^{(t,k-1)}-\eta_{w,t}\nabla_{w}\mathcal{L}_{i}^{(t,k-1)}$

11: end for

12: Compute $\lambda_{i,t}$ from (13)

13: Send $w_{i,t}=w_{i}^{(t,E)}$ and $\lambda_{i,t}$ to the server

14: end for

15: Update $\eta_{w,t}$ using a LR scheduler

16: Aggregation phase to obtain $w_{t}$ and $\lambda_{t}$ following (15)

17:end for

###  IV-B Theoretical Convergence Guarantee

The proposed algorithm can be viewed as solving a minimax problem with
$F_{S}(w,\lambda)=\sum_{i=1}^{N}p_{i}\mathcal{L}_{S}(w,\lambda_{t-1},\mathcal{D}_{i})$.
We provide the upper bound of the convergence rate based on how close the
empirical primal risk function
$R_{S}(w_{t}):=\max_{\lambda}\mathcal{L}(w_{t},\lambda)$ is to the optimal.
Before presenting the result, we list several definitions and key assumptions.

######  Definition 1.

Define a function $h:\mathcal{W}\times\Lambda\rightarrow\mathbb{R}$.
$h(\cdot,\cdot)$ is $L$-smooth if it is continuously differentiable and there
exists a constant $L>0$ such that for any $w,w^{\prime}\in\mathcal{W}$,
$\lambda,\lambda^{\prime}\in\mathbb{R}$, and $\xi\in\mathcal{D}$,

$\displaystyle\left\|\begin{pmatrix}\nabla_{w}h(w,\lambda;\xi)-\nabla_{w}h(w^{\prime},\lambda^{\prime};\xi)\\\
\nabla_{\lambda}h(w,\lambda;\xi)-\nabla_{\lambda}h(w^{\prime},\lambda^{\prime};\xi)\end{pmatrix}\right\|\leq
L\left\|\begin{pmatrix}w-w^{\prime}\\\
\lambda-\lambda^{\prime}\end{pmatrix}\right\|.$

######  Definition 2.

$h(w,\cdot)$ is $\rho$-strongly convex if for all $w\in\mathcal{W}$ and
$\lambda,\lambda^{\prime}\in\Lambda$, $h(w,\lambda)\geq
h(w,\lambda^{\prime})+\langle\nabla_{\lambda}h(w,\lambda^{\prime}),\lambda-\lambda^{\prime}\rangle+\frac{\rho}{2}\|\lambda-\lambda^{\prime}\|^{2}$.

######  Definition 3.

$h(w,\cdot)$ is $\rho$-strongly concave if $-h(w,\cdot)$ is $\rho$-strongly
convex.

######  Assumption 1 (Bounded Variance).

For randomly drawn batch samples $\xi$ and for all $i\in[N]$, the gradients
$\nabla_{w}F_{i,S}(w,\lambda;\xi)$ and
$\nabla_{\lambda}F_{i,S}(w,\lambda;\xi)$ have bounded variances $B_{w}$ and
$B_{\lambda}$ respectively. If
$g_{i,w}(w,\lambda|\xi):=\nabla_{w}F_{i,S}(w,\lambda;\xi)$ is the unbiased
local estimator of the gradient,
$\mathbb{E}_{\xi}[\|g_{i,w}(w,\lambda|\xi)-\nabla_{w}F_{i,S}(w,\lambda)\|^{2}]\leq
B_{w}^{2}$, and the case for $\lambda$ is similar but bounded by
$B_{\lambda}^{2}$.

######  Assumption 2 (Bounded Dissimilarity).

There exist $\beta\geq 1$ and $\kappa\geq 0$ such that

$\displaystyle\sum_{i=1}^{N}p_{i}\|\nabla_{w}F_{i,S}(w,\lambda)\|^{2}\leq\beta^{2}\|\sum_{i=1}^{N}p_{i}\nabla_{w}F_{i,S}(w,\lambda)\|^{2}+\kappa^{2}.$
(16)

######  Definition 4.

A function $h(\cdot,\lambda)$ satisfies the PL condition if for all $\lambda$,
there exists a constant $\mu>0$ such that, for any $w\in\mathcal{W}$,
$\frac{1}{2}\|\nabla
h(w)\|^{2}\geq\mu(h(w)-\min_{w^{\prime}\in\mathcal{W}}h(w^{\prime}))$.

Definition 1, Assumption 1 and Assumption 2 are commonly used for federated
learning [20]. Assumption 2 is satisfied when the gradient clipping method is
employed. Lastly, the PL-condition of $\mathcal{L}(\cdot,\lambda)$ is shown to
hold on a large class of neural networks [21].

For simplicity, we assume full participation and the same number of local
iterations for each client. The minimum empirical primal risk is
$R^{*}_{S}=\min_{w}R_{S}(w)$. The upper bound of the convergence rate of FFALM
is given by the following theorem.

######  Theorem 1.

Define $\kappa=\frac{L}{\mu}$. Let $\eta_{w,t}=\frac{2\beta^{2}}{\mu t}$ and
$\eta_{\lambda,t}=\mathcal{O}(\frac{1}{t^{2/3}})$. Given that Assumption 1 and
Assumption 2 hold, each $F_{i,S}(w,\lambda)$ is $L$-smooth, each
$F_{i,S}(\cdot,\lambda)$ satisfies $\mu$-PL condition, and each
$F_{i,S}(w,\cdot)$ is $\rho$-strongly concave, we have

$\displaystyle\mathbb{E}R_{S}(w_{T+1})-R_{S}^{*}=\mathcal{O}\biggl{(}\frac{\Gamma+B_{w}^{2}+B_{\lambda}^{2}}{T^{2/3}}\biggr{)},$

after $T$ communication rounds, where
$\Gamma:=F_{S}^{*}-\sum_{i=1}^{N}p_{i}F_{i,S}^{*}$,
$F_{S}^{*}:=\min_{w}\max_{\lambda}F_{S}(w,\lambda)$ and
$F_{i,S}^{*}:=\min_{w}\max_{\lambda}F_{i,S}(w,\lambda)$.

###### Proof:

See Appendix A. ∎

Note that $\Gamma$ quantifies statistical heterogeneity of the FL system. In
the case of strong non-iid, the saddle solution of the global risk function is
significantly different from the weighted sum of each saddle local risks.

##  V Empirical Study

(a)

(b)

Figure 2: Learning curves on validation set for (a) CelebA dataset and (b)
UTKFace dataset

In this section, we evaluate the effectiveness of FFALM based on three
important performance metrics: the prediction accuracy, DPD, and EOD on real-
world datasets. We provide the results and comparison with other existing FL
algorithms.

###  V-A Datasets and Tasks

We want to investigate how the classification model trained in FL, which is
designed to further increase the prediction accuracy on specific domain,
influences the fairness performance. In particular, two datasets from CV
domain are used in this study: CelebA [22], and UTKFace [23], and ResNet-18
[24] models are used for both datasets. The task of CelebA dataset is a binary
classification for predicting attractiveness in images with gender as the
sensitive attribute, and the task of UTKFace dataset is to determine whether
the age of a person from an image is above or below 20 with skin color as the
sensitive attribute. The favorable labels for CelebA and UTKFace are
attractive and age below 20 respectively. Each dataset is split into three
categories: training set, validation set, and testing set.

###  V-B FL Setting

There are 10 clients participating in FL training. We synthetically simulate
statistical heterogeneity by introducing label skews, which can be implemented
using Dirichlet distribution parameterized by $\alpha$ on the proportion of
samples for a given class and client on a given centralized samples [25]. The
case of severe data heterogeneity is investigated in this experiment by
setting $\alpha=0.3$. The experiment for each instance is repeated 10 times
with different seeds. The FL training ends after 70 communication rounds.

###  V-C Baselines

The following are the baselines used for the comparison study.

  1. 1.

FedAvg. It is the universal baseline in FL which trains the model locally
without considering fairness and aggregates all model updates by weighted
average.

  2. 2.

FairFed [12]. The server receives the local DP metrics, and based on them and
the global trend, the server adjusts the value of $p_{i}$ adaptively before
averaging the model updates.

  3. 3.

FPFL [14]. It enforces fairness by solving the constrained optimization on the
sample loss function $L_{S}$ with two constraints. These constraints ensure
that the absolute difference between the overall loss and the loss of each
sensitive group
($\delta_{\mu_{0},i}=[L(w^{t},\mathcal{D}_{i})-\mu(w^{t},\mathcal{D}_{i}^{s_{0}})]_{+}$
and
$\delta_{\mu_{1},i}=[L(w^{t},\mathcal{D}_{i})-\mu(w^{t},\mathcal{D}_{i}^{s_{1}})]_{+}$,
where $[x]_{+}:=\max(0,x)$) does not deviate from a specific threshold. We set
this threshold to be zero. Hence, we reformulate it as a local constrained
optimization with

$\displaystyle g(\mathcal{D}_{i};\lambda,w)$
$\displaystyle=\lambda_{0}\delta_{\mu_{0},i}+\lambda_{1}\delta_{\mu_{1},i}+\frac{\beta}{2}(\delta_{\mu_{0},i}^{2}+\delta_{\mu_{1},i}^{2}).$

TABLE I: Comprehensive list of hyperparameter values used in the experiments on CelebA and UTKFace datasets for baselines and FFALM. Hyperparameters |  Algorithms |  Values
---|---|---
Batch size |  all |  128
Gradient clipping on $w$ |  all |  1.0
LR of $w$ decay step size |  all |  50
LR or $w$ decay step factor |  all |  0.5
$\eta_{w,0}$ |  all |  0.05
$\beta$ |  FairFed |  0.5
$\beta$ |  FPFL |  5.0
$\eta_{\lambda,t}$ |  FPFL |  0.5
$b$ |  FFALM |  1.05
$\beta$ |  FFALM |  2.0
$\eta_{\lambda,0}$ |  FFALM |  2.0
$\lambda_{0}$ |  FFALM and FPFL |  0.0

###  V-D Implementation Details

The hyperparameters used in the experiment for FFALM and all baselines are
shown in Table LABEL:table:hyperparameters. Following the setup from augmented
Lagrangian method, we slowly increase the learning rate of $\lambda$ by a
factor $b$ per communication round for FFALM. All algorithms have the same
learning rate scheduler of $w$ which is realized by a constant-step decay
factor.

TABLE II: Comparison of the performance of the proposed algorithm across
baselines on CelebA and UTKFace datasets. $\uparrow$ indicates the larger the
value the better and $\downarrow$ indicates the smaller the value the better.

Algorithm | Acc $\uparrow(\%)$ |  $\Delta_{DP}$ $\downarrow(\%)$ |  $\Delta_{EO}$ $\downarrow(\%)$
---|---|---|---
CelebA
FedAvg | $74.66$ | $39.54$ | $21.11$
FairFed | $73.27$ | $37.36$ | $19.50$
FPFL | $74.58$ | $30.87$ | $9.71$
FFALM | $74.06$ | $28.92$ | $6.82$
UTKFace
FedAvg | $86.22$ | $13.07$ | $19.12$
FairFed | $79.54$ | $10.80$ | $17.92$
FPFL | $86.73$ | $7.97$ | $15.85$
FFALM | $86.02$ | $6.60$ | $13.87$

###  V-E Results

The training curves of all algorithms are shown in Figure 2. The curves show
how accuracy, DPD, and EOD change as the communication round increases. It can
be observed that FFALM shows similar improvement in accuracy as FedAvg while
improving the fairness performance.

The experimental results evaluated on the testing set are presented in Table
III. For celebA dataset, FFALM improves DPD by almost $11\%$ and EOD by
roughly $14\%$ compared to FedAvg. FFALM outperforms other baselines in
fairness performance with minimal accuracy loss. For UTKFace dataset, FFALM
improves DP difference performance by about $6\%$, and reduces EO gap by $6\%$
compared to FedAvg. The overall fairness improvement is also apparent on FFALM
in UTKFace dataset compared with different baselines.

##  VI Conclusion

In this paper, we proposed FFALM, an FL algorithm based on augmented
Lagrangian framework to handle group fairness issues. It leveraged accuracy
parity constraint for smooth loss formulation of minimax FL framework. It was
shown that the theoretical convergence rate of FFALM is
$\mathcal{O}(\frac{1}{T^{2/3}})$. Experiment results on CelebA and UTKFace
datasets demonstrated the effectiveness of the proposed algorithm in improving
fairness with negligible accuracy drop under severe statistical heterogeneity.

##  Appendix A Proof of Theorem 1

We introduce and prove several lemmas before proving Theorem 4.2. The proof is
similar to Theorem 3 in [26] except we have additional lemmas to cover the
local training phases and statistical heterogeneity that exist in FL.

######  Lemma 1 ([27]).

Assume $L$-smoothness for $F_{i,S}$ with $i\in[N]$ and $F_{i,S}(w,\cdot)$ is
$\rho$-strongly concave. Let $\lambda$ be bounded. Then the function
$R_{S}(w)$ is $L+\frac{L^{2}}{\rho}$-smooth and $\nabla
R_{S}(w)=\nabla_{w}F_{S}(w,\hat{\lambda}_{S}(w))$, where
$\hat{\lambda}_{S}(w)=\arg\max_{\lambda}F_{S}(w,\lambda)$. Moreover,
$\hat{\lambda}_{S}(w)$ is $\frac{L}{\rho}$-Lipschitz continuous.

######  Lemma 2.

Assume $F_{i,S}$ with $i\in[N]$ is $L$-smooth and satisfies and all
$F_{i,S}(w,\cdot)$ with $i\in[N]$ is $\rho$-strongly concave. Assume that the
FL system data dissimilarity can be bounded, i.e.
$\sum_{i=1}^{N}p_{i}\|\nabla_{w}F_{i,S}(w,\lambda)\|^{2}\leq\beta^{2}\|\sum_{i=1}^{N}p_{i}\nabla_{w}F_{i,S}(w,\lambda)\|^{2}+\kappa^{2}$.
With $\Gamma:=F^{*}_{S}-\sum_{i=1}^{N}p_{i}F^{*}_{i,S}$, where
$F_{S}^{*}:=\min_{w}\max_{\lambda}F_{S}(w,\lambda)$ and
$F_{i,S}^{*}:=\min_{w}\max_{\lambda}F_{i,S}(w,\lambda)$, we have

$\displaystyle\|\nabla
R_{S}(w)\|^{2}\geq\frac{2\mu}{\beta^{2}}(R_{S}(w)-\min_{w}R_{S}(w)+\Gamma-\kappa^{2}).$
(18)

###### Proof:

$\displaystyle\|\nabla_{w}F_{S}(w,\hat{\lambda}_{S}(w))\|^{2}=\|\sum_{i=1}^{N}p_{i}\nabla_{w}F_{i,S}(w,\hat{\lambda}_{S}(w))\|^{2}$
$\displaystyle\geq\frac{\sum_{i=1}^{N}p_{i}\|\nabla_{w}F_{i,S}(w,\hat{\lambda}_{S}(w))\|^{2}-\kappa^{2}}{\beta^{2}}$
$\displaystyle\geq\frac{\sum_{i=1}^{N}2\mu
p_{i}(F_{i,S}(w,\hat{\lambda}_{S}(w))-\min_{w^{\prime}\in\mathcal{W}}F_{i,S}(w^{\prime},\hat{\lambda}_{S}(w)))-\kappa^{2}}{\beta^{2}}$
$\displaystyle\geq\frac{\sum_{i=1}^{N}2\mu
p_{i}(F_{i,S}(w,\hat{\lambda}_{S}(w))-\min_{w^{\prime}\in\mathcal{W}}F_{i,S}(w^{\prime},\hat{\lambda}_{S}(w^{\prime})))-\kappa^{2}}{\beta^{2}}$
$\displaystyle=\frac{2\mu(F_{S}(w,\hat{\lambda}_{S}(w))-\sum_{i=1}^{N}p_{i}\min_{w^{\prime}\in\mathcal{W}}F_{i,S}(w^{\prime},\hat{\lambda}_{S}(w^{\prime})))-\kappa^{2}}{\beta^{2}}$
$\displaystyle=\frac{2\mu}{\beta^{2}}(F_{S}(w,\hat{\lambda}_{S}(w))-\min_{w}F_{S}(w,\hat{\lambda}_{S}(w))+\Gamma-\kappa^{2}).$

The first inequality is obtained by Jensen inequality, and the second
inequality comes from the $L$-smoothness property. Since
$\nabla_{w}R_{S}(w)=\nabla_{w}F_{S}(w,\hat{\lambda}_{S}(w))$ by Lemma 1, the
proof is complete. ∎

######  Lemma 3.

Assume Assumption 2 holds for each $F_{i,S}$ and each $F_{i,S}$ with $i\in[N]$
is $L$-smooth. If the number of local iterations is $E$, we have

$\displaystyle\|\nabla
R_{S}(w)-\sum_{i=1}^{N}p_{i}\sum_{k=1}^{E}\nabla_{w}F_{i,S}(w_{i}^{(t,k)},\lambda_{t})\|^{2}\leq
L^{2}[\||\hat{\lambda}_{S}(w_{t})-\lambda_{t}\|^{2}]$
$\displaystyle+\frac{2L^{2}E}{\mu}(R_{S}(w_{t})-R_{S}^{*}+\Gamma).$ (19)

###### Proof:

Using Jensen inequality and $L$-smoothness, we get

$\displaystyle\|\nabla
R_{S}(w)-\sum_{i=1}^{N}p_{i}\sum_{k=1}^{E}\nabla_{w}F_{i,S}(w_{i}^{(t,k)},\lambda_{t})\|^{2}$
$\displaystyle\leq\sum_{i=1}^{N}p_{i}\sum_{k=1}^{E}[\|\nabla_{w}F_{i,S}(w_{i}^{(t,0)},\hat{\lambda}_{S}(w_{t}))-\nabla_{w}F_{i,S}(w_{i}^{(t,k)},\lambda_{t})\|^{2}]$
$\displaystyle\leq\sum_{i=1}^{N}p_{i}L^{2}\sum_{k=1}^{E}([\||\hat{\lambda}_{S}(w_{t})-\lambda_{t}\|^{2}]+[\|w^{(t,0)}-w_{i}^{(t,k)}\|^{2}].$
$\displaystyle\leq
L^{2}E\||\hat{\lambda}_{S}(w_{t})-\lambda_{t}\|^{2}+L^{2}E\sum_{i=1}^{N}p_{i}\|w^{(t,0)}-w_{i}^{*}\|^{2}$
$\displaystyle\leq
L^{2}E\||\hat{\lambda}_{S}(w_{t})-\lambda_{t}\|^{2}+\frac{2L^{2}E}{\mu}\sum_{i=1}^{N}p_{i}(F_{i,S}(w_{t},\hat{\lambda}_{S}(w_{t}))$
$\displaystyle-F_{i,S}(w_{i}^{*},\hat{\lambda}_{S}(w_{i}^{*})))$
$\displaystyle\leq
L^{2}E\||\hat{\lambda}_{S}(w_{t})-\lambda_{t}\|^{2}+\frac{2L^{2}E}{\mu}(R_{S}(w_{t})-R_{S}^{*}+\Gamma)$
(20)

The second inequality is from $L$-smoothness. The fourth inequality is from
Appendix A (PL $\rightarrow$ QG) in [28]. ∎

######  Lemma 4.

Rewrite the combined primal update (local update + aggregation) of the global
model at round $t$ ($w_{t}$) as
$w_{t+1}=w_{t}-\eta_{w,t}\sum_{i=1}^{N}p_{i}\sum_{k=1}^{E}g_{w}(w^{(t,k)}_{i},\lambda_{t};\xi_{i}^{(t,k)})$,
where $E$ is the number of local iterations. Assume that Assumption 1 and the
conditions from Lemma 3, 2, and 1 hold. We have

$\displaystyle\mathbb{E}_{t}[R_{S}(w_{t+1})-R^{*}_{S}]$
$\displaystyle\leq(1-\frac{\mu}{\beta^{2}}\eta_{w,t}+\frac{L^{2}E\eta_{w,t}}{\mu})(R_{S}(w_{t})-R^{*}_{S})$
(21)
$\displaystyle+\frac{L^{2}\eta_{w,t}}{2}\mathbb{E}[\||\hat{\lambda}_{S}(w_{t})-\lambda_{t}\|^{2}]$
$\displaystyle+\frac{\eta^{\prime}_{w,t}}{2}E^{2}B_{w}^{2}+\frac{\mu}{\beta^{2}}\eta^{\prime}_{w,t}\kappa^{2}$
$\displaystyle+(\frac{L^{2}E\eta_{w,t}}{\mu}-\frac{\mu}{\beta^{2}}\eta_{w,t})\Gamma.$

###### Proof:

For simplicity, denote $w_{t}=w^{(t,0)}_{i}=w^{(t,0)}$. We start from the
smoothness of $R_{S}$.

$\displaystyle R_{S}(w_{t+1})-R^{*}_{S}\leq
R_{S}(w_{t})-R^{*}_{S}+\langle\nabla_{w}R_{S}(w_{t}),w_{t+1}-w_{t}\rangle$
$\displaystyle+\frac{L+\frac{L^{2}}{\rho}}{2}\|w_{t+1}-w_{t}\|^{2}$
$\displaystyle\leq
R_{S}(w_{t})-R^{*}_{S}-\eta_{w,t}\langle\nabla_{w}R_{S}(w_{t}),\sum_{i=1}^{N}p_{i}\sum_{k=1}^{E}g_{w}(w^{(t,k)}_{i},\lambda_{t};\xi_{i}^{(t,k)})\rangle$
$\displaystyle+\frac{(L+\frac{L^{2}}{\rho})\eta_{w,t}^{2}}{2}\left\|\sum_{i=1}^{N}p_{i}\sum_{k=1}^{E}g_{w}(w^{(t,k)}_{i},\lambda_{t};\xi_{i}^{(t,k)})\right\|^{2}.$

We denote $\mathbb{E}_{t}$ as the conditional expectation over all samples
$\\{\xi_{i}^{(t,k)}\\}$ with given $w_{t}$ and $\lambda_{t}$. Taking this
conditional expectation on both sides we get,

$\displaystyle\mathbb{E}_{t}[R_{S}(w_{t+1})-R^{*}_{S}]\leq
R_{S}(w_{t})-R^{*}_{S}$
$\displaystyle-\eta_{w,t}\langle\nabla_{w}R_{S}(w_{t}),\sum_{i=1}^{N}p_{i}\sum_{k=1}^{E}\nabla
F_{i,S}(w^{(t,k)},\lambda_{t})\rangle$
$\displaystyle+\frac{(L+\frac{L^{2}}{\rho})\eta_{w,t}^{2}}{2}\mathbb{E}_{t}\|\sum_{i=1}^{N}p_{i}\sum_{k=1}^{E}g_{w}(w^{(t,k)}_{i},\lambda_{t};\xi_{i}^{(t,k)})$
$\displaystyle-\sum_{i=1}^{N}p_{i}\sum_{k=1}^{E}\nabla
F_{i,S}(w^{(t,k)}_{i},\lambda_{t})+\sum_{i=1}^{N}p_{i}\sum_{k=1}^{E}\nabla
F_{i,S}(w^{(t,k)}_{i},\lambda_{t})\|^{2}$ $\displaystyle\leq
R_{S}(w_{t})-R^{*}_{S}-\eta_{w,t}\langle\nabla_{w}R_{S}(w_{t}),\sum_{i=1}^{N}p_{i}\sum_{k=1}^{E}\nabla
F_{i,S}(w_{i}^{(t,k)},\lambda_{t})\rangle$
$\displaystyle+\frac{(L+\frac{L^{2}}{\rho})\eta_{w,t}^{2}}{2}\left\|\sum_{i=1}^{N}p_{i}\sum_{k=1}^{E}\nabla
F_{i,S}(w^{(t,k)}_{i},\lambda_{t})\right\|^{2}$
$\displaystyle+\frac{(L+\frac{L^{2}}{\rho})\eta_{w,t}^{2}}{2}E^{2}B_{w}^{2}.$
(22)

The last inequality is obtained by zeroing uncorrelated gradient estimators
between two different clients or different local iterations, Jensen
inequality, and
$\mathbb{E}_{t}\left\|\sum_{i=1}^{N}p_{i}\sum_{k=1}^{E}(g_{w}(w^{(t,k)}_{i},\lambda_{t};\xi_{i}^{(t,k)})-\nabla
F_{i,S}(w^{(t,k)}_{i},\lambda_{t}))\right\|^{2}\leq E^{2}B_{w}^{2}$. Further
simplifying the term with $\eta_{w,t}\leq\frac{1}{(L+L^{2}/\rho)}$ and
$\eta^{\prime}_{w,t}=(L+L^{2}/\rho)\eta_{w,t}^{2}$, and using Lemma 2, we
obtain

$\displaystyle\mathbb{E}_{t}[R_{S}(w_{t+1})-R^{*}_{S}]=R_{S}(w_{t})-R^{*}_{S}$
$\displaystyle+\frac{\eta_{w,t}}{2}\|\nabla_{w}R_{S}(w)-\sum_{i=1}^{N}p_{i}\sum_{k=1}^{E}\nabla
F_{i,S}(w^{(t,k)},\lambda_{t})\|^{2}$
$\displaystyle-\frac{\eta_{w,t}}{2}\|\nabla_{w}R_{S}(w)\|^{2}\
+\frac{\eta_{w,t}}{2}E^{2}B_{w}^{2}$
$\displaystyle\leq(1-\frac{\mu}{\beta^{2}}\eta_{w,t})(R_{S}(w_{t})-R^{*}_{S})+\frac{\eta_{w,t}}{2}\|\nabla_{w}R_{S}(w_{t})$
$\displaystyle-\sum_{i=1}^{N}p_{i}\sum_{k=1}^{E}\nabla
F_{i,S}(w^{(t,k)},\lambda_{t})\|^{2}$
$\displaystyle+\frac{\eta^{\prime}_{w,t}}{2}E^{2}B_{w}^{2}-\frac{\mu}{\beta^{2}}\eta_{w,t}\Gamma+\frac{\mu}{\beta^{2}}\eta^{\prime}_{w,t}\kappa^{2}$
$\displaystyle\leq(1-\frac{\mu}{\beta^{2}}\eta_{w,t}+\frac{L^{2}E\eta_{w,t}}{\mu})(R_{S}(w_{t})-R^{*}_{S})+\frac{L^{2}\eta_{w,t}}{2}\mathbb{E}[\||\hat{\lambda}_{S}(w_{t})-\lambda_{t}\|^{2}]$
$\displaystyle+\frac{\eta^{\prime}_{w,t}}{2}E^{2}B_{w}^{2}+\frac{\mu}{\beta^{2}}\eta^{\prime}_{w,t}\kappa^{2}+(\frac{L^{2}E\eta_{w,t}}{\mu}-\frac{\mu}{\beta^{2}}\eta_{w,t})\Gamma.$

The second inequality is from Lemma 3. ∎

######  Lemma 5.

Rewrite the combined dual update as
$\lambda_{t+1}=\lambda_{t}+\eta_{\lambda,t}\sum_{i=1}^{N}p_{i}g_{\lambda}(w^{(t,k)}_{i},\lambda_{t};\xi_{i}^{(t,k)})$.
Assuming the conditions from Lemma 1 and Assumption 1 hold, we have

$\displaystyle\mathbb{E}_{t}\|\lambda_{t+1}-\hat{\lambda}_{S}(w_{t+1})\|^{2}\leq(1+\frac{1}{\epsilon})\frac{L^{2}\eta_{w,t}^{2}}{\rho^{2}}(E^{2}B_{w}^{2}$
$\displaystyle+2L^{2}\mathbb{E}[\||\hat{\lambda}_{S}(w_{t})-\lambda_{t}\|^{2}])$
$\displaystyle+(1+\frac{1}{\epsilon})(\frac{4L^{2}\eta_{w,t}}{\rho^{2}}+\frac{4L^{4}E\eta_{w,t}^{2}}{\rho^{2}\mu})(R_{S}(w_{t})-R_{S}(w^{*}))$
$\displaystyle+(1+\frac{1}{\epsilon})\frac{4L^{4}E\eta_{w,t}^{2}}{\rho^{2}\mu}\Gamma$
$\displaystyle+(1+\epsilon)(1-\rho\eta_{\lambda,t})\|\lambda_{t}-\hat{\lambda}_{S}(w_{t})\|^{2}+(1+\epsilon)\eta_{\lambda,t}^{2}B_{\lambda}^{2}$
(23)

for any $\epsilon>0$.

###### Proof:

By Young’s inequality, for any $\epsilon>0$, we have

$\displaystyle\|\lambda_{t+1}-\hat{\lambda}_{S}(w_{t+1})\|^{2}$
$\displaystyle\leq(1+\epsilon)\|\lambda_{t+1}-\hat{\lambda}_{S}(w_{t})\|^{2}$
(24)
$\displaystyle+(1+\frac{1}{\epsilon})\|\hat{\lambda}_{S}(w_{t})-\hat{\lambda}_{S}(w_{t+1})\|^{2}.$

For the second term, using the fact that $\hat{\lambda}_{S}(w)$ is
$\frac{L}{\rho}$-Lipschitz (Lemma 1) and applying the expectation to get

$\displaystyle\mathbb{E}_{t}[\|\hat{\lambda}_{S}(w_{t+1})-\hat{\lambda}_{S}(w_{t})\|^{2}]\leq\frac{L^{2}}{\rho^{2}}\mathbb{E}_{t}[\|w_{t+1}-w_{t}\|^{2}]$
$\displaystyle=\frac{L^{2}\eta_{w,t}^{2}}{\rho^{2}}\mathbb{E}_{t}[\|\sum_{i=1}^{N}p_{i}\sum_{k=1}^{E}g_{w}(w^{(t,k)}_{i},\lambda_{t};\xi_{i}^{(t,k)})\|^{2}]$
$\displaystyle=\frac{L^{2}\eta_{w,t}^{2}}{\rho^{2}}\mathbb{E}_{t}\bigl{\|}\sum_{i=1}^{N}p_{i}\sum_{k=1}^{E}(g_{w}(w^{(t,k)}_{i}-\nabla_{w}F_{i,S}(w^{(t,k)}_{i},\lambda_{t}))$
$\displaystyle+\sum_{i=1}^{N}p_{i}\sum_{k=1}^{E}\nabla_{w}F_{i,S}(w^{(t,k)}_{i},\lambda_{t})\|^{2}$
$\displaystyle\leq\frac{L^{2}\eta_{w,t}^{2}}{\rho^{2}}(EB_{w}^{2}+\bigl{\|}\sum_{i=1}^{N}p_{i}\sum_{k=1}^{E}\nabla_{w}F_{i,S}(w^{(t,k)}_{i},\lambda_{t})\|^{2})$
$\displaystyle\leq\frac{2L^{2}\eta_{w,t}^{2}}{\rho^{2}}\left\|\nabla
R_{S}(w_{t})-\sum_{i=1}^{N}\frac{p_{i}}{E}\sum_{k=1}^{E}\nabla
F_{i}(w^{(t,k)}_{i},\lambda_{t})\right\|^{2}$
$\displaystyle+\frac{2L^{2}\eta_{w,t}^{2}}{\rho^{2}}\|\nabla
R_{S}(w_{t})\|^{2}+\frac{L^{2}\eta_{w,t}^{2}EB_{w}^{2}}{\rho^{2}}$
$\displaystyle\leq\frac{2L^{4}\eta_{w,t}^{2}}{\rho^{2}}\mathbb{E}[\||\hat{\lambda}_{S}(w_{t})-\lambda_{t}\|^{2}]+\frac{4L^{4}E\eta_{w,t}^{2}}{\rho^{2}\mu}(R_{S}(w_{t})-R_{S}^{*}-\Gamma)$
$\displaystyle+\frac{2L^{2}\eta_{w,t}^{2}}{\rho^{2}}\|\nabla
R_{S}(w_{t})\|^{2}+\frac{L^{2}\eta_{w,t}^{2}EB_{w}^{2}}{\rho^{2}}$
$\displaystyle\leq\frac{2L^{4}\eta_{w,t}^{2}}{\rho^{2}}\mathbb{E}[\||\hat{\lambda}_{S}(w_{t})-\lambda_{t}\|^{2}]+\frac{4L^{4}E\eta_{w,t}^{2}}{\rho^{2}\mu}(R_{S}(w_{t})-R_{S}^{*}-\Gamma)$
$\displaystyle+\frac{L^{2}\eta_{w,t}^{2}B_{w}^{2}}{\rho^{2}}+\frac{4L^{2}\eta_{w,t}}{\rho^{2}}(R_{S}(w_{t})-R_{S}(w^{*})).$
(25)

The first inequality is obtained from Assumption 1, the second inequality is
obtained from triangle inequality and independence, and the third inequality
is obtained from Lemma 3.

For the first term, we have

$\displaystyle\mathbb{E}_{t}[\|\lambda_{t+1}-\hat{\lambda}_{S}(w_{t})\|^{2}]$
$\displaystyle\leq\mathbb{E}_{t}[\|\lambda_{t}+\eta_{\lambda,t}\sum_{i=1}^{N}p_{i}g_{\lambda}(w^{(t,E)}_{i},\lambda_{t};\xi_{i}^{(t,k)})-\hat{\lambda}_{S}(w_{t})\|^{2}]$
$\displaystyle\leq\|\lambda_{t}-\hat{\lambda}_{S}(w_{t})\|^{2}+2\eta_{\lambda,t}\mathbb{E}_{t}[\langle\lambda_{t}-\hat{\lambda}_{S}(w_{t}),$
$\displaystyle\eta_{\lambda,t}\sum_{i=1}^{N}p_{i}g_{\lambda}(w^{(t,E)}_{i},\lambda_{t};\xi_{i}^{(t,k)})\rangle]$
$\displaystyle+\eta_{\lambda,t}^{2}\mathbb{E}_{t}[\|\sum_{i=1}^{N}p_{i}(g_{\lambda}(w^{(t,E)}_{i},\lambda_{t};\xi_{i}^{(t,k)})-\nabla_{\lambda}F_{i,S}(w_{i}^{(t,E)},\lambda_{t}))$
$\displaystyle+\nabla_{\lambda}F_{i,S}(w_{i}^{(t,E)},\lambda_{t})\|^{2}]$
$\displaystyle\leq(1-\rho\eta_{\lambda,t})\|\lambda_{t}-\hat{\lambda}_{S}(w_{t})\|^{2}+\eta_{\lambda,t}^{2}\|\nabla_{\lambda}F_{S}(w^{(t,E)})\|^{2}+\eta_{\lambda,t}^{2}B_{\lambda}^{2}$
$\displaystyle+2\eta_{\lambda,t}\sum_{i=1}^{N}p_{i}(F_{i,S}(w_{i}^{(t,E)},\hat{\lambda}_{S}(w_{t}))-F_{i,S}(w_{i}^{(t,E)},\lambda_{t}))$
$\displaystyle\leq(1-\rho\eta_{\lambda,t})\|\lambda_{t}-\hat{\lambda}_{S}(w_{t})\|^{2}-\frac{\eta_{\lambda,t}}{L}\|\nabla_{\lambda}F_{S}(w_{i}^{(t,E)},\lambda_{t})\|^{2}$
$\displaystyle+\eta_{\lambda,t}^{2}\|\nabla_{\lambda}F_{S}(w^{(t,E)})\|^{2}+\eta_{\lambda,t}^{2}B_{\lambda}^{2}$
$\displaystyle\leq(1-\rho\eta_{\lambda,t})\|\lambda_{t}-\hat{\lambda}_{S}(w_{t})\|^{2}+\eta_{\lambda,t}^{2}B_{\lambda}^{2}.$
(26)

The fourth inequality is from the $\rho$-strongly concave property of
$F_{i,S}(w,\cdot)$. The last inequality is obtained by
$\eta_{\lambda,t}\leq\frac{1}{L}$ By substituting (25) and (26) into (24), we
get the final result. ∎

######  Lemma 6.

Define $a_{t}=\mathbb{E}_{t}[R_{S}(w_{t})-R^{*}_{S}]$ and
$b_{t}=\mathbb{E}_{t}\|\lambda_{t}-\hat{\lambda}_{S}(w_{t})\|^{2}$. For any
non-increasing sequence $\\{\nu_{t}>0\\}$ and any $\epsilon>0$, we have the
following relations,

$\displaystyle a_{t+1}+\nu_{t+1}b_{t+1}\leq
k_{1,t}a_{t}+k_{2,t}\nu_{t}b_{t}+c_{t},$

where

$\displaystyle k_{1,t}$
$\displaystyle=1-\frac{\mu\eta_{w,t}}{\beta^{2}}+\frac{L^{2}E\eta_{w,t}}{\mu}+\nu_{t}(1+\frac{1}{\epsilon})(\frac{4L^{2}\eta_{w,t}}{\rho^{2}}+\frac{4L^{4}E\eta_{w,t}^{2}}{\rho^{2}\mu})$
$\displaystyle k_{2,t}$
$\displaystyle=\frac{L^{2}\eta_{w,t}}{2\nu_{t}}+(1+\epsilon)(1-\rho\eta_{\lambda,t})$
$\displaystyle c_{t}$
$\displaystyle=\frac{\eta_{w,t}}{2}\zeta_{t}+\frac{\eta^{\prime}_{w,t}}{2}E^{2}B_{w}^{2}+\frac{\mu}{\beta^{2}}\eta^{\prime}_{w,t}\kappa^{2}$
$\displaystyle+\nu_{t}(1+\frac{1}{\epsilon})\frac{L^{2}\eta_{w,t}^{2}}{\rho^{2}}EB_{w}^{2}+\nu_{t}(1+\epsilon)\eta_{\lambda,t}^{2}B_{\lambda}^{2}$
$\displaystyle+(\frac{L^{2}E\eta_{w,t}}{\mu}-\frac{\mu}{\beta^{2}}\eta_{w,t}+\nu_{t}(1+\frac{1}{\epsilon})\frac{4L^{4}E\eta_{w,t}^{2}}{\rho^{2}\mu})\Gamma.$

###### Proof:

The result follows from combining Lemma 4 and Lemma 5. ∎

We can proceed to prove Theorem 1. From Lemma 6, by choosing
$\epsilon=\frac{\rho\eta_{\lambda,t}}{2(1-\rho\eta_{\lambda,t})}$, we have

$\displaystyle k_{2,t}$
$\displaystyle=\frac{L^{2}\eta_{w,t}}{2\nu_{t}}+1-\frac{\rho\eta_{\lambda,t}}{2}.$

Then, set $\nu_{t}=\frac{2L^{2}\eta_{w,t}}{\rho\eta_{\lambda,t}}$. we have

$\displaystyle k_{1,t}$
$\displaystyle=1-\frac{\mu\eta_{w,t}}{\beta^{2}}+\frac{L^{2}E\eta_{w,t}}{\mu}+\frac{16L^{2}\eta_{w,t}^{2}}{\rho^{4}\eta_{\lambda,t}^{2}}(1+\frac{L^{2}E\eta_{w,t}}{\mu})$
(27) $\displaystyle k_{2,t}$ $\displaystyle=1-\frac{\rho\eta_{\lambda,t}}{4}.$
(28)

Let
$\eta^{2}_{\lambda,t}\geq\frac{16L^{2}\eta_{w,t}}{\rho^{4}(\frac{\mu}{2\beta^{2}}-\frac{L^{2}E}{\mu})}(1+\frac{L^{2}E\eta_{w,t}}{\mu})$
and $\eta_{\lambda,t}\geq\frac{2\mu}{\beta^{2}\rho}\eta_{w,t}$. We get

$\displaystyle
a_{t+1}+\nu_{t+1}b_{t+1}\leq(1-\frac{\mu}{2\beta^{2}}\eta_{w,t})(a_{t}+\nu_{t}b_{t})+\frac{\eta^{\prime}_{w,t}}{2}E^{2}B_{w}^{2}+\frac{\mu}{\beta^{2}}\eta^{\prime}_{w,t}\kappa^{2}$
$\displaystyle+\frac{4L^{2}\eta_{w,t}^{3}}{\rho^{3}\eta^{2}_{\lambda,t}}EB_{w}^{2}+\frac{L^{2}\eta_{\lambda,t}\eta_{w,t}B_{\lambda}^{2}}{\rho(1-\rho\eta_{\lambda,t})}$
$\displaystyle+(\frac{L^{2}E\eta_{w,t}}{\mu}-\frac{\mu}{\beta^{2}}\eta_{w,t}+\frac{16L^{6}E\eta_{w,t}^{3}}{\rho^{4}\mu\eta_{\lambda,t}^{2}})\Gamma.$
(29)

Choose $\eta_{w,t}=\frac{2\beta^{2}}{\mu t}$ and
$\eta_{\lambda,t}=\mathcal{O}(\frac{1}{t^{2/3}})$ and multiply both sides of
(29) by $t$ to get

$\displaystyle t(a_{t+1}+\nu_{t+1}b_{t+1})$
$\displaystyle\leq(t-1)(a_{t}+\nu_{t}b_{t})+\frac{2\kappa^{2}(L+\frac{L^{2}}{\rho})}{\beta^{2}t}$
$\displaystyle+\frac{E^{2}B_{w}^{2}(L+\frac{L^{2}}{\rho})}{\mu
t}+(\Gamma+B_{w}^{2}+B_{\lambda}^{2})\mathcal{O}(\frac{1}{t^{2/3}})$
$\displaystyle+\mathcal{O}(\frac{1}{t^{4/3}})$ (30)

Applying the inequality inductively from $t=1$ to $T$ and divide both sides by
$T$, we get

$\displaystyle\mathbb{E}[R_{S}(w_{T+1})-R^{*}_{S}]$
$\displaystyle\leq\frac{2\kappa^{2}(L+\frac{L^{2}}{\rho})}{\beta^{2}T^{2}}+\frac{E^{2}B_{w}^{2}(L+\frac{L^{2}}{\rho})}{\mu
T^{2}}$
$\displaystyle+(\Gamma+B_{w}^{2}+B_{\lambda}^{2})\mathcal{O}(\frac{1}{T^{2/3}})+\mathcal{O}(\frac{1}{T^{4/3}}).$
(31)

∎

##  Appendix B Additional Experiments on Different Levels of Heterogeneity

TABLE III: Comparison of the performance of the proposed algorithm across
baselines on CelebA and UTKFace datasets on other levels of heterogeneity
($\alpha=1.0$ and $\alpha=5.0$). $\uparrow$ indicates the larger the value the
better and $\downarrow$ indicates the smaller the value the better.

Algorithm | Acc $\uparrow(\%)$ |  $\Delta_{DP}$ $\downarrow(\%)$ |  $\Delta_{EO}$ $\downarrow(\%)$
---|---|---|---
CelebA ($\alpha=1.0$)
FedAvg | $76.76$ | $41.53$ | $22.66$
FairFed | $75.36$ | $39.37$ | $21.99$
FPFL | $77.54$ | $36.41$ | $16.94$
FFALM | $77.47$ | $35.00$ | $14.00$
CelebA ($\alpha=5.0$)
FedAvg | $77.36$ | $41.78$ | $22.79$
FairFed | $77.1$ | $41.40$ | $22.24$
FPFL | $79.28$ | $42.99$ | $22.48$
FFALM | $78.89$ | $41.55$ | $23.08$
UTKFace ($\alpha=1.0$)
FedAvg | $91.95$ | $15.59$ | $15.94$
FairFed | $89.14$ | $13.52$ | $13.23$
FPFL | $87.56$ | $11.08$ | $4.38$
FFALM | $90.39$ | $13.89$ | $3.87$
UTKFace ($\alpha=5.0$)
FedAvg | $92.02$ | $15.74$ | $9.35$
FairFed | $91.25$ | $15.5$ | $10.09$
FPFL | $86.74$ | $17.38$ | $7.37$
FFALM | $88.48$ | $11.78$ | $8.43$

## References

  * [1] 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 AISTATS (A. Singh and J. Zhu, eds.), vol. 54, pp. 1273–1282, PMLR, 2017.
  * [2] D. Pedreshi, S. Ruggieri, and F. Turini, “Discrimination-aware data mining,” in Proceedings of the 14th SIGKDD, ACM, 2008.
  * [3] European Parliament and Council of the European Union, “Regulation (EU) 2016/679 of the European Parliament and of the Council.”
  * [4] A. Howard and J. Borenstein, “The ugly truth about ourselves and our robot creations: The problem of bias and social inequity,” Science and Engineering Ethics, vol. 24, no. 5, pp. 1521–1536, 2017.
  * [5] P. Kairouz and et al., “Advances and open problems in federated learning,” Found. Trends Mach. Learn., vol. 14, no. 1–2, p. 1–210, 2021.
  * [6] S. Caton and C. Haas, “Fairness in machine learning: A survey,” ACM Comput. Surv., aug 2023.
  * [7] J. Nocedal and S. J. Wright, Numerical Optimization.  New York, NY, USA: Springer, 2e ed., 2006.
  * [8] T. Wang, J. Zhao, M. Yatskar, K.-W. Chang, and V. Ordonez, “Balanced datasets are not enough: Estimating and mitigating gender bias in deep image representations,” in ICCV, 2019.
  * [9] H. J. Ryu, H. Adam, and M. Mitchell, “Inclusivefacenet: Improving face attribute detection with race and gender diversity,” arXiv preprint arXiv:1405.6245, 2017.
  * [10] V. V. Ramaswamy, S. S. Y. Kim, and O. Russakovsky, “Fair attribute classification through latent space de-biasing,” in IEEE/CVF Conference on CVPR, 2021.
  * [11] V. S. Lokhande, A. K. Akash, S. N. Ravi, and V. Singh, “FairALM: Augmented lagrangian method for training fair models with little regret,” in Computer Vision – ECCV 2020, pp. 365–381, Springer International Publishing, 2020.
  * [12] Y. H. Ezzeldin, S. Yan, C. He, E. Ferrara, and A. S. Avestimehr, “FairFed: Enabling group fairness in federated learning,” Proceedings of the AAAI Conference on Artificial Intelligence, vol. 37, no. 6, pp. 7494–7502, 2023.
  * [13] S. Cui, W. Pan, J. Liang, C. Zhang, and F. Wang, “Addressing algorithmic disparity and performance inconsistency in federated learning,” in Advances in Neural Information Processing Systems (M. Ranzato, A. Beygelzimer, Y. Dauphin, P. Liang, and J. W. Vaughan, eds.), vol. 34, pp. 26091–26102, Curran Associates, Inc., 2021.
  * [14] B. R. Gálvez, F. Granqvist, R. van Dalen, and M. Seigel, “Enforcing fairness in private federated learning via the modified method of differential multipliers,” in NeurIPS 2021 Workshop Privacy in Machine Learning, 2021.
  * [15] Y. Zeng, H. Chen, and K. Lee, “Improving fairness via federated learning,” arXiv preprint arXiv:2110.15545, 2022.
  * [16] T. Calders and S. Verwer, “Three naive bayes approaches for discrimination-free classification,” Data Mining and Knowledge Discovery, vol. 21, no. 2, pp. 277–292, 2010.
  * [17] M. B. Zafar, I. Valera, M. Gomez Rodriguez, and K. P. Gummadi, “Fairness beyond disparate treatment & disparate impact: Learning classification without disparate mistreatment,” in Proceedings of the 26th International Conference on World Wide Web, WWW ’17, (Republic and Canton of Geneva, CHE), p. 1171–1180, International World Wide Web Conferences Steering Committee, 2017.
  * [18] M. Hardt, E. Price, and N. Srebro, “Equality of opportunity in supervised learning,” in Advances in Neural Information Processing Systems, vol. 29, Curran Associates, Inc., 2016.
  * [19] M. I. J. Peter L Bartlett and J. D. McAuliffe, “Convexity, classification, and risk bounds,” Journal of the American Statistical Association, vol. 101, no. 473, pp. 138–156, 2006.
  * [20] J. Wang and et al., “A field guide to federated optimization,” arXiv preprint arXiv:2107.06917, 2021.
  * [21] C. Liu, L. Zhu, and M. Belkin, “Toward a theory of optimization for over-parameterized systems of non-linear equations: the lessons of deep learning,” arXiv preprint arXiv:2003.00307, 2020.
  * [22] Z. Liu, P. Luo, X. Wang, and X. Tang, “Deep learning face attributes in the wild,” in Proceedings of ICCV, 2015.
  * [23] Z. Zhifei, S. Yang, and Q. Hairong, “Age progression/regression by conditional adversarial autoencoder,” in IEEE Conference on CVPR, IEEE, 2017.
  * [24] K. He, X. Zhang, S. Ren, and J. Sun, “Deep residual learning for image recognition,” arXiv preprint arXiv:1512.03385, 2015.
  * [25] H. Hsu, H. Qi, and M. Brown, “Measuring the effects of non-identical data distribution for federated visual classification,” arXiv preprint arXiv:1909.06335, 2019.
  * [26] Z. Yang, S. Hu, Y. Lei, K. R. Vashney, S. Lyu, and Y. Ying, “Differentially private sgda for minimax problems,” in Proceedings of the Thirty-Eighth Conference on Uncertainty in Artificial Intelligence (J. Cussens and K. Zhang, eds.), vol. 180 of Proceedings of Machine Learning Research, pp. 2192–2202, PMLR, 2022.
  * [27] T. Lin, C. Jin, and M. Jordan, “On gradient descent ascent for nonconvex-concave minimax problems,” in Proceedings of the 37th International Conference on Machine Learning (H. D. III and A. Singh, eds.), vol. 119 of Proceedings of Machine Learning Research, pp. 6083–6093, PMLR, 2020.
  * [28] H. Karimi, J. Nutini, and M. Schmidt, “Linear convergence of gradient and proximal-gradient methods under the polyak-łojasiewicz condition,” in Machine Learning and Knowledge Discovery in Databases - European Conference, ECML PKDD 2016, Riva del Garda, Italy, September 19-23, 2016, Proceedings, Part I (P. Frasconi, N. Landwehr, G. Manco, and J. Vreeken, eds.), vol. 9851 of Lecture Notes in Computer Science, pp. 795–811, Springer, 2016.


$\epsilon^{\prime}$. Observe that

$\displaystyle u\equiv\langle
0^{m+n_{\mathcal{G}}+1+n_{\mathcal{M}^{-1}}}|(U^{\dagger}_{initial}\otimes\mathbf{1}^{n_{\mathcal{G}}+1+n_{\mathcal{M}^{-1}}})(U^{\dagger}_{\mathcal{M}^{-1}}\otimes\mathbf{1}^{n_{\mathcal{G}}+1})U_{\mathcal{G}^{\prime}}(U_{\mathcal{M}^{-1}}\otimes\mathbf{1}^{n_{\mathcal{G}}+1})(U_{initial}\otimes\mathbf{1}^{n_{\mathcal{G}}+1+n_{\mathcal{M}^{-1}}})|0^{m+n_{\mathcal{G}}+1+n_{\mathcal{M}^{-1}}}\rangle$
$\displaystyle=\langle
0^{m+n_{\mathcal{M}^{-1}}}|(U^{\dagger}_{initial}\otimes\mathbf{1}^{n_{\mathcal{M}^{-1}}})U^{\dagger}_{\mathcal{M}^{-1}}(|0^{n_{\mathcal{M}^{-1}}}\rangle\langle
0^{n_{\mathcal{M}^{-1}}}|\otimes\mathcal{G})U_{\mathcal{M}^{-1}}(U_{initial}\otimes\mathbf{1}^{n_{\mathcal{M}^{-1}}})|0^{m+n_{\mathcal{M}^{-1}}}\rangle$
$\displaystyle=\langle\psi_{0}|(0^{n_{\mathcal{M}^{-1}}}|U^{\dagger}_{\mathcal{M}^{-1}}|0^{n_{\mathcal{M}^{-1}}}\rangle)\mathcal{G}(0^{n_{\mathcal{M}^{-1}}}|U_{\mathcal{M}^{-1}}|0^{n_{\mathcal{M}^{-1}}}\rangle)|\psi_{0}\rangle=\langle\psi_{0}|\mathcal{L}^{\dagger}\mathcal{G}\mathcal{L}|\psi_{0}\rangle,$
(137)

where $\mathcal{L}\equiv\langle
0^{n_{\mathcal{M}^{-1}}}|U_{\mathcal{M}^{-1}}|0^{n_{\mathcal{M}^{-1}}}\rangle$
and $|\psi_{0}\rangle=U_{initial}|0^{m}\rangle$. Since block-access to
$\mathcal{M}^{-1}$ can be made with small enough error $\delta$ as will be
seen later, then

$\displaystyle\alpha^{2}_{\mathcal{M}^{-1}}u\approx\langle\psi_{0}|(\mathcal{M}^{-1})^{\dagger}\mathcal{G}\mathcal{M}^{-1}|\psi_{0}\rangle\equiv\Upsilon.$
(138)

The benefit of approximating $\Upsilon$ using $u$ is that $u$ is written in
terms of the expectation value with respect to unitary operators, so one is
able to more elegantly employ the amplitude estimation algorithm [41] to gain
a quadratic speedup with respect to error in computing $\Upsilon$. This allows
an optimal estimation of the expectation value.

The amplitude estimation algorithm [41] is an algorithm for estimating the
value $u\equiv|\langle 0^{r}|U^{\dagger}VU|0^{r}\rangle|_{u}$ to
$\epsilon_{u}$-additive precision with success probability at least $2/3$,
where $U$, $V$ are $r$-qubit unitary black-boxes. Here we have

$\displaystyle
U\equiv(U_{\mathcal{M}^{-1}}\otimes\mathbf{1}^{n_{\mathcal{G}}+1})(U_{initial}\otimes\mathbf{1}^{n_{\mathcal{G}}+1+n_{\mathcal{M}^{-1}}}),$
$\displaystyle V\equiv U_{\mathcal{G}^{\prime}},$ $\displaystyle r\equiv
m+n_{\mathcal{G}}+1+n_{\mathcal{M}^{-1}}.$ (139)

The key is to construct a unitary operation $S$ whose eigenvalue is
$\exp(i\theta)$ where $u=|\cos(\theta/2)|$. Then a quantum phase estimation
algorithm is used to extract $u$. Defining $|\phi_{0}\rangle\equiv
U|0^{r}\rangle$ and $|\phi_{1}\rangle\equiv VU|0^{r}\rangle$, our desired
quantity is then the inner product $u=\langle\phi_{0}|\phi_{1}\rangle$. The
aim is then to construct $S$ as a rotation operator with eigenvalue
$\exp(i\theta)$, that rotates $|\phi_{0}\rangle$ to $|\phi_{1}\rangle$ and the
two states are separated by angle
$2\theta=4\cos^{-1}(\langle\phi_{0}|\phi_{1}\rangle)$. Just like in Grover’s
search algorithm, one can construct the rotation operator as a combination of
two reflection operators $S=S_{0}S_{1}$ where
$S_{0}=\mathbf{1}^{r}-2|\phi_{0}\rangle\langle\phi_{0}|=UP_{0}U^{\dagger}$ and
$S_{1}=\mathbf{1}^{r}-2|\phi_{1}\rangle\langle\phi_{1}|=VUP_{0}U^{\dagger}V^{\dagger}$
where $P_{0}=\mathbf{1}^{r}-2|\phi_{0}\rangle\langle\phi_{0}|$. This means
that one can construct the unitary operator

$\displaystyle S=UP_{0}U^{\dagger}VUP_{0}U^{\dagger}V^{\dagger}$ (140)

by concatenating the unitary black-boxes $U_{\mathcal{G}}$,
$U_{\mathcal{M}^{-1}}$ and $U_{initial}$ and their adjoints, which we already
assume one has access to. Then it is straightforward to apply the standard
quantum phase estimation algorithm using controlled-$S$, i.e.,
$|0\rangle\langle 0|\otimes\mathbf{1}+|1\rangle\langle 1|\otimes S$, which can
be easily constructed from controlled-$U$, controlled-$V$ and are in turn
constructed from controlled-$U_{\mathcal{M}^{-1}}$ and
controlled-$U_{\mathcal{G}}$. The output of the quantum phase estimation
algorithm is then a $2\epsilon_{u}$-additive error estimate of $\theta$, which
leads to an $\epsilon_{u}$-additive error estimate of $u=|\cos(\theta/2)|$.

Note that $u=|u|$ for our problem since it has the interpretation of being
proportional to quantum fidelity, so we don’t need to be concerned about
possible negative expectation values $\langle 0|U^{\dagger}VU|0\rangle$ that
can happen for general $U,V$. This means one also does not require the extra
steps in [36] to deal with possible negative expectation values.

To obtain an $\epsilon_{u}$-additive error estimate to $u$, one can use the
following lemma.

######  Lemma 35.

([36, 41]) An amplitude estimation algorithm exists that gives a
$\epsilon_{u}$-additive estimate of $u$ that makes
$\mathcal{O}(1/\epsilon_{u})$ queries to $U$, $V$ and
$\mathcal{O}(r/\epsilon_{u})$ additional $2$-qubit gates.

###### Proof.

See [41] and Lemma II.15 from [36]. ∎

Since $\Upsilon\approx\alpha^{2}_{\mathcal{M}^{-1}}u$, one also needs to
identify $\alpha_{\mathcal{M}^{-1}}$ to find the total query complexity. To
demonstrate an appropriate value, one first requires three other lemmas.

######  Lemma 36.

([54, 36]) Given block-access
$(\alpha_{\mathcal{M}},n_{\mathcal{M}},\delta_{\mathcal{M}},U_{\mathcal{M}})$
to matrix $\mathcal{M}$, and error $\sigma>0$ and a polynomial
$\mathcal{P}:\mathbb{R}\rightarrow\mathbb{R}$ of degree $D(\mathcal{P})$
satisfying $|\mathcal{P}(x)|\leq 1/2$ for all $x\in[-1,1]$, then block access
$(1,n_{\mathcal{M}},4D(\mathcal{P})\sqrt{\delta_{\mathcal{M}}/\alpha_{\mathcal{M}}}+\sigma,U_{\mathcal{P}(\mathcal{M}/\alpha_{\mathcal{M}})})$
to $\mathcal{P}(\mathcal{M}/\alpha_{\mathcal{M}})$ can be constructed where
each $U_{\mathcal{P}(\mathcal{M}/\alpha_{\mathcal{M}})}$ makes
$2D(\mathcal{P})+1$ queries to $U_{\mathcal{M}}$ and
$\mathcal{O}(n_{\mathcal{M}}D(\mathcal{P}))$ additional $2$-qubit gates.

###### Proof.

See Theorem 56 [54] and Theorem II.6 from [36]. ∎

######  Lemma 37.

([54]) For any $\chi,\zeta\in(0,1/2]$, a polynomial
$\mathcal{P}:\mathbb{R}\rightarrow\mathbb{R}$ with odd degree
$\mathcal{O}(\log(1/\zeta)/\chi)$ exists such that for all
$x\in[-1,1]\backslash[-\chi,\chi]$, $|\mathcal{P}(x)|\leq 1$ and
$|\mathcal{P}(x)-\chi/(2x)|<\zeta$.

###### Proof.

See Corollary 67 in [54]. ∎

Now we can prove the following lemma and show how the choice
$\alpha_{\mathcal{M}^{-1}}=4\theta\kappa/\|\mathcal{M}\|$ for any $\theta>1$
is possible.

######  Lemma 38.

([43]) If given $m$ qubits and block access
$(\alpha_{\mathcal{M}},n_{\mathcal{M}},\delta_{\mathcal{M}},U_{\mathcal{M}})$
to $\mathcal{M}$ with $\|\mathcal{M}\|\leq\alpha_{\mathcal{M}}$, then block
access
$(\alpha_{\mathcal{M}^{-1}}=4\kappa\theta/\|\mathcal{M}\|,n_{\mathcal{M}^{-1}}=n_{\mathcal{M}}+2,\delta,U_{\mathcal{M}^{-1}})$
to $\mathcal{M}^{-1}$ can be constructed, for any $\theta>1$, such that
$U_{\mathcal{M}^{-1}}$ makes
$\mathcal{O}((\alpha_{\mathcal{M}}\kappa\theta/\|\mathcal{M}\|)\log(\theta\kappa/(\|\mathcal{M}\|\delta)))$
queries to $U_{\mathcal{M}}$ and $\mathcal{O}((\alpha_{\mathcal{M}}\kappa
n_{\mathcal{M}}\theta/\|\mathcal{M}\|)\log(\kappa\theta/(\|\mathcal{M}\|\delta)))$
additional $2$-qubit gates, with
$\delta\geq(2\theta\kappa\zeta+4\theta\kappa\sigma+16\theta^{2}\kappa^{2}v\log(1/\zeta)\sqrt{\delta_{\mathcal{M}}\alpha_{\mathcal{M}}})/\|\mathcal{M}\|\
$, where appropriate $\sigma>0$, $\zeta\in[0,1/2)$ can be chosen, $v>0$ and
$\delta\in(0,2]$.

###### Proof.

This proof is an extension of Corollary IV. 15 in [36] for the case of
$\delta_{\mathcal{M}}>0$ and similar steps are used. Here in order to use
Lemma 36 where $x\in[-1,1]\backslash[-\chi,\chi]$ and we want $x$ to represent
the spectrum of the operator $\mathcal{M}/\alpha_{\mathcal{M}}$, its spectrum
should also lie within $[-1,1]\backslash[-\chi,\chi]$. One can just look at
the positive eigenvalues without losing generality. If $\lambda_{min}$ and
$\lambda_{max}$ represent the smallest and largest absolute values of the
eigenvalues of $\mathcal{M}$, then the above requirement demands

$\displaystyle\chi<\frac{\lambda_{min}}{\alpha_{\mathcal{M}}}\leq\text{Spec}\left(\frac{\mathcal{M}}{\alpha_{\mathcal{M}}}\right)\leq\frac{\lambda_{max}}{\alpha_{\mathcal{M}}}\leq
1\,.$ (141)

Since by definition $\|\mathcal{M}\|=\lambda_{max}$ and
$\|\mathcal{M}^{-1}\|=1/\lambda_{min}$, the last inequality gives the
condition $\|\mathcal{M}\|\leq\alpha_{\mathcal{M}}$, which we note to be
sufficient for the inequality
$\|\mathcal{M}\|\leq\alpha_{\mathcal{M}}+\delta_{\mathcal{M}}$ demanded by the
definition of the block access to $\mathcal{M}$. Then the first equality in
combination with the definition $\|\mathcal{M}\|\|\mathcal{M}^{-1}\|=\kappa$
gives

$\displaystyle\chi<\frac{\|\mathcal{M}\|}{\alpha_{\mathcal{M}}\kappa}.$ (142)

When $\|\mathcal{M}\|\leq\alpha_{\mathcal{M}}$, $\chi<1/\kappa$ suffices. (In
our application to PDEs and ODEs, $\|\mathcal{M}\|<2$). One can choose for
instance $\chi=\|\mathcal{M}\|/(\theta\alpha_{\mathcal{M}}\kappa)$ for any
constant $\theta>1$. Then interpreting $x$ in Lemma 37 to be the spectrum of
$\mathcal{M}/\alpha_{\mathcal{M}}$ one can replace the inequality

$\displaystyle|\chi/4x-\mathcal{P}(x)/2|\leq\zeta/2$ (143)

by

$\displaystyle\|\mathcal{M}^{-1}-(4\theta\kappa/\|\mathcal{M}\|)\mathcal{P}(\mathcal{M}/\alpha_{\mathcal{M}})/2\|\leq
2\theta\kappa\zeta/\|\mathcal{M}\|$ (144)

By defining the block encoding
$U_{\mathcal{M}^{-1}}=U_{\mathcal{P}(\mathcal{M}/\alpha_{\mathcal{M}})/2}$
then the definition of $\delta$ requires

$\displaystyle\|\mathcal{M}^{-1}-(4\theta\kappa/\|\mathcal{M}\|)\langle
0^{n_{\mathcal{M}}+2}|U_{\mathcal{P}(\mathcal{M}/\alpha_{\mathcal{M}})/2}|0^{n_{\mathcal{M}}+2}\rangle\|\leq\delta.$
(145)

Lemmas 36 and 37 imply

$\displaystyle\|\langle
0^{n_{\mathcal{M}}+2}|U_{\mathcal{P}(\mathcal{M}/\alpha_{\mathcal{M}})/2}|0^{n_{\mathcal{M}}+2}\rangle-\mathcal{P}(\mathcal{M}/\alpha_{\mathcal{M}})/2\|\leq
4D(\mathcal{P})\sqrt{\delta_{\mathcal{M}}/\alpha_{\mathcal{M}}}+\sigma$ (146)

where it is sufficient to choose
$D(\mathcal{P})=v\log(1/\zeta)/\chi=v\theta\alpha_{\mathcal{M}}\kappa\log(1/\zeta)/\|\mathcal{M}\|$
for some constant $v>0$. Putting together Eqs. (144), (145) and (146) gives

$\displaystyle\|\mathcal{M}^{-1}-(4\theta\kappa/\|\mathcal{M}\|)\langle
0^{n_{\mathcal{M}}+2}|U_{\mathcal{P}(\mathcal{M}/\alpha_{\mathcal{M}})/2}|0^{n_{\mathcal{M}}+2}\rangle\|\leq\frac{1}{\|\mathcal{M}\|}(2\theta\kappa\zeta+4\theta\kappa\sigma+16\theta^{2}\kappa^{2}v\log(1/\zeta)\sqrt{\delta_{\mathcal{M}}\alpha_{\mathcal{M}}})\leq\delta.$
(147)

For example, in the limit $\delta_{\mathcal{M}}=0$, one can choose the
parameters $\zeta=\delta\|\mathcal{M}\|/(8\kappa)$ and
$\sigma=\delta\|\mathcal{M}\|/(16\kappa)$. However, we are interested in the
case $\delta_{\mathcal{M}}>0$. Then it is sufficient to choose
$\zeta=\delta\|\mathcal{M}\|(1-\theta^{\prime})/(4\theta\kappa)$,
$\sigma=\delta\|\mathcal{M}\|(1-\theta^{\prime})/(8\theta\kappa)$ and
$\sqrt{\delta_{\mathcal{M}}\alpha_{\mathcal{M}}}\leq\delta\|\mathcal{M}\|\theta^{\prime}/(16\theta^{2}v\kappa^{2}\log(8\kappa/(\|\mathcal{M}\|\delta)))$
for any constant $0<\theta^{\prime}<1$. The latter implies
$\delta_{\mathcal{M}}<\delta^{2}\|\mathcal{M}\|^{2}(\theta^{\prime})^{2}/(16^{2}\alpha_{\mathcal{M}}\kappa^{4}\log^{2}(\kappa/\delta))$.
From Lemma 32 and ignoring constant factors
$\theta,\theta^{\prime},v,\|\mathcal{M}\|=O(1)$ this gives rise to an
additional
$O(\log^{2.5}(\alpha_{\mathcal{M}}^{2}\kappa^{4}\log^{2}(\kappa/\delta)/\delta^{2}))$
2-qubit gates to create block access to $\mathcal{M}$ from sparse access.
Later from Eq. (F) we see that $\delta\sim\epsilon^{\prime}/\kappa^{2}$ can be
chosen where $\epsilon^{\prime}$ is the final error in $\Upsilon$, so this
gives the additional gate cost
$O(\log^{2.5}(s^{2}\|\mathcal{M}\|_{max}^{2}\kappa^{8}\log^{2}(\kappa^{3}/\epsilon^{\prime})/(\epsilon^{\prime})^{2})<O(\log^{2.5}(s\|\mathcal{M}\|_{max}\kappa^{4}/\epsilon^{\prime}))$.

Since one can choose
$\mathcal{D}(\mathcal{P})=v\theta\kappa\alpha_{\mathcal{M}}\log(4\theta\kappa/(\|\mathcal{M}\|\delta(1-\theta^{\prime})))/\|\mathcal{M}\|$
and then ignoring all constants except $\|\mathcal{M}\|$ and $\theta$, the
rest of the proof follows using Lemma 36. ∎

From the above lemma one sees one can set
$\alpha_{\mathcal{M}^{-1}}=4\theta\kappa/\|\mathcal{M}\|$ for any $\theta>1$
and $n_{\mathcal{M}^{-1}}=n_{\mathcal{M}}+2$. Let the total error in
estimating $\Upsilon$ be $\epsilon^{\prime}$. There are two sources of error:
one in the amplitude estimation algorithm that outputs $\tilde{u}$, which
approximates $u$ with error $\epsilon_{u}$ and the other error is in the
block-encoding of $\mathcal{M}^{-1}$, which has error $\delta$. This means

$\displaystyle|u-\tilde{u}|\leq\epsilon_{u}$ (148)

where $\tilde{u}$ is the estimate of $u$ from the amplitude estimation
algorithm and

$\displaystyle\|\alpha_{\mathcal{M}^{-1}}\mathcal{L}-\mathcal{M}^{-1}\|\leq\delta.$
(149)

Then the total error

$\displaystyle|\alpha^{2}_{\mathcal{M}^{-1}}u-\Upsilon|=|\alpha^{2}_{\mathcal{M}^{-1}}\langle\psi_{0}|\mathcal{L}^{\dagger}\mathcal{G}\mathcal{L}|\psi_{0}\rangle-\langle\psi_{0}|(\mathcal{M}^{-1})^{\dagger}\mathcal{G}\mathcal{M}^{-1}|\psi_{0}\rangle|$
$\displaystyle=|\alpha^{2}_{\mathcal{M}^{-1}}\|\mathcal{F}\mathcal{L}|\psi_{0}\rangle\|^{2}-\|\mathcal{F}\mathcal{M}^{-1}|\psi_{0}\rangle\|^{2}|$
$\displaystyle\leq\|\alpha^{2}_{\mathcal{M}^{-1}}\mathcal{F}\mathcal{L}|\psi_{0}\rangle-\mathcal{F}\mathcal{M}^{-1}|\psi_{0}\rangle\|^{2}$
$\displaystyle\leq\|\mathcal{F}\|^{2}\|\alpha^{2}_{\mathcal{M}^{-1}}\mathcal{L}|\psi_{0}\rangle-\mathcal{M}^{-1}|\psi_{0}\rangle\|^{2}$
$\displaystyle\leq\|\mathcal{F}\|^{2}\delta^{2}=\delta^{2}$ (150)

where in the second line we used
$\mathcal{G}=\mathcal{F}^{\dagger}\mathcal{F}$ since $\mathcal{G}$ is positive
semi-definite. In the last line we used $\|\mathcal{F}\|=\|\mathcal{G}\|=1$.
This means that the total error in $\Upsilon$ can be written as

$\displaystyle|\alpha^{2}_{\mathcal{M}^{-1}}\tilde{u}-\Upsilon|\leq|\alpha^{2}_{\mathcal{M}^{-1}}\tilde{u}-\alpha^{2}_{\mathcal{M}^{-1}}u|+|\alpha^{2}_{\mathcal{M}^{-1}}u-\Upsilon|\leq\alpha^{2}_{\mathcal{M}^{-1}}\epsilon_{u}+\delta^{2}\leq\epsilon^{\prime}$
(151)

Then a choice of
$\alpha^{2}_{\mathcal{M}^{-1}}\epsilon_{u}\sim\delta^{2}\sim\epsilon^{\prime}/2$
is sufficient. Since $\alpha_{\mathcal{M}^{-1}}=4\theta\kappa/\|\mathcal{M}\|$
one can choose

$\displaystyle\delta=\sqrt{\epsilon^{\prime}/2}$
$\displaystyle\epsilon_{u}=\epsilon^{\prime}\|\mathcal{M}\|/(32\theta^{2}\kappa^{2}).$
(152)

For instance, it is possible to set $\theta^{2}=33/32$ and let
$\epsilon_{u}=\epsilon^{\prime}\|\mathcal{M}\|/(33\kappa^{2})$. Using Lemma
38, inserting Eq. (F) and ignoring all constants except for
$\|\mathcal{M}\|=O(1)$ for convenience, it is sufficient for the amplitude
estimation algorithm to make
$\mathcal{O}(1/\epsilon_{u})\sim\mathcal{O}(\kappa^{2}/(\|\mathcal{M}\|\epsilon^{\prime}))$
queries to $U_{\mathcal{G}^{\prime}}$ and $U_{initial}$. To query
$U_{\mathcal{M}}$, one must multiply the amplitude estimation cost with the
query cost for $\mathcal{M}\rightarrow\mathcal{M}^{-1}$, hence a total
$\mathcal{O}(\alpha_{\mathcal{M}}\kappa\log(\kappa/\delta)/\epsilon_{u})\sim\mathcal{O}(\alpha_{\mathcal{M}}\kappa^{3}\log(\kappa^{2}/(\epsilon^{\prime})/(\|\mathcal{M}\|\epsilon^{\prime}))$
queries to $U_{\mathcal{M}}$, and
$\mathcal{O}(r/\epsilon_{u})\sim\mathcal{O}(\kappa^{2}(m+n_{\mathcal{G}}+1+n_{\mathcal{M}^{-1}})/(\|\mathcal{M}\|\epsilon^{\prime}))$
additional $2$-qubit gates.

Since we are given the block encoding to $\mathcal{G}$ rather than
$\mathcal{G}^{\prime}$, Lemma 34 requires only one query to $U_{\mathcal{G}}$
and $\mathcal{O}(n^{2}_{\mathcal{M}})\sim\mathcal{O}(1)$ additional $2$-qubit
gates, since $n_{\mathcal{M}}=2$ from Lemma 32.

From Lemmas 36 and 38, to construct block access to $\mathcal{M}^{-1}$ from
block access to $\mathcal{M}$ requires an additional
$\mathcal{O}(n_{\mathcal{M}}\kappa\alpha_{\mathcal{M}}\log(\kappa^{3}/\epsilon^{\prime}))$
$2$-qubit gates. Although from the proof of Lemma 38 one sees that
$O(\log^{2.5}(s\|\mathcal{M}\|_{max}\kappa^{4}/\epsilon^{\prime}))$ additional
2-qubit gates are required to construct block access from sparse access to
$\mathcal{M}$, these are all logarithmic factors, which we will ignore in the
final expression.

Putting all these results together with
$\alpha_{\mathcal{M}}=s\|\mathcal{M}\|_{max}$,
$n_{\mathcal{M}^{-1}}=n_{\mathcal{M}}+2$ and ignoring all constants except
$\|\mathcal{M}\|$, one finds that to approximate $\Upsilon$ to precision
$\epsilon^{\prime}$, one needs
$\mathcal{O}(\kappa^{2}/(\|\mathcal{M}\|\epsilon^{\prime}))$ queries to
$U_{\mathcal{G}}$ and $U_{initial}$,
$\mathcal{O}(s\|\mathcal{M}\|_{max}\kappa^{3}\log(\kappa^{2}/\epsilon^{\prime})/(\|\mathcal{M}\|\epsilon^{\prime}))$
queries to sparse oracles for $\mathcal{M}$, and
$\mathcal{O}((\kappa^{2}/(\|\mathcal{M}\|\epsilon^{\prime}))(m+n_{\mathcal{G}}+1+n_{\mathcal{M}}+2+s\|\mathcal{M}\|_{max}n_{\mathcal{M}}\kappa\log(\kappa^{2}/(\|\mathcal{M}\|\epsilon^{\prime}))))$
additional $2$-qubit gates. Since $n_{\mathcal{M}}=2$, $n_{\mathcal{G}}=2m+1$
for our scenario, this leads to
$\mathcal{O}((\kappa^{2}/(\|\mathcal{M}\|\epsilon^{\prime}))(3m+6+2s\|\mathcal{M}\|_{max}\kappa\log(\kappa^{2}/(\|\mathcal{M}\|\epsilon^{\prime}))))\sim\mathcal{O}((\kappa^{2}/(\|\mathcal{M}\|\epsilon^{\prime}))(m+s\|\mathcal{M}\|_{max}\kappa\log(\kappa^{2}/(\|\mathcal{M}\|\epsilon^{\prime}))))$
where we ignore the constant factor terms.

##  Appendix G Proof of Theorem 17

To estimate $\langle G(t_{n},x)\rangle$ to precision $\epsilon$, Lemma 15
requires the estimation of $\sqrt{\Upsilon}$ to precision
$\epsilon_{G}\sim\epsilon/(n_{G}n_{\psi_{0}})$. Using the relation
$\Delta(x^{2})=2x\Delta x$, where $\Delta x$ is the error in $x$ and
$\Delta(x^{2})$ is the error in $x^{2}$, taking $x=\sqrt{\Upsilon}$, one finds
that the additive error $\epsilon^{\prime}$ in $\Upsilon$ has size
$\epsilon^{\prime}\sim\sqrt{\Upsilon}\epsilon_{G}$. To identify how
$\sqrt{\Upsilon}$ would scale with $N$, from Eq. (54) one observes that
$O(1)=G^{n}_{\omega}(\boldsymbol{j})$ implies $\sqrt{\Upsilon}\sim
1/(n_{G}n_{\psi_{0}})$ which gives
$\epsilon^{\prime}\sim\epsilon_{G}/(n_{G}n_{\psi_{0}})\sim\epsilon/(n_{G}n_{\psi_{0}})^{2}\sim\epsilon/n_{\psi_{0}}^{2}$
when we suppress the $n_{G}=O(1)$ factor.

Then following the quantum algorithm outlined in Lemma 16 with
$2^{m}=N_{t}N^{2d}$, we include $s=O(d)$ and $\kappa\leq O(N_{t})$ from Lemma
10. Since $N_{t}=T/\Delta t$, where from the stability condition we have
$\Delta t\sim 1/(Nd)$, this gives $\kappa\leq O(TNd)$. Inserting this and
$\epsilon^{\prime}\sim\epsilon/n_{\psi_{0}}^{2}$ into Lemma 16, one directly
finds that one needs to make
$\mathcal{O}(T^{2}d^{2}N^{2}n^{2}_{\psi_{0}}/\epsilon)$ queries to
$L(n,\boldsymbol{j})$ and $U_{initial}$ and
$\mathcal{O}((n_{\psi_{0}}^{2}/\epsilon)d^{4}N^{3}T^{3}\log(d^{2}N^{2}T^{2}n^{2}_{\psi_{0}}/\epsilon))$
queries to sparse oracles for $\mathcal{M}$. One also needs an additional
$\mathcal{O}\left((n_{\psi_{0}}^{2}d^{2}N^{2}T^{2}/\epsilon)(\log(TNd)+d\log(N)+d^{2}NT\log(T^{2}d^{2}N^{2}n_{\psi_{0}}^{2}/\epsilon)\right)$
two-qubit gates. The largest of these terms is
$\mathcal{O}\left((n^{2}_{\psi_{0}}d^{4}N^{3}T^{3}/\epsilon)\log(d^{2}N^{2}T^{2}n^{2}_{\psi_{0}}/\epsilon)\right)$.

To express this entirely in terms of the natural parameters $\epsilon$, $T$,
$M$ and $d$ only, we note that in the quantum algorithm we require
$\epsilon\sim\epsilon_{\text{CL}}$ from Lemma 15 and here $N=N_{\text{CL}}$,
so from Lemma 13 we have $N\sim d/\epsilon^{3}$. This easily gives us
$\mathcal{Q}=\mathcal{O}(n^{2}_{\psi_{0}}T^{3}d^{7}(1/\epsilon)^{10}\log(n^{2}_{\psi_{0}}T^{2}d^{4}(1/\epsilon)^{7})$.

The above analysis is only for computing $\Upsilon$ to the required accuracy
to estimate $\langle G(t_{n},x)\rangle$. However, to obtain our observable it
is insufficient to compute $\Upsilon$ alone since the normalisation constants
need to be considered to derive $\langle
G^{n}_{\omega}(\boldsymbol{j})\rangle$ from Eq. (54). The normalisation
constant $n_{\mathcal{G}}$ can be computed easily on a classical device since
the function $G(p)$ is known and is generally a simple polynomial. Since the
initial data is also known, it is also sufficient to use a classical device to
first compute $n_{\psi_{0}}$. Alternatively, if one doesn’t wish to compute
$n_{\psi_{0}}$, one can instead estimate the observable in Definition 11,
where $G_{O}(t_{n},x)\equiv\langle G(t_{n},x)\rangle/\langle
G^{\prime}(t_{n},x)\rangle$ with $G^{\prime}(p)\equiv 1$, where
$n_{G^{\prime}}=1$ exactly. This means we estimate

$\displaystyle G_{O}(t_{n},x)\equiv\frac{\langle
G^{n}_{\omega}(\boldsymbol{j})\rangle}{\mathbf{1}^{n}_{\omega}(\boldsymbol{j})}=n_{G}\frac{|\sqrt{\Upsilon}|}{|\sqrt{\Upsilon_{+}}|}$
(153)

where
$\Upsilon_{+}\equiv\sqrt{\langle\psi_{0}|(\mathcal{M}^{-1})^{\dagger}\mathcal{H}\mathcal{M}^{-1}|\psi_{0}\rangle}$
and $\mathcal{H}=|H\rangle\langle H|$ and $|H\rangle$ is the state with equal
superposition across the $\boldsymbol{l}$ basis and can be easily created by
applying Hadamard gates on $|0\rangle$. Since $G_{O}(t_{n},x)=O(1)$ just like
$\langle G(t_{n},x)\rangle=O(1)$, the error
$|G_{O}(t_{n},x)-\tilde{G}_{O}(t_{n},x)|$ is also $O(\epsilon)$.

##  Appendix H Basic summary of quantum algorithm to estimate $\Upsilon$

Here we present a condensed version of the main steps of the quantum algorithm
in Theorems 17 and 30 (where we use $\mathcal{M}=\mathcal{M}_{\text{ODE}}$).
For query and gate complexity at each step refer to the lemmas referenced:

  1. 1.

(Step 0) Inputs: sparse access
$(s,\|\mathcal{M}\|_{max},O_{\mathcal{M}},O_{F})$ to $m$-qubit operator
$\mathcal{M}$, access to $U_{initial}$ and access to unitary
$L(n,\boldsymbol{j})$;

  2. 2.

(Step 1) Using sparse access to $\mathcal{M}$ to construct block access
$(s\|\mathcal{M}\|_{max},2,\delta_{\mathcal{M}},U_{\mathcal{M}})$ to
$\mathcal{M}$, using Lemma 32, so
$\alpha_{\mathcal{M}}=s\|\mathcal{M}\|_{max}$ and $n_{\mathcal{M}}=2$;

  3. 3.

(Step 2) Using block access to
$(s\|\mathcal{M}\|_{max},2,\delta_{\mathcal{M}},U_{\mathcal{M}})$ to
$\mathcal{M}$ to construct block access
$(4\theta\kappa/\|\mathcal{M}\|,4,\delta,U_{\mathcal{M}^{-1}})$ to
$\mathcal{M}^{-1}$ where
$\alpha_{\mathcal{M}^{-1}}=4\theta\kappa/\|\mathcal{M}\|$,
$n_{\mathcal{M}^{-1}}=n_{\mathcal{M}}+2=4$, from Lemma 38;

  4. 4.

(Step 3) Starting from access to unitary $L(n,\boldsymbol{j})$ can create
block access $(1,n_{\mathcal{G}},0,U_{\mathcal{G}})$ to $\mathcal{G}$ with
construction given in Lemma 33 with choice $n_{\mathcal{G}}=2m+1$;

  5. 5.

(Step 4) Using block access $(1,n_{\mathcal{G}},0,U_{\mathcal{G}})$ to
$\mathcal{G}$ to construct block access $(1,2m+2,0,U_{\mathcal{G}^{\prime}})$
to $\mathcal{G}^{\prime}=|0^{4}\rangle\langle 0^{4}|\otimes\mathcal{G}$, where
$\alpha_{\mathcal{G}^{\prime}}=1$,
$n_{\mathcal{G}^{\prime}}=n_{\mathcal{G}}+1=2m+2$, from Lemma 34;

  6. 6.

(Step 5) Given $U_{\mathcal{M}^{-1}}$ and $U_{\mathcal{G}^{\prime}}$, and also
assume access to their adjoints and their respective controlled- unitaries.
From this to construct controlled-$S$ and the corresponding quantum phase
estimation circuit, where $S$ is in Eqs. (F) and (140). Extract $u$;

  7. 7.

Output: Multiply $u$ by $\alpha^{2}_{\mathcal{M}^{-1}}$ where
$\alpha_{\mathcal{M}^{-1}}=4\kappa\theta/\|\mathcal{M}\|$. This estimates
$\Upsilon$ to precision $\epsilon^{\prime}$.

Figure 2: Summary of algorithm to estimate $\Upsilon$ to error
$\epsilon^{\prime}$

##  Appendix I Proof of Lemma 23

Let $\chi(q)$ be a smooth test function, and consider

$\int_{\mathbb{R}^{d}}\chi(q)\Phi(t,q)dq=\frac{1}{M}\sum_{k=1}^{M}\chi(X^{[k]}(t)).$

Taking the time derivative of this equation and using (12), one deduces

$\displaystyle\quad\int_{\mathbb{R}^{d}}\chi(q)\partial_{t}\Phi(t,q)dq=\frac{1}{M}\sum_{k=1}^{M}\partial_{t}\chi(X^{[k]})=\frac{1}{M}\sum_{k=1}^{M}\nabla_{q}\chi(X^{[k]})\cdot\partial_{t}X^{[k]}(t)=\frac{1}{M}\sum_{k=1}^{M}\nabla_{q}\chi(X^{[k]})\cdot
F(X^{[k]})$
$\displaystyle=\frac{1}{M}\int_{\mathbb{R}^{d}}\nabla_{q}\chi(q)\cdot
F(q)\,\sum_{k=1}^{M}\delta(q-X^{[k]}(t))\,dq=\int_{\mathbb{R}^{d}}\nabla_{q}\chi(q)\cdot
F(q)\Phi(t,q)\,dq=-\int_{\mathbb{R}^{d}}\chi(q)\nabla_{q}\cdot[F(q)\Phi(t,q)]\,dq.$
(154)

The uniqueness of the solution is also classical, see [48]. We omit the
details.

##  Appendix J Proof of Lemma 26

Consider the problem

$\displaystyle\frac{\partial X(t)}{\partial t}=F(X),\qquad
X\in\mathbb{R}^{D}\,,$ (155) $\displaystyle X(s)=x,\quad s\in[0,T].$ (156)

Let the solution to (155) be $X(t;x,s)$. Define the Jacobian determination of
the map from $x$ to $X$ be

$J(t;x,s)={\text{det}}\left(\frac{\partial X_{i}}{\partial
x_{j}}(t;x,s)\right)$

then classical result (see for example [48]) shows that

$J(t;x,s)>0,\qquad J(t;x,s)=\exp\int_{s}^{t}\nabla\cdot
F(X(\sigma;x,s))d\sigma$

and, by the method of characteristics,

$\Phi(t,p)=\Phi(0,X(0;p,t))J(0;p,t)$

Now,

$\displaystyle A_{p}(t)$
$\displaystyle=\int_{\mathbb{R}^{D}}A(p)\Phi(t,p)\,dp=\int_{\mathbb{R}^{D}}A(p)\Phi(0,X(0;p,t))J(0;p,t)\,dp$
$\displaystyle=\int_{\mathbb{R}^{D}}A(p)\Phi(0,X(0;p,t))J(0;p,t)\,dp=\frac{1}{M}\sum_{k=1}^{M}\int_{\mathbb{R}^{D}}A(p)\delta(X(0;p,t)-X_{0}^{[k]})J(0;p,t)\,dp$
$\displaystyle=\frac{1}{M}\sum_{k=1}^{M}\int_{\mathbb{R}^{D}}A(p)J(0;p,t)^{-1}\delta(p-X(t;X_{0}^{[k]},0))J(0;p,t)\,dp$
$\displaystyle=\frac{1}{M}\sum_{k=1}^{M}\int_{\mathbb{R}^{D}}A(p)\delta(p-X(t;X_{0}^{[k]},0))\,dp=\frac{1}{M}\sum_{k=1}^{M}A(X(t;X_{0}^{[k]},0))=A_{0}(t)\,.$

Note here, by the definition in (155), $X(t;X_{0}^{[k]},0)$ is the solution to
(12), hence the last equality holds.

##  Appendix K Discretised System of ODEs

We now discretise the linear PDE in Eq. (V.1) by finite difference schemes. As
an example we use the upwind scheme, which takes the following form:

$\displaystyle\frac{\partial\Phi(t_{n},q)}{\partial
t}\rightarrow\frac{\Phi_{n+1,\boldsymbol{j}}-\Phi_{n,\boldsymbol{j}}}{\Delta
t}$ $\displaystyle\frac{\partial}{\partial
q_{i}}(F_{i}(q)\Phi(t_{n},q))\rightarrow\frac{1}{h}\left[(F^{-}_{i}(q_{i+1/2})T^{+}_{i}\Phi_{n})_{\boldsymbol{j}}-(F^{+}_{i}(q_{i-1/2})T^{-}_{i}\Phi_{n,\boldsymbol{j}}+[(F^{+}_{i}(q_{j+1/2})-F^{-}_{i}(q_{j-1/2}))\Phi_{n}]_{\boldsymbol{j}}\right]$

Here $\boldsymbol{j}\equiv(j_{1},...,j_{D})$, $j_{i}=1,...,N$ for $i=1,...,D$,
$n=1,...,N_{t}$, $t_{n}=n\Delta t$ and $q_{i}=j_{i}h$, where $q_{i}$ are the
$i^{th}$ components of the vector $q$. In addition,
$T^{\pm}_{i}\Phi_{n,\boldsymbol{j}}=\Phi_{n,j_{1},...,j_{i\pm 1},...,j_{d}}$,
$(F_{i}(q_{i+1/2}))_{\boldsymbol{j}}=\frac{1}{2}[(F_{i}(q_{j_{1},...,j_{i+1},...,j_{d}})+F_{i}(q_{j_{1},...,j_{i},...,j_{d}})]$,
and $F^{+}(q)=\max(F(q),0)$ and $F^{-}(q)=\min(F(q),0)$.

Define $\lambda=\Delta t/h$. We require

$D\lambda=D\frac{\Delta t}{h}\leq 1$ (157)

for numerical stability.

Then the discretised version of Eq. (V.1) can be rewritten as

$\displaystyle\Phi_{n+1,\boldsymbol{j}}+\lambda\sum_{i=1}^{D}\left(F^{-}_{i}(q_{i+1/2})_{\boldsymbol{j}}T^{+}_{i}\Phi_{n,\boldsymbol{j}}-F^{+}_{i}(q_{i-1/2})_{\boldsymbol{j}}T^{-}_{i}\Phi_{n,\boldsymbol{j}}\right)+[1+\lambda\sum_{i=1}^{D}(F^{+}_{i}(q_{j+1/2})-F^{-}_{i}(q_{j-1/2}))_{\boldsymbol{j}}]\Phi_{n,\boldsymbol{j}}=0$
(158)

with the initial condition (for $n=0$)

$\displaystyle\Phi_{0,\boldsymbol{j}}=\frac{1}{M}\sum_{k=1}^{M}\Pi_{i=1}^{D}\delta_{\omega}(j_{i}h-(X_{0}^{[k]})_{i}).$
(159)

Then the discretised PDEs can be written as a matrix equation

$\displaystyle\mathcal{K}\begin{pmatrix}\Phi_{1,\boldsymbol{j}}\\\
\Phi_{2,\boldsymbol{j}}\\\ \vdots\\\ \Phi_{N_{t}-1,\boldsymbol{j}}\\\
\Phi_{N_{t},\boldsymbol{j}}\end{pmatrix}=\begin{pmatrix}\Phi_{0,\boldsymbol{j}}\\\
0\\\ \vdots\\\ 0\\\ 0\end{pmatrix}$ (160)

and $\mathcal{K}_{\text{ODE}}$ is a $(N_{t}N^{D})\times(N_{t}N^{D})$ Toeplitz
matrix of the form

$\displaystyle\mathcal{K}_{\text{ODE}}=\begin{pmatrix}I&0&0&\ldots&0&0\\\
K&I&0&\ldots&0&0\\\ \vdots&\vdots&\vdots&\vdots&\vdots&\vdots\\\
0&0&0&\ldots&K&I\end{pmatrix}$ (161)

where each $I$ is the $(N^{D})\times(N^{D})$ identity matrix and $K$ is the
$(N^{D})\times(N^{D})$ matrix

$\displaystyle
K=-I+\lambda\sum_{i=1}^{D}\left(F^{-}_{i}(q_{i+1/2})_{\boldsymbol{j}}T^{+}_{i}-F^{+}_{i}(q_{i-1/2})_{\boldsymbol{j}}T^{-}_{i}+(F^{+}_{i}(q_{j+1/2})-F^{-}_{i}(q_{j-1/2}))_{\boldsymbol{j}}I\right).$
(162)

We can then solve for $\Phi_{n,\boldsymbol{j}}$ by matrix inversion

$\displaystyle\begin{pmatrix}\Phi_{1,\boldsymbol{j}}\\\
\Phi_{2,\boldsymbol{j}}\\\ \vdots\\\ \Phi_{N_{t}-1,\boldsymbol{j}}\\\
\Phi_{N_{t},\boldsymbol{j}}\end{pmatrix}=\mathcal{K}^{-1}_{\text{ODE}}\begin{pmatrix}\Phi_{0,\boldsymbol{j}}\\\
0\\\ \vdots\\\ 0\\\ 0\end{pmatrix}$ (163)

Then given the discretised solutions $\Phi_{n,\boldsymbol{j}}$, we can
estimate the ensemble average by

$\displaystyle\langle
A_{\omega}^{n}\rangle=h^{D}\sum_{\boldsymbol{|j|=1}}^{N}A_{\boldsymbol{j}}\Phi_{n,\boldsymbol{j}}$
(164)

where $A_{\boldsymbol{j}}$ is the discretisation of $A(q)$, i.e.
$G(q)\rightarrow A(\boldsymbol{j}h)\equiv A_{\boldsymbol{j}}$. From Section
III.2.2, we see that we also use $\langle A^{\omega}_{n}\rangle$ to estimate
$A_{O}(t_{n})$ with error

$\displaystyle\epsilon_{c}\equiv|A_{o}(t_{n})-\langle
A^{\omega}_{n}\rangle|\leq C(Dh)^{1/3})$ (165)

for $\omega=(Dh)^{1/3}$.

Remark: The stability condition in (157) requires $\Delta t$ to be
proportional to $1/d$. This is the issue for an explicit scheme. This means
$N_{t}$, and consequently the size of matrix $\mathcal{K}$, will be $O(d)$
times larger than a scheme in which $\Delta t$ is independent of $d$. An
implicit scheme for (V.1)–which is usually not needed in classical algorithms
for a transport equation– will remove such an dependence. Since the HHL
algorithm is based on inverting the matrix $\mathcal{K}$, which is basically
implicit anyway, one can starts with an implicit scheme for equation (V.1)
[55]. Since the HHL algorithm just depends on logarithmic of the matrix size,
an $d$ times larger matrix size will not increase significantly the
computational cost, hence we do not explore this issue and leave it for a
future work.

Just like in previous analyses, since our aim is to approximate the observable
$A_{O}(t_{n})$ by using quantum algorithms that require a matrix inversion
subroutine $\mathcal{M}^{-1}\vec{y}$ where $\mathcal{M}$ is Hermitian, we can
define a new Hermitian matrix

$\displaystyle\mathcal{M}_{\text{ODE}}=\begin{pmatrix}0&\mathcal{K}_{\text{ODE}}\\\
\mathcal{K}_{\text{ODE}}^{\dagger}&0\end{pmatrix}$ (166)

which has the same sparsity and condition number as
$\mathcal{K}_{\text{ODE}}$.

The condition number of $\mathcal{M}_{\text{ODE}}$ is $\kappa\lesssim DNT$ and
sparsity is $s\sim D$, from Appendix C with the replacement $d\rightarrow
D/2$.

## References

  * Clader  _et al._ [2013] B. D. Clader, B. C. Jacobs, and C. R. Sprouse, Preconditioned quantum linear system algorithm, Physical review letters 110, 250504 (2013).
  * Childs  _et al._ [2021] A. M. Childs, J.-P. Liu, and A. Ostrander, High-precision quantum algorithms for partial differential equations, Quantum 5, 574 (2021).
  * Costa  _et al._ [2019] P. C. Costa, S. Jordan, and A. Ostrander, Quantum algorithm for simulating the wave equation, Physical Review A 99, 012323 (2019).
  * Linden  _et al._ [2020] N. Linden, A. Montanaro, and C. Shao, Quantum vs. classical algorithms for solving the heat equation, arXiv preprint arXiv:2004.06516 (2020).
  * Engel  _et al._ [2019] A. Engel, G. Smith, and S. E. Parker, Quantum algorithm for the vlasov equation, Physical Review A 100, 062315 (2019).
  * Cao  _et al._ [2013] Y. Cao, A. Papageorgiou, I. Petras, J. Traub, and S. Kais, Quantum algorithm and circuit design solving the poisson equation, New Journal of Physics 15, 013021 (2013).
  * Liu  _et al._ [2021] J.-P. Liu, H. Ø. Kolden, H. K. Krovi, N. F. Loureiro, K. Trivisa, and A. M. Childs, Efficient quantum algorithm for dissipative nonlinear differential equations, Proceedings of the National Academy of Sciences 118 (2021).
  * Lloyd  _et al._ [2020] S. Lloyd, G. De Palma, C. Gokler, B. Kiani, Z.-W. Liu, M. Marvian, F. Tennie, and T. Palmer, Quantum algorithm for nonlinear differential equations, arXiv preprint arXiv:2011.06571 (2020).
  * Leyton and Osborne [2008] S. K. Leyton and T. J. Osborne, A quantum algorithm to solve nonlinear differential equations, arXiv preprint arXiv:0812.4423 (2008).
  * Joseph [2020] I. Joseph, Koopman–von neumann approach to quantum simulation of nonlinear classical dynamics, Physical Review Research 2, 043102 (2020).
  * Dodin and Startsev [2021] I. Y. Dodin and E. A. Startsev, On applications of quantum computing to plasma simulations, Physics of Plasmas 28, 092101 (2021).
  * Jin and Li [2003] S. Jin and X. Li, Multi-phase computations of the semiclassical limit of the Schrödinger equation and related problems: Whitham vs. Wigner, Phys. D 182, 46 (2003).
  * Berry [2014] D. W. Berry, High-order quantum algorithm for solving linear differential equations, Journal of Physics A: Mathematical and Theoretical 47, 105301 (2014).
  * Berry  _et al._ [2017] D. W. Berry, A. M. Childs, A. Ostrander, and G. Wang, Quantum algorithm for linear differential equations with exponentially improved dependence on precision, Communications in Mathematical Physics 356, 1057 (2017).
  * Mishra [2009] S. Mishra, Uncertainty and sensitivity analysis techniques for hydrologic modeling, Journal of hydroinformatics 11, 282 (2009).
  * Fjordholm  _et al._ [2016] U. S. Fjordholm, S. Mishra, and E. Tadmor, On the computation of measure-valued solutions, Acta numerica 25, 567 (2016).
  * Ying and Candes [2006] L. Ying and E. J. Candes, The phase flow method, Journal of Computational Physics 220, 184 (2006).
  * Heller [2006] E. J. Heller, Guided gaussian wave packets, Accounts of Chemical Research 39, 127 (2006).
  * Lax [1973] P. D. Lax,  _Hyperbolic systems of conservation laws and the mathematical theory of shock waves_ (SIAM, 1973).
  * Crandall and Lions [1983] M. G. Crandall and P.-L. Lions, Viscosity solutions of hamilton-jacobi equations, Transactions of the American mathematical society 277, 1 (1983).
  * Osher and Sethian [1988] S. Osher and J. A. Sethian, Fronts propagating with curvature-dependent speed: algorithms based on Hamilton-Jacobi formulations, J. Comput. Phys. 79, 12 (1988).
  * Lasry and Lions [2007] J.-M. Lasry and P.-L. Lions, Mean field games, Japanese journal of mathematics 2, 229 (2007).
  * Whitham [2011] G. B. Whitham,  _Linear and nonlinear waves_ , Vol. 42 (John Wiley & Sons, 2011).
  * Sparber  _et al._ [2003] C. Sparber, P. Markowich, and N. Mauser, Wigner functions versus WKB-methods in multivalued geometrical optics, Asymptot. Anal. 33, 153 (2003).
  * Jin and Osher [2003] S. Jin and S. Osher, A level set method for the computation of multivalued solutions to quasi-linear hyperbolic PDEs and Hamilton-Jacobi equations, Commun. Math. Sci. 1, 575 (2003).
  * Fomel and Sethian [2002] S. Fomel and J. A. Sethian, Fast-phase space computation of multiple arrivals, Proceedings of the National Academy of Sciences 99, 7329 (2002).
  * Engquist and Runborg [2003] B. Engquist and O. Runborg, Computational high frequency wave propagation, Acta numerica 12, 181 (2003).
  * Jin  _et al._ [2011] S. Jin, P. Markowich, and C. Sparber, Mathematical and computational methods for semiclassical schrödinger equations, Acta Numerica 20, 121 (2011).
  * Jin  _et al._ [2005] S. Jin, H. Liu, S. Osher, and Y.-H. R. Tsai, Computing multivalued physical observables for the semiclassical limit of the Schrödinger equation, J. Comput. Phys. 205, 222 (2005).
  * Glassner [1989] A. S. Glassner,  _An introduction to ray tracing_ (Morgan Kaufmann, 1989).
  * Benamou [1996] J.-D. Benamou, Big ray tracing: Multivalued travel time field computation using viscosity solutions of the eikonal equation, Journal of Computational Physics 128, 463 (1996).
  * Engquist and Runborg [1996] B. Engquist and O. Runborg, Multi-phase computations in geometrical optics, Journal of Computational and Applied Mathematics 74, 175 (1996).
  * Cheng  _et al._ [2003] L.-T. Cheng, H. Liu, and S. Osher, Computational high-frequency wave propogation using the level-set method with applications to the semi-classical limit of the schrödinger equations, Communications in Mathematical Sciences 1, 593 (2003).
  * Dutt  _et al._ [2000] A. Dutt, L. Greengard, and V. Rokhlin, Spectral deferred correction methods for ordinary differential equations, BIT Numerical Mathematics 40, 241 (2000).
  * Childs  _et al._ [2017] A. M. Childs, R. Kothari, and R. D. Somma, Quantum algorithm for systems of linear equations with exponentially improved dependence on precision, SIAM Journal on Computing 46, 1920 (2017).
  * Alase  _et al._ [2021] A. Alase, R. R. Nerem, M. Bagherimehrab, P. Høyer, and B. C. Sanders, Tight bound for estimating expectation values from a system of linear equations, arXiv preprint arXiv:2111.10485 (2021).
  * Harrow  _et al._ [2009] A. W. Harrow, A. Hassidim, and S. Lloyd, Quantum algorithm for linear systems of equations, Physical review letters 103, 150502 (2009).
  * Berry  _et al._ [2015] D. W. Berry, A. M. Childs, and R. Kothari, Hamiltonian simulation with nearly optimal dependence on all parameters, in  _2015 IEEE 56th Annual Symposium on Foundations of Computer Science_ (IEEE, 2015) pp. 792–809.
  * Buhrman  _et al._ [2001] H. Buhrman, R. Cleve, J. Watrous, and R. De Wolf, Quantum fingerprinting, Physical Review Letters 87, 167902 (2001).
  * Aharonov  _et al._ [2009] D. Aharonov, V. Jones, and Z. Landau, A polynomial quantum algorithm for approximating the jones polynomial, Algorithmica 55, 395 (2009).
  * Knill  _et al._ [2007] E. Knill, G. Ortiz, and R. D. Somma, Optimal quantum measurements of expectation values of observables, Physical Review A 75, 012328 (2007).
  * Low and Chuang [2019] G. H. Low and I. L. Chuang, Hamiltonian simulation by qubitization, Quantum 3, 163 (2019).
  * Gilyén  _et al._ [2019] A. Gilyén, Y. Su, G. H. Low, and N. Wiebe, Quantum singular value transformation and beyond: exponential improvements for quantum matrix arithmetics, in  _Proceedings of the 51st Annual ACM SIGACT Symposium on Theory of Computing_ (2019) pp. 193–204.
  * Tornberg and Engquist [2004] A.-K. Tornberg and B. Engquist, Numerical approximations of singular source terms in differential equations, Journal of Computational Physics 200, 462 (2004).
  * Lions and Paul [1993] P.-L. Lions and T. Paul, Sur les mesures de wigner, Revista matemática iberoamericana 9, 553 (1993).
  * Gérard  _et al._ [1997] P. Gérard, P. A. Markowich, N. J. Mauser, and F. Poupaud, Homogenization limits and wigner transforms, Communications on Pure and Applied Mathematics: A Journal Issued by the Courant Institute of Mathematical Sciences 50, 323 (1997).
  * Ryzhik  _et al._ [1996] L. Ryzhik, G. Papanicolaou, and J. B. Keller, Transport equations for elastic and other waves in random media, Wave motion 24, 327 (1996).
  * Raviart [1985] P.-A. Raviart, An analysis of particle methods, in  _Numerical methods in fluid dynamics_ (Springer, 1985) pp. 243–324.
  * Majda  _et al._ [2002] A. J. Majda, A. L. Bertozzi, and A. Ogawa, Vorticity and incompressible flow. cambridge texts in applied mathematics, Appl. Mech. Rev. 55, B77 (2002).
  * Monaghan [1992] J. J. Monaghan, Smoothed particle hydrodynamics, Annual review of astronomy and astrophysics 30, 543 (1992).
  * Bird [1994] G. A. Bird, Molecular gas dynamics and the direct simulation of gas flows, Molecular gas dynamics and the direct simulation of gas flows (1994).
  * Sonnendrücker and Kormann [2013] E. Sonnendrücker and K. Kormann, Numerical methods for vlasov equations, Lecture notes (2013).
  * LeVeque [2002] R. J. LeVeque,  _Finite volume methods for hyperbolic problems_ , Vol. 31 (Cambridge university press, 2002).
  * Gilyén  _et al._ [2018] A. Gilyén, Y. Su, G. H. Low, and N. Wiebe, Quantum singular value transformation and beyond: exponential improvements for quantum matrix arithmetics, arXiv preprint arXiv:1806.01838 (2018).
  * Jin  _et al._ [2022] S. Jin, N. Liu, and Y. Yu, Time complexity analysis of quantum difference methods for linear high dimensional and multiscale partial differential equations, arXiv preprint arXiv:2202.04537 (2022).