# Rendering Wireless Environments Useful for Gradient Estimators: A Zero-Order
Stochastic Federated Learning Method

Elissa Mhanna and Mohamad Assaad  Laboratoire des Signaux et Systèmes
Université Paris-Saclay, CNRS, CentraleSupélec
91190 Gif-sur-Yvette, France
Email: {elissa.mhanna<EMAIL_ADDRESS>
###### Abstract

Federated learning (FL) is a novel approach to machine learning that allows
multiple edge devices to collaboratively train a model without disclosing
their raw data. However, several challenges hinder the practical
implementation of this approach, especially when devices and the server
communicate over wireless channels, as it suffers from communication and
computation bottlenecks in this case. By utilizing a communication-efficient
framework, we propose a novel zero-order (ZO) method with a one-point gradient
estimator that harnesses the nature of the wireless communication channel
without requiring the knowledge of the channel state coefficient. It is the
first method that includes the wireless channel in the learning algorithm
itself instead of wasting resources to analyze it and remove its impact. The
two main difficulties of this work are that in FL, the objective function is
usually not convex, which makes the extension of FL to ZO methods challenging,
and that including the impact of wireless channels requires extra attention.
However, we overcome these difficulties and comprehensively analyze the
proposed zero-order federated learning (ZOFL) framework. We establish its
convergence theoretically, and we prove a convergence rate of
$O(\frac{1}{\sqrt[3]{K}})$ in the nonconvex setting. We further demonstrate
the potential of our algorithm with experimental results, taking into account
independent and identically distributed (IID) and non-IID device data
distributions.

##  I Introduction

ZO methods are a subfield of optimization that assume that first-order (FO)
information or access to function gradients is unavailable. ZO optimization is
based on estimating the gradient using function values queried at a certain
number of points. The number of function queries depends on the assumptions of
the problem. For example, in multi-point gradient estimates [1, 2], they
construct the gradient by performing the difference of function values
obtained at many random or predefined points. However, they assume that the
stochastic setting stays the same during all these queries. For example, for
functions $\theta\mapsto f(\theta,S)$ subject to a stochastic variable $S$,
two-point gradient estimates have the form,

$g=d\frac{f(\theta+\gamma\Phi,S)-f(\theta-\gamma\Phi,S)}{2\gamma}\Phi,$

with $\theta\in\mathbb{R}^{d}$ the optimization variable, $\gamma>0$ a small
value, and $\Phi$ a random vector with a symmetric distribution. By contrast,
one-point estimates that use only one function value [3, 4, 5], principally
obtained at a random point,

$g=\frac{d}{\gamma}f(\theta+\gamma\Phi,S)\Phi,$

assume that the settings are continuously changing during optimization. This
is an important property as it resonates with many realistic applications,
like when the optimization is performed in wireless environments or is based
on previous simulation results. Recently, an appeal to ZO optimization is
emerging in the machine-learning community, where optimizers are based on
gradient methods. Examples include reinforcement learning [6, 7], generating
contrastive explanations for black-box classification models [8], and
effecting adversarial perturbations on such models [9, 10].

On the other hand, with the massive amounts of data generated or accessed by
mobile devices, a growing research interest in both sectors of academia and
industry [11] is focused on FL [12, 13], as it’s a practical solution for
training models on such data without the need to log them to a server. A lot
of effort has been invested in developing first-order [12, 14, 15] and second-
order [16, 17] methods to improve the efficacy of FL. These methods typically
require access to the gradient or the Hessian of the local objective functions
in their implementation to solve the optimization problem. However, using and
exchanging such information raises many challenges, such as expensive
communication and computation and privacy concerns [18].

Figure 1: Federated learning over wireless networks.

There is more interest recently in learning over wireless environments [19,
20, 21, 22, 23, 24], with the increase of devices connected to servers through
cellular networks. In this paper, we are interested in this scenario
illustrated in Fig. 1. Similarly to the aforementioned work, we are examining
the case of analog communications between the devices and the server. However,
it is a challenging problem as when the information is sent over the wireless
channel, it becomes subject to a perturbation induced by the channel. This
perturbation is not limited to additive noise, as noise is, in fact, due to
thermal changes at the receiver. The channel acts as a filter for the
transmitted signal [25, 26],

$\hat{x}=Hx+n$ (1)

where $x$ and $\hat{x}\in\mathbb{R}^{d}$ are the sent and received signals,
respectively. $H\in\mathbb{R}^{d\times d}$ is the channel matrix, and
$n\in\mathbb{R}^{d}$ is the additive noise, both of which are stochastic,
constantly changing, and unknown. In general, all these entities are
considered to have complex values. However, as we are interested in sending
only real scalar values and we are not interested in decoding without error
the sent information (but rather on using a perturbing on loss function), the
phase shift introduced by the channel is considered as an additional
perturbation in the channel model. Since our model is based on perturbing the
loss function, considering only real values in the equation above is
sufficient. We elaborate further on the channel modeling and why we can
consider it real in Appendix A for the interested reader. In FL, $x$ may
denote the model or its gradients sent over the channel. To remove this
impact, every channel element must be analyzed and removed to retrieve the
sent information. This analysis is costly in computation and time resources.
Thus, here our objective is to study FL in wireless environments without
wasting resources.

Further, we are interested in exploring the potential of ZO optimization to
deal with some of the difficulties demonstrated by FL. We then consider an FL
setting where a central server coordinates with $N$ edge devices to solve an
optimization problem collaboratively. The data is private to every device and
the exchanges between the server and the devices is restricted to the
optimization parameters. To that end, let $\mathcal{N}=\\{1,...,N\\}$ be the
set of devices and $\theta\in\mathbb{R}^{d}$ denote the global model. We
define $F_{i}:\mathbb{R}^{d}\rightarrow\mathbb{R}$ as the loss function
associated with the local data stored on device $i$, $\forall
i\in\mathcal{N}$. The objective is to minimize the function
$F:\mathbb{R}^{d}\rightarrow\mathbb{R}$ that is composed of the said devices’
loss functions, such that

$\min_{\theta\in\mathbb{R}^{d}}F(\theta):=\sum_{i=1}^{N}F_{i}(\theta)\;\;\text{with}\;\;F_{i}(\theta)=\mathbb{E}_{S_{i}\sim
D_{i}}f_{i}(\theta,S_{i}).$ (2)

$S_{i}$ is an i.i.d. ergodic stochastic process following a local distribution
$D_{i}$. $S_{i}$ is used to model various stochastic perturbation, e.g. local
data distribution among others. We further consider the case where the devices
do not have access to their gradients for computational and communication
restraints, and they must estimate this gradient by querying their model only
once per update. They obtain a scalar value from this query (rather than a
long vector as in the gradient), that they must send back to the server, which
saves computation and communication resources since only a scalar is computed
instead of a long vector of gradient/model as in the standard gradient
approach.

###  I-A Motivation for Our Work

Communication Bottleneck. In general, the main idea of FL is that the devices
receive the model from the server, make use of their data to update the
gradient, and then send back their gradients without ever disclosing their
data. The server then updates the model using the collected and averaged
gradients, and the process repeats. Since the gradients have the same
dimension as the model, in every uplink step, there are $Nd$ values that need
to be uploaded, which forms the fundamental communication bottleneck in FL. To
deal with this issue, some propose local multiple gradient descent steps to be
done by the devices before sending their gradients back to the server to save
communication resources [27], or allow partial device participation at every
iteration [28, 29], or both [12]. Others propose lossy compression of the
gradient before uploading it to the server. For example, all [30, 31, 32]
suggest stochastic unbiased quantization approaches, where gradients are
approximated with a finite set of discrete values for efficiency. [33]
proposes the quantization of gradient differences of the current and previous
iterations, allowing the update to incorporate new information, while [34]
proposes the sparsification of this difference. Sparsification means that if a
vector component is not large enough, it will not be transmitted.

Channel Impact. In FL over wireless channels, there is a problem with channel
knowledge. When the devices upload their gradient $g\in\mathbb{R}^{d}$ to the
server, the server receives $Hg+n$ as shown in equation (1). In [19] and [35]
and all references within, they assume that they can remove the impact of the
channel. However, as the channel matrix $H$ coefficients follow a stochastic
process and there are two unknown received entities, the channel $H$ and the
gradient, the knowledge of the gradient requires estimating the channel
coefficients at each iteration of the FL. This requires computation resources,
and more importantly, it requires resources to exchange control/reference
signals between the devices and the server at each time/iteration to estimate
the channel coefficients $H$. Alternatively, our work offers a much simpler
approach. We do not waste resources trying to analyze the channel. We use the
channel in the implementation itself. It is part of the learning. We harness
it to construct our gradient estimate without the need to remove its impact,
saving both computation and communication resources.

Computation Demands. Unlike standard methods that rely on the computational
capabilities of participating devices, our approach is less demanding. Devices
simply receive the global model, query it with their data, and send back the
scalar loss, eliminating the need for "backward pass" computation. Only the
"forward pass" is performed.

Black-Box Optimization in FL. One motivation for employing ZO methods is
black-box problems [35] when gradient information cannot be acquired or is
complicated to compute. For example, in hyperparameter tuning, gradients
cannot be calculated, as there isn’t an analytic relationship between the loss
function and the hyperparameters [36].

###  I-B Challenges and Contribution

Addressing nonconvexity in FL is challenging. Our ZO method must handle
nonconvexity, noise, and stochasticity efficiently, which can slow down
convergence in gradient techniques. Additionally, the channel’s introduction
adds uncertainty and constraints on the number of communication exchanges. We
need to ensure consistent and reliable performance, considering unknown
probability distributions and fewer function evaluations. Unlike convex cases,
nonconvex optimization does not allow easy quantification of optimization
progress. Verifying gradient convergence becomes intricate due to biased
gradient estimates. Moreover, unbounded variance in one-point estimates can
lead to significant gradient deviations. These challenges involve technical
and intuitive complexities we navigate.

In this work, we overcome these difficulties and propose a new communication-
efficient algorithm in the nonconvex setting. This algorithm differs from the
standard gradient method as it entails two reception-update steps instead of
one, and it is not a simple extension of FO to ZO where the devices still have
to upload their full model/gradient, as is the case in [35]. By limiting the
exchange to scalar-valued updates, we counter the communication bottleneck,
and we save up to a factor of $O(d)$ per communication round, in comparison to
standard methods, in terms of total exchanges of variables between the devices
and the server, saving a lot of execution time. We harness the nature of
autocorrelated channels for truly "blind" reception of the data. We prove the
convergence theoretically with a one-point estimate and provide experimental
evidence. An important distinction worth noting is that standard ZO methods
establish convergence by focusing on the expected convergence of the exact
gradient. In contrast to prior research, our approach goes further in the
proof. We demonstrate the convergence of the exact gradient itself almost
surely, not solely its expected value. The key element in this proof is
employing Doob’s martingale inequality to constrain the stochastic error
resulting from the estimated gradient. We finally extend the analysis to non-
symmetrical channel models, i.e., channels without zero-mean, and thus provide
a practical algorithm for general settings.

##  II Algorithm

This section illustrates our proposed zero-order stochastic federated learning
algorithm with a new gradient estimator (ZOFL).

###  II-A The 1P-ZOFL Algorithm

Algorithm 1 The 1P-ZOFL Algorithm

Input: Initial model $\theta_{0}\in\mathbb{R}^{d}$, the initial step-sizes
$\alpha_{0}$ and $\gamma_{0}$, and the channels’ standard deviation
$\sigma_{h}$

1: for $k=0,2,4,...$ do

2: The server receives $\sum_{j=1}^{N}\frac{h_{j,k}}{\sigma_{h}^{2}}+n_{j,k}$

3: The server broadcasts
$\theta_{k}+\gamma_{k}\Phi_{k}\sum_{j=1}^{N}\Big{(}\frac{h_{j,k}}{\sigma_{h}^{2}}+n_{j,k}\Big{)}$
to all devices

4: The server
receives$\sum_{i=1}^{N}h_{i,k+1}\tilde{f}_{i}\Big{(}\theta_{k}+\gamma_{k}\Phi_{k}\sum_{j=1}^{N}\big{(}\frac{h_{j,k}}{\sigma_{h}^{2}}+n_{j,k}\big{)},S_{i,k+1}\Big{)}+n_{i,k+1}$

5: The server multiplies the received scalar sum by $\Phi_{k}$ to assemble
$g_{k}$ given in (3)

6: The server updates $\theta_{k+1}=\theta_{k}-\alpha_{k}g_{k}$

7: end for

We consider an intermediary wireless environment between the server and each
device $i$ for $i\in\mathcal{N}$ as shown in Fig. 1. Wireless channels
introduce a stochastic scaling on the sent signal as elaborated in equation
(1). As we only send a scalar value over the channel at a time, our channel
has only one scalar coefficient in addition to a scalar noise. Channel
coefficients are usually autocorrelated from one timeslot to the next. Let
$h_{i,k}$ denote the channel scaling affecting the sent signal from device $i$
to the server at timeslot $k$, independent from all other devices’ channels.
We assume $h_{i,k}$ to be a zero-mean random variable with standard deviation
$\sigma_{h}$, $\forall i\in\mathcal{N},\forall k\in\mathbb{N}^{+}$, and
$n_{i,k}$ an additive noise on the transmitted signal. Assuming that the
channel is time-correlated for two consecutive iterations $k$ and $k+1$, such
that the autocovariance is $\mathbb{E}[h_{i,k}h_{i,k+1}]=K_{hh}$, $\forall
i\in\mathcal{N}$, $\forall k\in\mathbb{N}^{+}$, we present our learning method
in Algorithm 1:

The devices must carry out two communication steps. In the first, every device
sends the value $\frac{1}{\sigma_{h}^{2}}$ to the server. According to
equation (1), the server receives $\frac{h_{j,k}}{\sigma_{h}^{2}}+n_{j,k}$
from every device $j$. Hence, it receives the sum in step $2$. Afterward, the
server uses the values received to adjust the model and broadcasts it to the
devices. When device $i$ receives the model, it receives
$h_{i,k+1}^{DL}[\theta_{k}+\gamma_{k}\Phi_{k}\sum_{j=1}^{N}(\frac{h_{j,k}}{\sigma_{h}^{2}}+n_{j,k})]+n_{i,k+1}^{DL}$,
and to simplify notation, we let the stochastic vector
$[h_{i,k+1}^{DL},n_{i,k+1}^{DL}]$ be included within the big vector
$S_{i,k+1}$ of stochastic perturbations. Device $i$ then queries this received
model to obtain the stochastic loss $f_{i}$. Then the devices send $f$ to the
server in the second communication step, and according to equation (1), the
server receives the quantity indicated in step $4$. Finally, the server
assembles the gradient estimate and is able to update $\theta$ according to
step $7$. All transmissions are subject to channel scaling and additive noise.
We designate them by $h$ and $n$ in the device-to-server transmission. In the
server-to-devices one, we designate them by $S$. We let
$\tilde{f}_{i}=\frac{f_{i}}{\sigma_{h}^{2}}$ be the normalized loss function
and define $\alpha_{k}$ and $\gamma_{k}$ as two step-sizes and
$\Phi_{k}\in\mathbb{R}^{d}$ as a perturbation vector generated by the server
that has the same dimension as that of the model.

We emphasize here that $g_{k}$ (in step $6$) is the gradient estimate in this
case, and one can see that the impact of the channel is included in the
gradient estimate and hence in the learning. The major advantage of this
algorithm is that each device sends only two scalar values. This is stark
improvement in communication efficiency over standard FL algorithms that
require each device to send back the whole model or local gradient of
dimension $d$. In effect, it is resource draining and can be unrealistic to
assume it is possible.

###  II-B The Estimated Gradient

We provide here analysis of our ZO gradient estimate. We propose the one-point
estimate:

$g_{k}=\Phi_{k}\sum_{i=1}^{N}\Big{[}h_{i,k+1}\tilde{f}_{i}\Big{(}\theta_{k}^{\prime},S_{i,k+1}\Big{)}+n_{i,k+1}\Big{]},$
(3)

where
$\theta_{k}^{\prime}=\theta_{k}+\gamma_{k}\Phi_{k}\sum_{j=1}^{N}\big{(}\frac{h_{j,k}}{\sigma_{h}^{2}}+n_{j,k}\big{)}$.
The values $h_{i,k}$, $h_{i,k+1}$, and the noise remain unknown. This saves
computation complexity and is very communication efficient as it transcends
the need to send pilot signals to estimate the channel continuously.

We next consider the following assumptions on the additive noise, the
perturbation vector, and the local loss functions.

###### Assumption 1

$n_{i,k}$ is assumed to be a zero-mean uncorrelated noise with bounded
variance, meaning $E(n_{i,k})=0$ and $E(n_{i,k}^{2})=\sigma_{n}^{2}<\infty$,
$\forall i\in\mathcal{N}$, $\forall k\in\mathbb{N}^{+}$. For any timeslot $k$,
$E(n_{i,k}n_{j,k})=0$ if $i\neq j$. For any device $i$,
$E(n_{i,k}n_{i,k^{\prime}})=0$ if $k\neq k^{\prime}$.

###### Assumption 2

Let $\Phi_{k}=(\phi_{k}^{1},\phi_{k}^{2},\ldots,\phi_{k}^{d})^{T}$. At each
iteration $k$, the server generates its $\Phi_{k}$ vector independently from
other iterations. In addition, the elements of $\Phi_{k}$ are assumed i.i.d
with $\mathbb{E}(\phi_{k}^{d_{1}}\phi_{k}^{d_{2}})=0$ for $d_{1}\neq d_{2}$
and there exists $\alpha_{2}>0$ such that
$\mathbb{E}(\phi_{k}^{d_{j}})^{2}=\alpha_{2}$, $\forall{d_{j}}$, $\forall k$.
We further assume there exists a constant $\alpha_{3}>0$ where
$\|\Phi_{k}\|\leq\alpha_{3}$, $\forall k$.

###### Example 1

An example of a perturbation vector satisfying Assumption 2, is picking every
dimension of $\Phi_{k}$ from $\\{-\frac{1}{\sqrt{d}},\frac{1}{\sqrt{d}}\\}$
with equal probability. Then, $\alpha_{2}=\frac{1}{d}$ and $\alpha_{3}=1$.

###### Assumption 3

All loss functions $\theta\mapsto f_{i}(\theta,S_{i})$ are Lipschitz
continuous with Lipschitz constant $L_{S_{i}}$,
$|f_{i}(\theta,S_{i})-f_{i}(\theta^{\prime},S_{i})|\leq
L_{S_{i}}\|\theta-\theta^{\prime}\|$, $\forall i\in\mathcal{N}$. In addition,
$\mathbb{E}_{S_{i}}f_{i}(\theta,S_{i})<\infty,\forall i\in\mathcal{N}$.

Let
$\mathcal{H}_{k}=\\{\theta_{0},S_{0},\theta_{1},S_{1},...,\theta_{k},S_{k}\\}$
denote the history sequence, then the following two lemmas characterize our
gradient estimate.

###### Lemma 1

Let Assumptions 1 and 2 be satisfied and define the scalar value
$c_{1}=\alpha_{2}\frac{K_{hh}}{\sigma_{h}^{4}}$, then the gradient estimator
is biased w.r.t. the objective function’s exact gradient $\nabla F(\theta)$.
Concretely, $\mathbb{E}[g_{k}|\mathcal{H}_{k}]=c_{1}\gamma_{k}(\nabla
F(\theta_{k})+b_{k})$ $\forall k\in\mathbb{N}^{+}$, where $b_{k}$ is the bias
term.

Proof: Refer to Appendix B-A.

###### Lemma 2

Let Assumptions 1-3 and the inequality $\|\theta_{k}\|<\infty$ hold almost
surely. There exist a bounded constant $c_{2}>0$, such that
$\mathbb{E}[\|g_{k}\|^{2}|\mathcal{H}_{k}]\leq c_{2}$ almost surely.

Proof: Refer to Appendix B-B.

##  III Convergence analysis

We begin by assuming that a global minimizer $\theta^{*}\in\mathbb{R}^{d}$
exists such that
$\min_{\theta\in\mathbb{R}^{d}}F(\theta)=F(\theta^{*})>-\infty$ and $\nabla
F(\theta^{*})=0$.

###### Assumption 4

We assume the existence and the continuity of $\nabla F_{i}(\theta)$ and
$\nabla^{2}F_{i}(\theta)$, and that there exists a constant $\alpha_{1}>0$
such that $\|\nabla^{2}F_{i}(\theta)\|_{2}\leq\alpha_{1}$,$\forall
i\in\mathcal{N}$.

###### Lemma 3

By Assumption 4, we know that the objective function $\theta\longmapsto
F(\theta)$ is $L$-smooth for some positive constant $L$, $\|\nabla
F(\theta)-\nabla F(\theta^{\prime})\|\leq
L\|\theta-\theta^{\prime}\|,\;\forall\theta,\theta^{\prime}\in\mathbb{R}^{d},$
or equivalently, $F(\theta)\leq F(\theta^{\prime})+\langle\nabla
F(\theta^{\prime}),\theta-\theta^{\prime}\rangle+\frac{L}{2}\|\theta-\theta^{\prime}\|^{2}$.

###### Lemma 4

By Assumptions 1, 2, and 4, we can find a scalar value $c_{3}>0$ such that
$\|b_{k}\|\leq c_{3}\gamma_{k}$.

Proof: Refer to Appendix B-C.

###  III-A 1P-ZOFL convergence

In Lemma 1, we see that in expectation, our estimator deviates from the
gradient direction by the bias term. To provide that this term does not grow
larger and preferably grows smaller as the algorithm evolves, we impose that
$\gamma_{k}$ vanishes. Additionally, to ensure that the expected norm squared
of the estimator, as shown in Lemma 2, does not accumulate residual constant
terms, we impose that the step size $\alpha_{k}$ vanishes. The series
properties in the following assumption come from the recursive analysis of the
algorithm.

###### Assumption 5

Both the step sizes $\alpha_{k}$ and $\gamma_{k}$ vanish to zero as
$k\rightarrow\infty$ and the following series composed of them satisfy the
convergence assumptions $\sum_{k=0}^{\infty}\alpha_{k}\gamma_{k}=\infty$,
$\sum_{k=0}^{\infty}\alpha_{k}\gamma_{k}^{3}<\infty$, and
$\sum_{k=0}^{\infty}\alpha_{k}^{2}<\infty$.

###### Example 2

To satisfy Assumption 5, we consider the following form of the step sizes,
$\alpha_{k}=\alpha_{0}(1+k)^{-\upsilon_{1}}$ and
$\gamma_{k}=\gamma_{0}(1+k)^{-\upsilon_{2}}$ with
$\upsilon_{1},\upsilon_{2}>0$. Then, it is sufficient to find $\upsilon_{1}$
and $\upsilon_{2}$ such that $0<\upsilon_{1}+\upsilon_{2}\leq 1$,
$\upsilon_{1}+3\upsilon_{2}>1$, and $\upsilon_{1}>0.5$.

We next define the stochastic error $e_{k}$ as the difference between the
value of a single realization of $g_{k}$ and its conditional expectation given
the history sequence, i.e., $e_{k}=g_{k}-\mathbb{E}[g_{k}|\mathcal{H}_{k}]$.

The study of this noise and how it evolves is essential for the analysis of
the algorithm as it gives access to the exact gradient when examining the
algorithm’s convergence behavior and permits us to prove that, in fact, the
exact gradient converges to zero and not just the expectation of the exact
gradient. This is a stronger convergence property, and it has not been done
before in ZO nonconvex optimization to the best of our knowledge. The trick is
to show that $e_{k}$ is a martingale difference sequence and to apply Doob’s
martingale inequality to derive the following lemma.

###### Lemma 5

If all Assumptions 1-5 hold and $\|\theta_{k}\|<\infty$ almost surely, then
for any constant $\nu>0$, we have
$\mathbb{P}(\lim_{K\rightarrow\infty}\sup_{K^{\prime}\geq
K}\|\sum_{k=K}^{K^{\prime}}\alpha_{k}e_{k}\|\geq\nu)=0$.

Proof: Refer to Appendix C-A.

The smoothness inequality allows for the first main result, leading to the
second in the following theorem.

###### Theorem 1

When Assumptions 1-5 hold, we have $\sum_{k}\alpha_{k}\gamma_{k}\|\nabla
F(\theta_{k})\|^{2}<+\infty$ and $\lim_{k\rightarrow\infty}\|\nabla
F(\theta_{k})\|=0$ almost surely, meaning that the algorithm converges.

Proof: Refer to Appendix C-B.

Proof sketch: We substitute the algorithm’s updates in the second inequality
of Lemma 3 and replace the estimate by its expectation plus the stochastic
error. We then perform a recursive addition over the iterations $k>0$. With
Lemma 5, the conditions on the step sizes, and the upper bound estimate’s
squared norm, we are able to find an upper bound on
$\sum_{k}\alpha_{k}\gamma_{k}\|\nabla F(\theta_{k})\|^{2}$ when $k$ grows to
$\infty$. The next step is to consider the hypothesis
$\lim_{k\rightarrow\infty}\sup\|\nabla F(\theta_{k})\|\geq\rho$, for $\rho>0$,
and prove that it contradicts with the first result.

Define $\delta_{k}=F(\theta_{k})-F(\theta^{*})$. We next find an upper bound
on the convergence rate of Algorithm 1.

###### Theorem 2

Consider in addition to the assumptions in Theorem 1, that the step sizes are
those of Example 2 with $\upsilon_{3}=\upsilon_{1}+\upsilon_{2}<1$. Then, we
can write

$\begin{split}\frac{\sum_{k}\alpha_{k}\gamma_{k}\mathbb{E}\big{[}\|\nabla
F(\theta_{k})\|^{2}\big{]}}{\sum_{k}\alpha_{k}\gamma_{k}}\leq\frac{A(1-\upsilon_{3})}{(K+2)^{1-\upsilon_{3}}-1}.\\\
\end{split}$ (4)

with
$A=\frac{2\delta_{0}}{c_{1}\alpha_{0}\gamma_{0}}+c_{3}^{2}\gamma_{0}^{2}\big{(}\frac{\upsilon_{1}+3\upsilon_{2}}{\upsilon_{1}+3\upsilon_{2}-1}\big{)}+\frac{Lc_{2}\alpha_{0}}{c_{1}\gamma_{0}}\big{(}\frac{2\upsilon_{1}}{2\upsilon_{1}-1}\big{)}$.

Proof: Refer to Appendix C-C.

In Theorem 2, we see that the optimal choice of the exponents for the time-
varying component $O(\frac{1}{K^{1-\upsilon_{1}-\upsilon_{2}}})$, is
$\upsilon_{1}=\frac{1}{2}$ and $\upsilon_{2}=\frac{1}{6}$ for a rate of
$O(\frac{1}{\sqrt[3]{K}})$. However, to prevent the constant component from
growing too large, it is recommended to choose slightly larger exponents of
$\upsilon_{1}=\frac{1}{2}+\frac{\epsilon}{2}$ and
$\upsilon_{2}=\frac{1}{6}+\frac{\epsilon}{2}$, where $\epsilon$ is a small
strictly positive value. This will result in a rate of
$O\big{(}\frac{1}{K^{\frac{1}{3}-\epsilon}}\big{)}$.

###  III-B Non-Symmetrical Channels Case

Assuming a non-symmetrical channel model with $\mathbb{E}[h_{i,k}]=\mu_{h}$
and $\sigma_{h}^{2}=\mathbb{E}[h_{i,k}^{2}]-\mu_{h}^{2}$, $\forall i,\forall
k$, we provide how our gradient estimates and algorithms can be adjusted in
Appendix D to account for this case. In fact, non-symmetrical channel models
(e.g., Rician) offer a simplification of both analysis and implementation in
comparison to symmetrical models (e.g., Rayleigh), as the non-zero mean no
longer cancels out the gradient, and the design is further independent of the
autocorrelation of the channels.

##  IV Experimental results

For the experimental results, we test our algorithms in nonconvex binary image
classification problems, and we compare them against the original FL algorithm
FedAvg [12] with exact gradient and one local update per round. However, we do
not consider the effect of the channel or any noise/stochasticity for the
FedAvg algorithm. All experiments are done for $100$ devices and data batches
of $10$ images per user per round. Every communication round in the graphs
include all steps $2$ through $7$ for Algorithm 1.

We classify photos of the two digits “$0$” and “$1$” from the MNIST dataset
[37] using a nonconvex logistic regression model with a regularization
parameter of $0.001$. All images are divided equally among the devices and are
considered to be preprocessed by being compressed to have dimension $d=10$
using a lossy autoencoder. We run our code on $50$ simulations with different
random model initializations testing the accuracy in every iteration against
an independent test set. The graphs in Fig. 2 are averaged over all these
simulations. For the non-IID data distribution, we first sort the images
according to their labels and then divide them among the devices. $\Phi_{k}$
is generated according to Example 1. All channels are generated using the
normal distribution with autocovariance $K_{hh}=\frac{1}{2}$. The noise is
Gaussian with $0$ mean and variance $\sigma_{n}^{2}=\frac{1}{4}$. The step
size for FedAvg is taken as $\eta=0.15$. For 1P-ZOFL,
$\alpha_{k}=0.5(1+k)^{-0.51}$ and $\gamma_{k}=2.5(1+k)^{-0.18}$. Our algorithm
performs consistently well with all the different random variations
influencing every simulation in line with our theoretical results. Considering
non-IID data distribution seems to be better impacted by our ZO algorithm
slightly at the beginning without a major effect on the final result.

In Fig. 3, we test our algorithm with different noise variance
$\sigma_{n}^{2}$ and notice that while the algorithm is slowed slightly at the
beginning, it does not affect the convergence. When $\sigma_{n}^{2}=2.25$, we
change $\alpha_{0}=0.1$ and $\gamma_{0}=0.8$. When $\sigma_{n}^{2}=10.0489$,
$\alpha_{0}=0.07$ and $\gamma_{0}=0.3$. We provide all other details alongside
a discussion of our algorithm’s performance in relation to the noise variance
in Appendix E.

Figure 2: Accuracy evolution of 1P-ZOFL vs. FedAvg for IID data and non-IID
distribution.

Figure 3: Accuracy evolution of 1P-ZOFL for
$\sigma_{n}^{2}=\\{0.25,1,2.25,10.0489\\}$.

##  V Conclusion

This work considers a learning problem over wireless channels and proposes a
new zero-order federated learning method with a one-point gradient estimator.
We limit the communication to scalar-valued feedback from the devices and
incorporate the wireless channel in the learning algorithm. We provide
theoretical and experimental evidence for convergence and find an upper bound
on the convergence rate.

## References

  * [1] J. C. Duchi, M. I. Jordan, M. J. Wainwright, and A. Wibisono, “Optimal rates for zero-order convex optimization: The power of two function evaluations,” _IEEE Transactions on Information Theory_ , vol. 61, no. 5, pp. 2788–2806, 2015.
  * [2] A. Agarwal, O. Dekel, and L. Xiao, “Optimal algorithms for online convex optimization with multi-point bandit feedback,” in _COLT_ , 2010.
  * [3] A. Flaxman, A. T. Kalai, and H. B. McMahan, “Online convex optimization in the bandit setting: gradient descent without a gradient,” _CoRR_ , vol. cs.LG/0408007, 2004. [Online]. Available: http://arxiv.org/abs/cs.LG/0408007
  * [4] W. Li and M. Assaad, “Distributed stochastic optimization in networks with low informational exchange,” _IEEE Transactions on Information Theory_ , vol. 67, no. 5, pp. 2989–3008, 2021.
  * [5] E. Mhanna and M. Assaad, “Zero-order one-point estimate with distributed stochastic gradient-tracking technique,” 2022. [Online]. Available: https://arxiv.org/abs/2210.05618
  * [6] A. Vemula, W. Sun, and J. Bagnell, “Contrasting exploration in parameter and action space: A zeroth-order optimization perspective,” in _Proceedings of the Twenty-Second International Conference on Artificial Intelligence and Statistics_ , vol. 89. PMLR, 16–18 Apr 2019, pp. 2926–2935.
  * [7] D. Malik, A. Pananjady, K. Bhatia, K. Khamaru, P. Bartlett, and M. Wainwright, “Derivative-free methods for policy optimization: Guarantees for linear quadratic systems,” in _Proceedings of the Twenty-Second International Conference on Artificial Intelligence and Statistics_ , vol. 89. PMLR, 16–18 Apr 2019, pp. 2916–2925.
  * [8] A. Dhurandhar, T. Pedapati, A. Balakrishnan, P.-Y. Chen, K. Shanmugam, and R. Puri, “Model agnostic contrastive explanations for structured data,” 2019\.
  * [9] A. Ilyas, L. Engstrom, A. Athalye, and J. Lin, “Black-box adversarial attacks with limited queries and information,” in _Proceedings of the 35th International Conference on Machine Learning_ , vol. 80. PMLR, 10–15 Jul 2018, pp. 2137–2146.
  * [10] X. Chen, S. Liu, K. Xu, X. Li, X. Lin, M. Hong, and D. Cox, “Zo-adamm: Zeroth-order adaptive momentum method for black-box optimization,” in _Advances in Neural Information Processing Systems_ , vol. 32. Curran Associates, Inc., 2019.
  * [11] K. Bonawitz, H. Eichner, W. Grieskamp, D. Huba, A. Ingerman, V. Ivanov, C. Kiddon, J. Konečný, S. Mazzocchi, B. McMahan, T. Van Overveldt, D. Petrou, D. Ramage, and J. Roselander, “Towards federated learning at scale: System design,” in _Proceedings of Machine Learning and Systems_ , vol. 1, 2019, pp. 374–388.
  * [12] B. McMahan, E. Moore, D. Ramage, S. Hampson, and B. A. y. Arcas, “Communication-Efficient Learning of Deep Networks from Decentralized Data,” in _Proceedings of the 20th International Conference on Artificial Intelligence and Statistics_ , vol. 54. PMLR, 20–22 Apr 2017, pp. 1273–1282.
  * [13] P. Kairouz, H. B. McMahan, B. Avent, A. Bellet, M. Bennis, A. N. Bhagoji, K. Bonawit, Z. Charles, G. Cormode, R. Cummings, R. G. L. D’Oliveira, H. Eichner, S. El Rouayheb, D. Evans, J. Gardner, Z. Garrett, A. Gascón, B. Ghazi, P. B. Gibbons, M. Gruteser, Z. Harchaoui, C. He, L. He, Z. Huo, B. Hutchinson, J. Hsu, M. Jaggi, T. Javidi, G. Joshi, M. Khodak, J. Konecný, A. Korolova, F. Koushanfar, S. Koyejo, T. Lepoint, Y. Liu, P. Mittal, M. Mohri, R. Nock, A. Özgür, R. Pagh, H. Qi, D. Ramage, R. Raskar, M. Raykova, D. Song, W. Song, S. U. Stich, Z. Sun, A. Theertha Suresh, F. Tramèr, P. Vepakomma, J. Wang, L. Xiong, Z. Xu, Q. Yang, F. X. Yu, H. Yu, and S. Zhao, 2021.
  * [14] X. Zhang, M. Hong, S. Dhople, W. Yin, and Y. Liu, “Fedpd: A federated learning framework with adaptivity to non-iid data,” _IEEE Transactions on Signal Processing_ , vol. 69, pp. 6055–6070, 2021.
  * [15] J. Wang, Q. Liu, H. Liang, G. Joshi, and H. V. Poor, “A novel framework for the analysis and design of heterogeneous federated learning,” _IEEE Transactions on Signal Processing_ , vol. 69, pp. 5234–5249, 2021.
  * [16] A. Elgabli, C. B. Issaid, A. S. Bedi, K. Rajawat, M. Bennis, and V. Aggarwal, “FedNew: A communication-efficient and privacy-preserving Newton-type method for federated learning,” in _Proceedings of the 39th International Conference on Machine Learning_ , vol. 162. PMLR, 17–23 Jul 2022, pp. 5861–5877.
  * [17] T. Li, A. K. Sahu, M. Zaheer, M. Sanjabi, A. Talwalkar, and V. Smithy, “Feddane: A federated newton-type method,” in _2019 53rd Asilomar Conference on Signals, Systems, and Computers_ , 2019, pp. 1227–1231.
  * [18] T. Li, A. K. Sahu, A. Talwalkar, and V. Smith, “Federated learning: Challenges, methods, and future directions,” _IEEE Signal Processing Magazine_ , vol. 37, no. 3, pp. 50–60, 2020.
  * [19] K. Yang, T. Jiang, Y. Shi, and Z. Ding, “Federated learning via over-the-air computation,” _IEEE Transactions on Wireless Communications_ , vol. 19, no. 3, pp. 2022–2035, 2020.
  * [20] M. M. Amiri and D. Gündüz, “Federated learning over wireless fading channels,” _IEEE Transactions on Wireless Communications_ , vol. 19, no. 5, pp. 3546–3557, 2020.
  * [21] T. Sery and K. Cohen, “On analog gradient descent learning over multiple access fading channels,” _IEEE Transactions on Signal Processing_ , vol. 68, pp. 2897–2911, 2020.
  * [22] H. Guo, A. Liu, and V. K. N. Lau, “Analog gradient aggregation for federated learning over wireless networks: Customized design and convergence analysis,” _IEEE Internet of Things Journal_ , vol. 8, no. 1, pp. 197–210, 2021.
  * [23] T. Sery, N. Shlezinger, K. Cohen, and Y. C. Eldar, “Over-the-air federated learning from heterogeneous data,” _IEEE Transactions on Signal Processing_ , vol. 69, pp. 3796–3811, 2021.
  * [24] Y. Sun, S. Zhou, Z. Niu, and D. Gündüz, “Dynamic scheduling for over-the-air federated edge learning with energy constraints,” _IEEE Journal on Selected Areas in Communications_ , vol. 40, no. 1, pp. 227–242, 2022.
  * [25] D. Tse and P. Viswanath, _Fundamentals of Wireless Communication_. Cambridge University Press, 2005.
  * [26] E. Björnson and L. Sanguinetti, “Making cell-free massive mimo competitive with mmse processing and centralized implementation,” _IEEE Transactions on Wireless Communications_ , vol. 19, no. 1, pp. 77–90, 2020.
  * [27] A. Khaled, K. Mishchenko, and P. Richtarik, “Tighter theory for local sgd on identical and heterogeneous data,” in _Proceedings of the Twenty Third International Conference on Artificial Intelligence and Statistics_ , vol. 108\. PMLR, 26–28 Aug 2020, pp. 4519–4529.
  * [28] T. Chen, G. Giannakis, T. Sun, and W. Yin, “Lag: Lazily aggregated gradient for communication-efficient distributed learning,” in _Advances in Neural Information Processing Systems_ , vol. 31. Curran Associates, Inc., 2018. [Online]. Available: https://proceedings.neurips.cc/paper/2018/file/feecee9f1643651799ede2740927317a-Paper.pdf
  * [29] M. M. Amiri, D. Gündüz, S. R. Kulkarni, and H. V. Poor, “Convergence of update aware device scheduling for federated learning at the wireless edge,” _IEEE Transactions on Wireless Communications_ , vol. 20, no. 6, pp. 3643–3658, 2021.
  * [30] J. Konečný, H. B. McMahan, F. X. Yu, P. Richtárik, A. T. Suresh, and D. Bacon, “Federated learning: Strategies for improving communication efficiency,” 2016. [Online]. Available: https://arxiv.org/abs/1610.05492
  * [31] S. Khirirat, H. R. Feyzmahdavian, and M. Johansson, “Distributed learning with compressed gradients,” 2018. [Online]. Available: https://arxiv.org/abs/1806.06573
  * [32] A. Elgabli, J. Park, A. S. Bedi, M. Bennis, and V. Aggarwal, “Q-gadmm: Quantized group admm for communication efficient decentralized machine learning,” in _ICASSP 2020 - 2020 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP)_ , 2020, pp. 8876–8880.
  * [33] K. Mishchenko, E. Gorbunov, M. Takáč, and P. Richtárik, “Distributed learning with compressed gradient differences,” 2019. [Online]. Available: https://arxiv.org/abs/1901.09269
  * [34] Y. Chen, R. S. Blum, M. Takáč, and B. M. Sadler, “Distributed learning with sparsified gradient differences,” _IEEE Journal of Selected Topics in Signal Processing_ , vol. 16, no. 3, pp. 585–600, 2022.
  * [35] W. Fang, Z. Yu, Y. Jiang, Y. Shi, C. N. Jones, and Y. Zhou, “Communication-efficient stochastic zeroth-order optimization for federated learning,” _IEEE Transactions on Signal Processing_ , vol. 70, pp. 5058–5073, 2022.
  * [36] Z. Dai, B. K. H. Low, and P. Jaillet, “Federated bayesian optimization via thompson sampling,” in _Advances in Neural Information Processing Systems_ , vol. 33. Curran Associates, Inc., 2020, pp. 9687–9699.
  * [37] Y. LeCun and C. Cortes, “The MNIST database of handwritten digits,” 2005.
  * [38] J. L. Doob, “Stochastic processes,” 1953.

##  Appendix A On the channel model

As described in Section 2.2 of Chapter 2 of [25] and considering a double-
sideband suppressed carrier amplitude modulation (DSB-SC) instead of
quadrature amplitude modulation (QAM):

Having a baseband signal $x$, to send it over the channel, we modulate
(multiply) it by $\sqrt{2}\cos 2\pi f_{c}t$ where $f_{c}$ is the carrier
frequency and $t$ is the time index.

When sent over the channel, the transmitted signal $x$ undergoes perturbation
and thus the received signal becomes:

$z=\sqrt{2}\sum_{i}a_{i}x\cos(2\pi f_{c}t+\varphi_{i}(t))+w(t),$ (5)

where $a_{i}$ is the amplitude attenuation of path $i$ and
$\varphi_{i}(t)=2\pi f_{l}t+\varphi_{l}$ is the phase shift incurred by
Doppler frequency shift $f_{l}$ and/or any time delay $\varphi_{l}$. $w(t)$ is
an additive noise.

By developing the cosine term in $z$, we obtain

$\begin{split}z=&x\sqrt{2}\underbrace{\sum_{i}a_{i}\cos(\varphi_{i}(t))}_{\text{in-
phase component, }I(t)}\cos(2\pi f_{c}t)\\\
&-x\sqrt{2}\underbrace{\sum_{i}a_{i}\sin(\varphi_{i}(t))}_{\text{quadrature
component, }Q(t)}\sin(2\pi f_{c}t)+w(t),\end{split}$ (6)

From Section $2.4.2$ of [25]: According to the central limit theorem, if there
is a large number of channel paths, the in-phase and quadrature components of
the received signal, which are not correlated with each other, will exhibit
distributions that resemble the normal (Gaussian) distribution. Specifically,
each component will have an average value of zero and a variance of
$\Sigma/2$, which is equivalent to $\sigma^{2}$. The magnitude of the
perturbation $\sqrt{I(t)^{2}+Q(t)^{2}}$ thus becomes Rayleigh distributed.
This is the Rayleigh fading model. In addition, when the line-of-sight path is
large and has a known magnitude, the probabilistic model becomes a Rician
fading.

Furthermore, as $I(t)$ and $Q(t)$ are orthogonal, an equivalent complex
channel model $\hat{h}(t)=I(t)+jQ(t)=a(t)e^{-j\varphi(t)}$ can be derived.
Since the carrier frequency $f_{c}$ is not involved in $\hat{h}(t)$, this
representation is valid at baseband level. Thus, the complex channel model is
usually used to represent the received signal $\hat{h}(t)x+n(t)$ at baseband
with $\hat{h}(t)$ a complex entity.

Continuing from (6), to demodulate $z$ and obtain the baseband received signal
$y$, $z$ is first multiplied by $\sqrt{2}\cos 2\pi f_{c}t$ then the result
goes through a low pass filter.

$\begin{split}&z\sqrt{2}\cos(2\pi f_{c}t)\\\
=&2x\sum_{i}a_{i}\cos(\varphi_{i}(t))\cos^{2}(2\pi f_{c}t)\\\
&-2x\sum_{i}a_{i}\sin(\varphi_{i}(t))\sin(2\pi f_{c}t)\cos(2\pi f_{c}t)\\\
&+\sqrt{2}w(t)\cos(2\pi f_{c}t)\end{split}$
$\begin{split}=&x\sum_{i}a_{i}\cos(\varphi_{i}(t))(1+\cos(4\pi f_{c}t))\\\
&-x\sum_{i}a_{i}\sin(\varphi_{i}(t))\sin(4\pi f_{c}t)+\sqrt{2}w(t)\cos(2\pi
f_{c}t)\\\ \end{split}$ (7)

After the low pass filter, we obtain the received baseband signal

$\begin{split}y=&x\sum_{i}a_{i}\cos(\varphi_{i}(t))+n(t)\\\
=&\mathfrak{R}[\hat{h}(t)]x+n(t)\\\ =&h(t)x+n(t)\end{split}$ (8)

where $n(t)$ is the baseband equivalent noise with a zero-mean Gaussian
distribution and IID components (Section $2.2.4$ of [25]) and
$\mathfrak{R}[\hat{h}(t)]=I(t)$ is the real part of the channel.

As we are interested to send real values over the wireless channel in this
paper, one can easily see how equation (8) is valid to use with a real channel
$h=\mathfrak{R}[\hat{h}]$ following a Gaussian distribution with zero mean and
variance equal to $\sigma^{2}$.

##  Appendix B On the Estimated Gradient

###  B-A Proof of Lemma 1: Biased Estimator

Let $g_{k}$ have the form in (3), then the conditional expectation given
$\mathcal{H}_{k}$ can be written as

$\begin{split}&\mathbb{E}[g_{k}|\mathcal{H}_{k}]\\\
=&\mathbb{E}\Big{[}\Phi_{k}\sum_{i=1}^{N}\Big{(}h_{i,k+1}\tilde{f}_{i}\big{(}\theta_{k}^{\prime},S_{i,k+1}\big{)}+n_{i,k+1}\Big{)}\Big{|}\mathcal{H}_{k}\Big{]}\\\
\overset{(a)}{=}&\mathbb{E}\Big{[}\Phi_{k}\sum_{i=1}^{N}h_{i,k+1}\tilde{F}_{i}\big{(}\theta_{k}+\gamma_{k}\Phi_{k}\sum_{j=1}^{N}\big{(}\frac{h_{j,k}}{\sigma_{h}^{2}}+n_{j,k}\big{)}\big{)}\Big{|}\mathcal{H}_{k}\Big{]}\\\
\overset{(b)}{=}&\mathbb{E}\Big{[}\Phi_{k}\Big{(}\sum_{i=1}^{N}h_{i,k+1}\tilde{F}_{i}(\theta_{k})\\\
&+\gamma_{k}\sum_{i=1}^{N}h_{i,k+1}\sum_{j=1}^{N}\Big{(}\frac{h_{j,k}}{\sigma_{h}^{2}}+n_{j,k}\Big{)}\Phi_{k}^{T}\nabla\tilde{F}_{i}(\theta_{k})\\\
&+\gamma_{k}^{2}\sum_{i=1}^{N}h_{i,k+1}\Big{(}\sum_{j=1}^{N}\frac{h_{j,k}}{\sigma_{h}^{2}}+n_{j,k}\Big{)}^{2}\Phi_{k}^{T}\nabla^{2}\tilde{F}_{i}(\breve{\theta}_{k})\Phi_{k}\Big{)}\Big{|}\mathcal{H}_{k}\Big{]}\\\
\overset{(c)}{=}&\mathbb{E}\Big{[}\Phi_{k}\bigg{(}\frac{\gamma_{k}}{\sigma_{h}^{2}}\sum_{i=1}^{N}h_{i,k+1}h_{i,k}\Phi_{k}^{T}\nabla\tilde{F}_{i}(\theta_{k})+\\\
&+\gamma_{k}^{2}\sum_{i=1}^{N}h_{i,k+1}\Big{(}\sum_{j=1}^{N}\frac{h_{j,k}}{\sigma_{h}^{2}}+n_{j,k}\Big{)}^{2}\Phi_{k}^{T}\nabla^{2}\tilde{F}_{i}(\breve{\theta}_{k})\Phi_{k}\Big{)}\Big{|}\mathcal{H}_{k}\Big{]}\\\
=&\frac{\gamma_{k}}{\sigma_{h}^{2}}\sum_{i=1}^{N}\mathbb{E}\Big{[}h_{i,k+1}h_{i,k}\Big{|}\mathcal{H}_{k}\Big{]}\mathbb{E}\Big{[}\Phi_{k}\Phi_{k}^{T}\Big{|}\mathcal{H}_{k}\Big{]}\nabla\tilde{F}_{i}(\theta_{k})+\\\
&\gamma_{k}^{2}\sum_{i=1}^{N}\mathbb{E}\Big{[}h_{i,k+1}\Big{(}\sum_{j=1}^{N}\frac{h_{j,k}}{\sigma_{h}^{2}}+n_{j,k}\Big{)}^{2}\Phi_{k}\Phi_{k}^{T}\nabla^{2}\tilde{F}_{i}(\breve{\theta}_{k})\Phi_{k}\Big{|}\mathcal{H}_{k}\Big{]}\end{split}$
$\begin{split}\overset{(d)}{=}&\gamma_{k}\alpha_{2}\frac{K_{hh}}{\sigma_{h}^{4}}\sum_{i=1}^{N}\nabla
F_{i}(\theta_{k})+\\\
&\gamma_{k}^{2}\sum_{i=1}^{N}\mathbb{E}\Big{[}h_{i,k+1}\Big{(}\sum_{j=1}^{N}\frac{h_{j,k}}{\sigma_{h}^{2}}+n_{j,k}\Big{)}^{2}\Phi_{k}\Phi_{k}^{T}\nabla^{2}\tilde{F}_{i}(\breve{\theta}_{k})\Phi_{k}\Big{|}\mathcal{H}_{k}\Big{]}\\\
\overset{(e)}{=}&c_{1}\gamma_{k}(\nabla F(\theta_{k})+b_{k}),\end{split}$ (9)

where $(a)$ is by the definition in (2) and due to Assumption 1, $(b)$ is by
Taylor expansion and mean-value theorem and considering $\breve{\theta}_{k}$
between $\theta_{k}$ and
$\theta_{k}+\gamma_{k}\Phi_{k}\sum_{j=1}^{N}\big{(}\frac{h_{j,k}}{\sigma_{h}^{2}}+n_{j,k}\big{)}$.
$(c)$ is since $\mathbb{E}[h_{i,k+1}]=0$ for the first element and
$\mathbb{E}[h_{i,k+1}h_{j,k}]=0$ when $i\neq j$ and the independence of noise
for the second element. $(d)$ is due to Assumption 2. In $(e)$, we let
$c_{1}=\alpha_{2}\frac{K_{hh}}{\sigma_{h}^{4}}$ and the bias
$b_{k}=\gamma_{k}\frac{\sigma_{h}^{2}}{\alpha_{2}K_{hh}}\sum_{i=1}^{N}\mathbb{E}\Big{[}h_{i,k+1}\Big{(}\sum_{j=1}^{N}\frac{h_{j,k}}{\sigma_{h}^{2}}+n_{j,k}\Big{)}^{2}\Phi_{k}\Phi_{k}^{T}\nabla^{2}\tilde{F}_{i}(\breve{\theta}_{k})\Phi_{k}\Big{|}\mathcal{H}_{k}\Big{]}$.

###  B-B Proof of Lemma 2: Expected Norm Squared of the Estimated Gradient

To bound the norm squared of the one-point gradient estimate, we start by
letting
$\theta_{k}^{\prime}=\theta_{k}+\gamma_{k}\Phi_{k}\sum_{j=1}^{N}\big{(}\frac{h_{j,k}}{\sigma_{h}^{2}}+n_{j,k}\big{)}$

$\begin{split}&\mathbb{E}[\|g_{k}\|^{2}|\mathcal{H}_{k}]\\\
=&\mathbb{E}\Big{[}\|\Phi_{k}\|^{2}\Big{(}\sum_{i=1}^{N}h_{i,k+1}\tilde{f}_{i}\big{(}\theta_{k}^{\prime},S_{i,k+1}\big{)}+n_{i,k+1}\Big{)}^{2}\Big{|}\mathcal{H}_{k}\Big{]}\\\
\overset{(a)}{\leq}&\alpha_{3}^{2}N\sum_{i=1}^{N}\mathbb{E}\Big{[}\Big{(}h_{i,k+1}\tilde{f}_{i}\big{(}\theta_{k}^{\prime},S_{i,k+1}\big{)}+n_{i,k+1}\Big{)}^{2}\Big{|}\mathcal{H}_{k}\Big{]}\\\
\overset{(b)}{=}&\alpha_{3}^{2}N\sum_{i=1}^{N}\mathbb{E}\Big{[}h_{i,k+1}^{2}\tilde{f}_{i}^{2}\big{(}\theta_{k}^{\prime},S_{i,k+1}\big{)}+n_{i,k+1}^{2}\Big{|}\mathcal{H}_{k}\Big{]}\\\
\overset{(c)}{\leq}&\alpha_{3}^{2}N\sum_{i=1}^{N}\mathbb{E}\Big{[}\frac{h_{i,k+1}^{2}}{\sigma_{h}^{4}}\Big{(}\|f_{i}(0,S_{i,k+1})\|+L_{S_{i,k+1}}\|\theta_{k}^{\prime}\|\Big{)}^{2}\Big{|}\mathcal{H}_{k}\Big{]}\\\
&+N^{2}\alpha_{3}^{2}\sigma_{n}^{2}\\\
\leq&\alpha_{3}^{2}N\sum_{i=1}^{N}\mathbb{E}\Big{[}\frac{h_{i,k+1}^{2}}{\sigma_{h}^{4}}\Big{(}\|f_{i}(0,S_{i,k+1})\|+L_{S_{i,k+1}}\|\theta_{k}\|\\\
&\hskip
14.22636pt+L_{S_{i,k+1}}\gamma_{k}\|\Phi_{k}\|\Big{|}\sum_{j=1}^{N}\frac{h_{j,k}}{\sigma_{h}^{2}}+n_{j,k}\Big{|}\Big{)}^{2}\Big{|}\mathcal{H}_{k}\Big{]}+N^{2}\alpha_{3}^{2}\sigma_{n}^{2}\\\
\overset{(d)}{\leq}&3\alpha_{3}^{2}N\sum_{i=1}^{N}\mathbb{E}\Big{[}\frac{h_{i,k+1}^{2}}{\sigma_{h}^{4}}\Big{(}\|f_{i}(0,S_{i,k+1})\|^{2}+L_{S_{i,k+1}}^{2}\|\theta_{k}\|^{2}\\\
&\hskip
14.22636pt+\alpha_{3}^{2}L_{S_{i,k+1}}^{2}\gamma_{k}^{2}\big{(}\sum_{j=1}^{N}\frac{h_{j,k}}{\sigma_{h}^{2}}+n_{j,k}\big{)}^{2}\Big{)}\Big{|}\mathcal{H}_{k}\Big{]}+N^{2}\alpha_{3}^{2}\sigma_{n}^{2}\\\
\overset{(e)}{\leq}&\frac{3N^{2}\alpha_{3}^{2}}{\sigma_{h}^{2}}\Big{(}\mu_{S}+L_{S}\|\theta_{k}\|^{2}\Big{)}+3N\alpha_{3}^{4}L_{S}\gamma_{k}^{2}\sum_{i=1}^{N}\mathbb{E}\Big{[}\frac{h_{i,k+1}^{2}}{\sigma_{h}^{4}}\times\\\
&\big{(}\sum_{j\neq
i}\frac{h_{j,k}^{2}}{\sigma_{h}^{4}}+\frac{h_{i,k}^{2}}{\sigma_{h}^{4}}+n_{j,k}^{2}+2\sum_{j<l}\frac{h_{j,k}h_{l,k}}{\sigma_{h}^{4}}\big{)}\Big{|}\mathcal{H}_{k}\Big{]}+N^{2}\alpha_{3}^{2}\sigma_{n}^{2}\\\
\end{split}$
$\begin{split}\overset{(f)}{=}&\frac{3N^{2}\alpha_{3}^{2}}{\sigma_{h}^{2}}\Big{(}\mu_{S}+L_{S}\|\theta_{k}\|^{2}\Big{)}\\\
&+3N^{2}\alpha_{3}^{4}L_{S}\gamma_{k}^{2}\big{(}\frac{(N-1)\sigma_{h}^{4}}{\sigma_{h}^{8}}+\frac{\sigma_{h}^{4}+2K_{hh}^{2}}{\sigma_{h}^{8}}+\frac{N\sigma_{n}^{2}}{\sigma_{h}^{2}}+0\big{)}\\\
&+N^{2}\alpha_{3}^{2}\sigma_{n}^{2}\\\ :=&c_{2},\end{split}$ (10)

where $(a)$ is by Assumption 2 and Cauchy-Schwartz $(\sum_{j=1}^{N}1\cdot
a_{i})^{2}\leq(\sum_{j=1}^{N}a_{i}^{2})\cdot(\sum_{j=1}^{N}1^{2})=N(\sum_{j=1}^{N}a_{i}^{2})$,
for real $a_{i}$. $(b)$ is due the independence of noise, $(c)$ is by
Assumption 3, and $(d)$ is again by Cauchy-Schwartz. In $(e)$, we let
$\mu_{S}=\max_{i}\mathbb{E}[\|f_{i}(0,S_{i,k+1})\|^{2}|\mathcal{H}_{k}]$ and
$L_{S}=\max_{i}\mathbb{E}[L_{S_{i,k+1}}^{2}|\mathcal{H}_{k}]$. $(f)$ is due to
the _normally-distributed_ channel random variables,
$\mathbb{E}[h_{i,k+1}^{2}\sum_{j\neq i}h_{j,k}^{2}]=(N-1)\sigma_{h}^{4}$,
$\mathbb{E}[h_{i,k+1}^{2}h_{i,k}^{2}|\mathcal{H}_{k}]=\sigma_{h}^{4}+2K_{hh}^{2}$,
and $\mathbb{E}[h_{i,k+1}^{2}\sum_{j<l}h_{j,k}h_{l,k}]$ is equal to zero as
one of $h_{j,k}$ and $h_{l,k}$ will always be independent of the other terms
in the product and has a mean equal to zero.

###  B-C Proof of Lemma 4: Norm of the Bias

Considering the form of the bias in (9), by Assumptions 1, 2 and 4,

$\begin{split}&\|b_{k}\|\\\
\overset{(a)}{\leq}&\gamma_{k}\frac{\sigma_{h}^{2}}{\alpha_{2}K_{hh}}\sum_{i=1}^{N}\mathbb{E}\Big{[}2N\Big{|}h_{i,k+1}\sum_{j=1}^{N}\Big{(}\frac{h_{j,k}^{2}}{\sigma_{h}^{4}}+n_{j,k}^{2}\Big{)}\Big{|}\times\\\
&\hskip
81.09052pt\|\Phi_{k}\|\|\Phi_{k}^{T}\|\|\nabla^{2}\tilde{F}_{i}(\breve{\theta}_{k})\|\|\Phi_{k}\|\Big{|}\mathcal{H}_{k}\Big{]}\\\
\leq&2\gamma_{k}\frac{N\alpha_{1}\alpha_{3}^{3}\sigma_{h}^{2}}{\alpha_{2}K_{hh}}\sum_{i=1}^{N}\mathbb{E}\Big{[}\Big{|}h_{i,k+1}\sum_{j=1}^{N}\Big{(}\frac{h_{j,k}^{2}}{\sigma_{h}^{4}}+n_{j,k}^{2}\Big{)}\Big{|}\Big{|}\mathcal{H}_{k}\Big{]}\\\
\overset{(b)}{\leq}&2\gamma_{k}\frac{N\alpha_{1}\alpha_{3}^{3}\sigma_{h}^{2}}{\alpha_{2}K_{hh}}\sum_{i=1}^{N}\Big{[}\sigma_{h}\sqrt{\frac{2}{\pi}}\Big{(}2K_{hh}+\sqrt{\sigma_{h}^{4}-K_{hh}^{2}}\Big{)}\\\
&\hskip
79.6678pt+(N-1)\sigma_{h}^{3}\sqrt{\frac{2}{\pi}}+N\sqrt{\frac{2}{\pi}}\sigma_{h}\sigma_{n}^{2}\Big{]}\\\
=&2\gamma_{k}\frac{N^{2}\alpha_{1}\alpha_{3}^{3}\sigma_{h}^{3}}{\alpha_{2}K_{hh}}\sqrt{\frac{2}{\pi}}\times\\\
&\hskip
22.76228pt\Big{[}\Big{(}2K_{hh}+\sqrt{\sigma_{h}^{4}-K_{hh}^{2}}\Big{)}+(N-1)\sigma_{h}^{2}+N\sigma_{n}^{2}\Big{]}\\\
\overset{(c)}{=}&c_{3}\gamma_{k},\end{split}$ (11)

where $(a)$ is due to Jensen’s inequality and Cauchy-Schwartz, $(b)$ is by
using the half-normal distribution for normal random variables in absolute
value explained in the following paragraph, and in $(c)$,
$c_{3}=2\frac{N^{2}\alpha_{1}\alpha_{3}^{3}\sigma_{h}^{3}}{\alpha_{2}K_{hh}}\sqrt{\frac{2}{\pi}}\Big{[}\Big{(}2K_{hh}+\sqrt{\sigma_{h}^{4}-K_{hh}^{2}}\Big{)}+(N-1)\sigma_{h}^{2}+N\sigma_{n}^{2}\Big{]}$.

Let $X$ and $Y$ be two variables representing time-correlated channel
realizations at times $k$ and $k^{\prime}$ respectively. Assume they follow
the $\mathcal{N}(0,\sigma)$ distribution and they have a correlation
coefficient $\varrho$. Then, we can write $Y=\varrho X+\sqrt{1-\varrho^{2}}Z$,
where $Z$ is independent of X and following the same distribution
$\mathcal{N}(0,\sigma)$. Then,
$\mathbb{E}[|YX^{2}|]\leq\mathbb{E}[\varrho|X^{3}|+\sqrt{1-\varrho^{2}}|ZX^{2}|]=2\varrho\sqrt{\frac{2}{\pi}}\sigma^{3}+\sqrt{1-\varrho^{2}}\sqrt{\frac{2}{\pi}}\sigma\times\sigma^{2}=(2\varrho+\sqrt{1-\varrho^{2}})\sqrt{\frac{2}{\pi}}\sigma^{3}.$
If we substitute $\sigma=\sigma_{h}$ and
$\varrho=\frac{K_{hh}}{\sigma_{h}^{2}}$, we obtain the previous inequality
$(b)$.

##  Appendix C 1P-ZOFL algorithm convergence

###  C-A Stochastic Noise

To prove Lemma 5, we begin by demonstrating that the sequence
$\\{\sum_{k=K}^{K^{\prime}}\alpha_{k}e_{k}\\}_{K^{\prime}\geq K}$ is a
martingale. To do so, we have to prove that for all $K\textquoteright\geq K$,
$X_{K\textquoteright}=\sum_{k=K}^{K^{\prime}}\alpha_{k}e_{k}$ satisfies the
following two conditions:

(i)
$\mathbb{E}[X_{K\textquoteright+1}|X_{K\textquoteright}]=X_{K\textquoteright}$

(ii) $\mathbb{E}[\|X_{K\textquoteright}\|^{2}]<\infty$

We know that
$\mathbb{E}[e_{k}]=\mathbb{E}[g_{k}-\mathbb{E}[g_{k}|\mathcal{H}_{k}]]=\mathbb{E}_{\mathcal{H}_{k}}\Big{[}\mathbb{E}\Big{[}g_{k}-\mathbb{E}[g_{k}|\mathcal{H}_{k}]\Big{|}\mathcal{H}_{k}\Big{]}\Big{]}=0$
by the law of total expectation. Hence,
$\mathbb{E}[X_{K\textquoteright+1}|X_{K\textquoteright}]=\mathbb{E}\Big{[}\alpha_{K\textquoteright+1}e_{K\textquoteright+1}+\sum_{k=K}^{K^{\prime}}\alpha_{k}e_{k}\Big{|}\sum_{k=K}^{K^{\prime}}\alpha_{k}e_{k}\Big{]}=0+\sum_{k=K}^{K^{\prime}}\alpha_{k}e_{k}=X_{K\textquoteright}.$

In addition, $e_{k}$ and $e_{k^{\prime}}$ are uncorrelated for any $k\neq
k^{\prime}$ since (assuming $k>k^{\prime}$)
$\mathbb{E}\big{[}e_{k}^{T}e_{k\textquoteright}\big{]}=\mathbb{E}\big{[}\mathbb{E}[e_{k}^{T}e_{k\textquoteright}|\mathcal{H}_{k}]\big{]}=\mathbb{E}\big{[}e_{k\textquoteright}\mathbb{E}[e_{k}^{T}|\mathcal{H}_{k}]\big{]}=0$.
Thus,

$\begin{split}\mathbb{E}(\|\sum_{k=K}^{K^{\prime}}\alpha_{k}e_{k}\|^{2})&=\mathbb{E}(\sum_{k=K}^{K^{\prime}}\sum_{k^{\prime}=K}^{K^{\prime}}\alpha_{k}\alpha_{k^{\prime}}\langle
e_{k},e_{k^{\prime}}\rangle)\\\
&\overset{(a)}{=}\mathbb{E}(\sum_{k=K}^{K^{\prime}}\|\alpha_{k}e_{k}\|^{2})\\\
&\leq\sum_{k=K}^{\infty}\mathbb{E}(\alpha_{k}^{2}\|g_{k}-\mathbb{E}[g_{k}|\mathcal{H}_{k}]\|^{2})\\\
&=\sum_{k=K}^{\infty}\alpha_{k}^{2}\mathbb{E}(\|g_{k}\|^{2})-\mathbb{E}_{\mathcal{H}_{k}}(\|\mathbb{E}[g_{k}|\mathcal{H}_{k}]\|^{2})\\\
&\leq\sum_{k=K}^{\infty}\alpha_{k}^{2}\mathbb{E}(\|g_{k}\|^{2})\overset{(b)}{\leq}c_{2}\sum_{k=K}^{\infty}\alpha_{k}^{2}\gamma_{k}^{2}\overset{(c)}{<}\infty,\\\
\end{split}$ (12)

where $(a)$ is due to the uncorrelatedness $\mathbb{E}[\langle
e_{k},e_{k^{\prime}}\rangle]=0$, $(b)$ is by Lemma 2, and $(c)$ is by
Assumption 5. Therefore, both (i) and (ii) are satisfied and we can say that
$\\{\sum_{k=K}^{K^{\prime}}\alpha_{k}e_{k}\\}_{K^{\prime}\geq K}$ is a
martingale. This permits us to use Doob’s martingale inequality [38]:

For any constant $\nu>0$,

$\begin{split}\mathbb{P}(\sup_{K^{\prime}\geq
K}\|\sum_{k=K}^{K^{\prime}}\alpha_{k}e_{k}\|\geq\nu)&\leq\frac{1}{\nu^{2}}\mathbb{E}(\|\sum_{k=K}^{K^{\prime}}\alpha_{k}e_{k}\|^{2})\\\
&\overset{(a)}{\leq}\frac{c_{2}}{\nu^{2}}\sum_{k=K}^{\infty}\alpha_{k}^{2}\gamma_{k}^{2},\\\
\end{split}$ (13)

where $(a)$ is following the exact same steps as (12).

Since $c_{2}$ is a bounded constant and
$\lim_{K\rightarrow\infty}\sum_{k=K}^{\infty}\alpha_{k}^{2}\gamma_{k}^{2}=0$
by Assumption 5, we get
$\lim_{K\rightarrow\infty}\frac{c_{2}}{\nu^{2}}\sum_{k=K}^{\infty}\alpha_{k}^{2}\gamma_{k}^{2}=0$
for any bounded constant $\nu$. Hence, the probability that
$\|\sum_{k=K}^{K^{\prime}}\alpha_{k}e_{k}\|\geq\nu$ also vanishes as
$K\rightarrow\infty$, which concludes the proof.

###  C-B Proof of Theorem 1: Convergence Analysis

By the $L$-smoothness assumption and the algorithm update step
$\theta_{k+1}=\theta_{k}-\alpha_{k}g_{k}$, we have

$\begin{split}&F(\theta_{k+1})\\\ \leq&F(\theta_{k})-\alpha_{k}\langle\nabla
F(\theta_{k}),g_{k}\rangle+\frac{\alpha_{k}^{2}L}{2}\|g_{k}\|^{2}\\\
=&F(\theta_{k})-\alpha_{k}\langle\nabla
F(\theta_{k}),e_{k}+\mathbb{E}[g_{k}|\mathcal{H}_{k}]\rangle+\frac{\alpha_{k}^{2}L}{2}\|g_{k}\|^{2}\\\
=&F(\theta_{k})-\alpha_{k}\langle\nabla F(\theta_{k}),e_{k}\rangle-
c_{1}\alpha_{k}\gamma_{k}\|\nabla F(\theta_{k})\|^{2}\\\
&-c_{1}\alpha_{k}\gamma_{k}\langle\nabla
F(\theta_{k}),b_{k}\rangle+\frac{\alpha_{k}^{2}L}{2}\|g_{k}\|^{2}\\\
\overset{(a)}{\leq}&F(\theta_{k})-\alpha_{k}\langle\nabla
F(\theta_{k}),e_{k}\rangle-c_{1}\alpha_{k}\gamma_{k}\|\nabla
F(\theta_{k})\|^{2}\\\ &+\frac{c_{1}\alpha_{k}\gamma_{k}}{2}\|\nabla
F(\theta_{k})\|^{2}+\frac{c_{1}\alpha_{k}\gamma_{k}}{2}\|b_{k}\|^{2}+\frac{\alpha_{k}^{2}L}{2}\|g_{k}\|^{2}\\\
=&F(\theta_{k})-\alpha_{k}\langle\nabla
F(\theta_{k}),e_{k}\rangle-\frac{c_{1}\alpha_{k}\gamma_{k}}{2}\|\nabla
F(\theta_{k})\|^{2}\\\
&+\frac{c_{1}\alpha_{k}\gamma_{k}}{2}\|b_{k}\|^{2}+\frac{\alpha_{k}^{2}L}{2}\|g_{k}\|^{2}\\\
\end{split}$ (14)

where $(a)$ is by $-\langle
a,b\rangle\leq\frac{1}{2}\|a\|^{2}+\frac{1}{2}\|b\|^{2}$.

By taking the telescoping sum, we get

$\begin{split}F(\theta^{*})\leq
F(\theta_{K+1})\leq&F(\theta_{0})-\frac{c_{1}}{2}\sum_{k=0}^{K}\alpha_{k}\gamma_{k}\|\nabla
F(\theta_{k})\|^{2}\\\ &-\sum_{k=0}^{K}\alpha_{k}\langle\nabla
F(\theta_{k}),e_{k}\rangle\\\
&+\frac{c_{1}}{2}\sum_{k=0}^{K}\alpha_{k}\gamma_{k}\|b_{k}\|^{2}+\frac{c_{2}L}{2}\sum_{k=0}^{K}\alpha_{k}^{2}\|g_{k}\|^{2}\\\
\end{split}$ (15)

Hence,

$\begin{split}\sum_{k=0}^{K}\alpha_{k}\gamma_{k}\|\nabla
F(\theta_{k})\|^{2}\leq&\frac{2}{c_{1}}(F(\theta_{0})-F(\theta^{*}))\\\
&-\frac{2}{c_{1}}\sum_{k=0}^{K}\alpha_{k}\langle\nabla
F(\theta_{k}),e_{k}\rangle\\\
&+\sum_{k=0}^{K}\alpha_{k}\gamma_{k}\|b_{k}\|^{2}+\frac{c_{2}L}{c_{1}}\sum_{k=0}^{K}\alpha_{k}^{2}\|g_{k}\|^{2}\\\
\end{split}$ (16)

By Assumption 3, $\|\nabla F(\theta_{k})\|$ is bounded for any
$\theta_{k}\in\mathbb{R}^{d}$ and by taking the summation in (13) between $0$
and $\infty$, we have

$\lim_{K\rightarrow\infty}\|\sum_{k=0}^{K}\alpha_{k}\langle\nabla
F(\theta_{k}),e_{k}\rangle\|<\infty.$ (17)

From Lemma 4, we know that $\|b_{k}\|^{2}\sim\gamma_{k}^{2}$. Hence, by
Assumption 5,
$\lim_{K\rightarrow\infty}\sum_{k=0}^{K}\alpha_{k}\gamma_{k}^{3}<\infty$.

From Lemma 2 and by looking closely at the use of the Lipschitz continuity
property in (10), we can say $\|g_{k}\|^{2}\leq c$ for some $c>0$. Thus, again
by Assumption 5,
$\lim_{K\rightarrow\infty}\sum_{k=0}^{K}\alpha_{k}^{2}<\infty.$ We conclude
that

$\lim_{K\rightarrow\infty}\sum_{k=0}^{K}\alpha_{k}\gamma_{k}\|\nabla
F(\theta_{k})\|^{2}<\infty.$ (18)

Moreover, since the series $\sum_{k}\alpha_{k}\gamma_{k}$ diverges by
Assumption 5, we have $\lim_{k\rightarrow\infty}\inf\|\nabla
F(\theta_{k})\|=0.$

To prove that $\lim_{k\rightarrow\infty}\|\nabla F(\theta_{k})\|=0$, we
consider the hypothesis H) $\lim_{k\rightarrow\infty}\sup\|\nabla
F(\theta_{k})\|\geq\rho$ for an arbitrary $\rho>0$.

Assume (H) to be true. Then, we can always find an arbitrary subsequence
$\big{(}\|\nabla F(\theta_{k_{l}})\|\big{)}_{l\in\mathbb{N}}$ of $\|\nabla
F(\theta_{k})\|$, such that $\|\nabla
F(\theta_{k_{l}})\|\geq\rho-\varepsilon$, $\forall l$, for
$\rho-\varepsilon>0$ and $\varepsilon>0$.

Then, by the $L$-smoothness property and applying the descent step of the
algorithm,

$\begin{split}\|\nabla F(\theta_{k_{l}+1})\|&\geq\|\nabla
F(\theta_{k_{l}})\|-\|\nabla F(\theta_{k_{l}+1})-\nabla F(\theta_{k_{l}})\|\\\
&\geq\rho-\varepsilon-L\|\theta_{k_{l}+1}-\theta_{k_{l}}\|\\\
&=\rho-\varepsilon-L\alpha_{k_{l}}\|g_{k_{l}}\|\\\
&\geq\rho-\varepsilon-L\sqrt{c}\alpha_{k_{l}},\end{split}$ (19)

Since $k_{l}\rightarrow\infty$ as $l\rightarrow\infty$, we can always find a
subsequence of $(k_{l_{p}})_{p\in\mathbb{N}}$ such that
$k_{l_{p+1}}-k_{l_{p}}>1$. As $\alpha_{k_{l}}$ is vanishing, we consider
$(k_{l})_{l\in\mathbb{N}}$ starting from
$\alpha_{k_{l}}<\frac{\rho-\varepsilon}{L\sqrt{c}}$. Thus,

$\begin{split}&\sum_{k=0}^{\infty}\alpha_{k+1}\gamma_{k+1}\|\nabla
F(\theta_{k+1})\|^{2}\\\
\geq&(\rho-\varepsilon)^{2}\sum_{k=0}^{\infty}\alpha_{k+1}\gamma_{k+1}-2(\rho-\varepsilon)L\sqrt{c}\sum_{k=0}^{\infty}\alpha_{k+1}\gamma_{k+1}\alpha_{k}\\\
&+L^{2}c\sum_{k=0}^{\infty}\alpha_{k+1}\gamma_{k+1}\alpha_{k}^{2}\\\
=&+\infty,\end{split}$ (20)

as the first series diverges, and the second and the third converge by
Assumption 5. This implies that the series
$\sum_{k}\alpha_{k}\gamma_{k}\|\nabla F(\theta_{k})\|^{2}$ diverges. This is a
contradiction as this series converges almost surely by (18). Therefore,
hypothesis (H) cannot be true and $\|\nabla F(\theta_{k})\|$ converges to zero
almost surely.

###  C-C Proof of Theorem 2: Convergence Rate

Starting again from the $L$-smoothness in Lemma 3 and the algorithm update
step $\theta_{k+1}=\theta_{k}-\alpha_{k}g_{k}$, we have

$\begin{split}F(\theta_{k+1})&\leq F(\theta_{k})-\alpha_{k}\langle\nabla
F(\theta_{k}),g_{k}\rangle+\frac{\alpha_{k}^{2}L}{2}\|g_{k}\|^{2}.\end{split}$
(21)

Taking the conditional expectation given $\mathcal{H}_{k}$,

$\begin{split}&F(\theta_{k+1})\\\
\leq&F(\theta_{k})-c_{1}\alpha_{k}\gamma_{k}\langle\nabla F(\theta_{k}),\nabla
F(\theta_{k})+b_{k}\rangle+\frac{\alpha_{k}^{2}Lc_{2}}{2}\\\
=&F(\theta_{k})-c_{1}\alpha_{k}\gamma_{k}\|\nabla
F(\theta_{k})\|^{2}-c_{1}\alpha_{k}\gamma_{k}\langle\nabla
F(\theta_{k}),b_{k}\rangle\\\ &+\frac{\alpha_{k}^{2}Lc_{2}}{2}\\\ \end{split}$
$\begin{split}\overset{(a)}{\leq}&F(\theta_{k})-c_{1}\alpha_{k}\gamma_{k}\|\nabla
F(\theta_{k})\|^{2}+\frac{c_{1}\alpha_{k}\gamma_{k}}{2}\|\nabla
F(\theta_{k})\|^{2}\\\
&+\frac{c_{1}\alpha_{k}\gamma_{k}}{2}\|b_{k}\|^{2}+\frac{\alpha_{k}^{2}Lc_{2}}{2}\\\
=&F(\theta_{k})-\frac{c_{1}\alpha_{k}\gamma_{k}}{2}\|\nabla
F(\theta_{k})\|^{2}+\frac{c_{1}\alpha_{k}\gamma_{k}}{2}\|b_{k}\|^{2}+\frac{\alpha_{k}^{2}Lc_{2}}{2}\end{split}$
(22)

where $(a)$ is by $-\langle
a,b\rangle\leq\frac{1}{2}\|a\|^{2}+\frac{1}{2}\|b\|^{2}$

By considering a large value $K>0$ and taking the telescoping sum of (22), we
get

$\begin{split}&\mathbb{E}[F(\theta_{K+1})|\mathcal{H}_{K}]\\\
\leq&F(\theta_{0})-\frac{c_{1}}{2}\sum_{k=0}^{K}\alpha_{k}\gamma_{k}\|\nabla
F(\theta_{k})\|^{2}+\frac{c_{1}}{2}\sum_{k=0}^{K}\alpha_{k}\gamma_{k}\|b_{k}\|^{2}\\\
&+\frac{Lc_{2}}{2}\sum_{k=0}^{K}\alpha_{k}^{2}.\end{split}$ (23)

Given the assumption that
$F(\theta^{*})=\min_{\theta\in\mathbb{R}^{d}}F(\theta)$ exists, we know that
$\delta_{k}=F(\theta_{k})-F(\theta^{*})\geq 0$. Then,

$\begin{split}0\leq\mathbb{E}[\delta_{K+1}|\mathcal{H}_{K}]\leq&\delta_{0}-\frac{c_{1}}{2}\sum_{k=0}^{K}\alpha_{k}\gamma_{k}\|\nabla
F(\theta_{k})\|^{2}\\\
&+\frac{c_{1}c_{3}}{2}\sum_{k=0}^{K}\alpha_{k}\gamma_{k}^{3}+\frac{Lc_{2}}{2}\sum_{k=0}^{K}\alpha_{k}^{2}.\end{split}$
(24)

Finally,

$\begin{split}\sum_{k}\alpha_{k}\gamma_{k}\mathbb{E}[\|\nabla
F(\theta_{k})\|^{2}]&\leq\frac{2}{c_{1}c_{3}}\delta_{0}+\sum_{k}\alpha_{k}\gamma_{k}^{3}+\frac{Lc_{2}}{c_{1}}\sum_{k}\alpha_{k}^{2}\\\
\end{split}$ (25)

Let $\alpha_{k}$ and $\gamma_{k}$ have the forms given in Example 2. We know
that, $\forall K>0$,

$\begin{split}\sum_{k=0}^{K}\alpha_{k}\gamma_{k}^{3}&=\alpha_{0}\gamma_{0}^{3}+\sum_{k=1}^{K}\alpha_{k}\gamma_{k}^{3}\\\
&\leq\alpha_{0}\gamma_{0}^{3}\Big{(}1+\int_{0}^{K}(x+1)^{-\upsilon_{1}-3\upsilon_{2}}dx\Big{)}\\\
&\leq\alpha_{0}\gamma_{0}^{3}\Big{(}1+\frac{1}{\upsilon_{1}+3\upsilon_{2}-1}\Big{)}\\\
&=\alpha_{0}\gamma_{0}^{3}\Big{(}\frac{\upsilon_{1}+3\upsilon_{2}}{\upsilon_{1}+3\upsilon_{2}-1}\Big{)}.\\\
\end{split}$ (26)

Similarly,
$\sum_{k=0}^{K}\alpha_{k}^{2}\leq\alpha_{0}^{2}\big{(}\frac{2\upsilon_{1}}{2\upsilon_{1}-1}\big{)}$.
Next, when further $0<\upsilon_{1}+\upsilon_{2}<1$,

$\begin{split}\sum_{k=0}^{K}\alpha_{k}\gamma_{k}&\geq\alpha_{0}\gamma_{0}\int_{0}^{K+1}(x+1)^{-\upsilon_{1}-\upsilon_{2}}dx\\\
&=\frac{\alpha_{0}\gamma_{0}}{(1-\upsilon_{1}-\upsilon_{2})}\Big{(}(K+2)^{1-\upsilon_{1}-\upsilon_{2}}-1\Big{)}.\end{split}$
(27)

Thus, making use of inequality (25)

$\frac{\sum_{k}\alpha_{k}\gamma_{k}\mathbb{E}[\|\nabla
F(\theta_{k})\|^{2}]}{\sum_{k}\alpha_{k}\gamma_{k}}\leq\frac{A(1-\upsilon_{1}-\upsilon_{2})}{(K+2)^{1-\upsilon_{1}-\upsilon_{2}}-1}$
(28)

with
$A=\frac{2}{c_{1}\alpha_{0}\gamma_{0}}\delta_{0}+c_{3}^{2}\gamma_{0}^{2}\big{(}\frac{\upsilon_{1}+3\upsilon_{2}}{\upsilon_{1}+3\upsilon_{2}-1}\big{)}+\frac{Lc_{2}\alpha_{0}}{c_{1}\gamma_{0}}\big{(}\frac{2\upsilon_{1}}{2\upsilon_{1}-1}\big{)}$

##  Appendix D Non-Symmetrical Channels

Assuming non-symmetrical channels with $\mathbb{E}[h_{i,k}]=\mu_{h}$ and
$\sigma_{h}^{2}=\mathbb{E}[h_{i,k}^{2}]-\mu_{h}^{2}$, $\forall i,\forall k$,
the one-point gradient estimate becomes

$g_{k}=\Phi_{k}\sum_{i=1}^{N}\bigg{[}h_{i,k+1}\tilde{f}_{i}\big{(}\theta_{k}+\gamma_{k}\Phi_{k},S_{i,k+1}\big{)}+n_{i,k+1}\bigg{]}.$
(29)

And 1P-ZOFL is updated to Algorithm 2.

Algorithm 2 The 1P-ZOFL algorithm with non-symmetrical channels

Input: Initial model $\theta_{0}\in\mathbb{R}^{d}$, the initial step-sizes
$\alpha_{0}$ and $\gamma_{0}$, and the channels’ standard deviation
$\sigma_{h}$

1: for $k=0,2,4,...$ do

2: The server broadcasts $\theta_{k}+\gamma_{k}\Phi_{k}$

3: The server receives
$\sum_{i=1}^{N}h_{i,k+1}\tilde{f}_{i}\big{(}\theta_{k}+\gamma_{k}\Phi_{k},S_{i,k+1}\big{)}+n_{i,k+1}$

4: The server multiplies the received scalar sum by $\Phi_{k}$ to assemble
$g_{k}$ in (29)

5: The server updates $\theta_{k+1}=\theta_{k}-\alpha_{k}g_{k}$

6: end for

We then analyze the properties of our modified gradient estimate:

$\begin{split}&\mathbb{E}[g_{k}|\mathcal{H}_{k}]\\\
=&\mathbb{E}\Big{[}\Phi_{k}\sum_{i=1}^{N}\Big{(}h_{i,k+1}\tilde{f}_{i}\big{(}\theta_{k}+\gamma_{k}\Phi_{k},S_{i,k+1}\big{)}+n_{i,k+1}\Big{)}\Big{|}\mathcal{H}_{k}\Big{]}\\\
=&\mathbb{E}\Big{[}\Phi_{k}\sum_{i=1}^{N}\mu_{h}\tilde{F}_{i}\big{(}\theta_{k}+\gamma_{k}\Phi_{k}\big{)}\Big{|}\mathcal{H}_{k}\Big{]}\\\
=&\mu_{h}\gamma_{k}\sum_{i=1}^{N}\mathbb{E}\Big{[}\Phi_{k}\Phi_{k}^{T}\nabla\tilde{F}_{i}(\theta_{k})+\gamma_{k}\Phi_{k}\Phi_{k}^{T}\nabla^{2}\tilde{F}_{i}(\breve{\theta}_{k})\Phi_{k}\Big{|}\mathcal{H}_{k}\Big{]}\\\
=&c_{1}\gamma_{k}(\nabla F(\theta_{k})+b_{k}),\\\ \end{split}$ (30)

with $c_{1}=\frac{\mu_{h}\alpha_{2}}{\sigma_{h}^{2}}$ and
$b_{k}=\frac{\gamma_{k}}{\alpha_{2}}\mathbb{E}\Big{[}\Phi_{k}\Phi_{k}^{T}\nabla^{2}F_{i}(\breve{\theta}_{k})\Phi_{k}\Big{|}\mathcal{H}_{k}\Big{]}$.

Then, $\|b_{k}\|\leq\frac{\gamma_{k}}{\alpha_{2}}\alpha_{3}^{3}\alpha_{1}$,
now $c_{3}=\frac{\alpha_{3}^{3}\alpha_{1}}{\alpha_{2}}$ and $\|b_{k}\|\leq
c_{3}\gamma_{k}$.

Let $\theta_{k}^{\prime}=\theta_{k}+\gamma_{k}\Phi_{k}$,

$\begin{split}&\mathbb{E}[\|g_{k}\|^{2}|\mathcal{H}_{k}]\\\
=&\mathbb{E}\bigg{[}\|\Phi_{k}\|^{2}\Big{(}\sum_{i=1}^{N}h_{i,k+1}\tilde{f}_{i}\bigg{(}\theta_{k}^{\prime},S_{i,k+1}\bigg{)}+n_{i,k+1}\Big{)}^{2}\Big{|}\mathcal{H}_{k}\bigg{]}\\\
=&\alpha_{3}^{2}N\sum_{i=1}^{N}\mathbb{E}\bigg{[}h_{i,k+1}^{2}\tilde{f}_{i}^{2}\bigg{(}\theta_{k}^{\prime},S_{i,k+1}\bigg{)}+n_{i,k+1}^{2}\Big{|}\mathcal{H}_{k}\bigg{]}\\\
\leq&\alpha_{3}^{2}N\sum_{i=1}^{N}\mathbb{E}\bigg{[}\frac{h_{i,k+1}^{2}}{\sigma_{h}^{4}}\Big{(}\|f_{i}(0,S_{i,k+1})\|+L_{S_{i,k+1}}\|\theta_{k}^{\prime}\|\Big{)}^{2}\Big{|}\mathcal{H}_{k}\bigg{]}\\\
&+N^{2}\alpha_{3}^{2}\sigma_{n}^{2}\\\ \end{split}$
$\begin{split}\leq&3\alpha_{3}^{2}N\sum_{i=1}^{N}\mathbb{E}\bigg{[}\frac{h_{i,k+1}^{2}}{\sigma_{h}^{4}}\Big{(}\|f_{i}(0,S_{i,k+1})\|^{2}+L_{S_{i,k+1}}^{2}\|\theta_{k}\|^{2}\\\
&\hskip
85.35826pt+\alpha_{3}^{2}L_{S_{i,k+1}}^{2}\gamma_{k}^{2}\Big{)}\Big{|}\mathcal{H}_{k}\bigg{]}+N^{2}\alpha_{3}^{2}\sigma_{n}^{2}\\\
=&\frac{3N^{2}\alpha_{3}^{2}(\sigma_{h}^{2}+\mu_{h}^{2})}{\sigma_{h}^{2}}\Big{(}\mu_{S}+L_{S}\|\theta_{k}\|^{2}+\alpha_{3}^{2}L_{S}\gamma_{k}^{2}\Big{)}+N^{2}\alpha_{3}^{2}\sigma_{n}^{2}\\\
:=&c_{2}.\end{split}$ (31)

##  Appendix E Experimental results details

###  E-A Autoencoder

We compress images from the MNIST dataset using a lossy autoencoder with
encoder-decoder architecture. The encoder reduces the input size from 784 to
10 through three linear layers with dimensions 784-512-128-10. The first two
layers have ELU activation functions. The decoder reverses this process, going
from 10 to 784 through three linear layers (10-128-512-784), with ELU
activation in the first two layers and a sigmoid activation in the last layer
to ensure output pixel intensities between 0 and 1. Training is done on the
MNIST dataset for 10 epochs, using mean squared error loss and the Adam
optimizer.

###  E-B Performance of 1P-ZOFL vs SNR

An important remark is that high SNR is generally needed when we must decode
the information in the received signal. In our case, nothing is decoded; there
is no channel estimation or gradient extraction from the received signal.
Rather, the received signal is fed directly into the learning (the channel is
part of the learning). Ultimately, the amount of noise present in the system
does not affect our algorithms’ convergence: Examining (10) and (11), we see
that the noise variance $\sigma_{n}^{2}$ might only increase the norms of the
estimate and bias, but if we refer to (14), we find that both terms are
multiplied by step sizes, and we can counter the noise effect by decreasing
the step sizes’ constant parts (i.e, in the terms
$\frac{c_{1}\alpha_{k}\gamma_{k}}{2}\|b_{k}\|^{2}$ and
$\frac{\alpha_{k}^{2}L}{2}\|g_{k}^{(1P)}\|^{2}$). For the different plots in
Fig. 3, we decrease $\alpha_{0}$ and $\gamma_{0}$ when we increase
$\sigma_{n}^{2}$. However, this change of step sizes generally affects the
rate of convergence as they enter in the structure of $g_{k}$ and $b_{k}$
themselves and as can be seen in (4), smaller $\alpha_{0}$ and $\gamma_{0}$
increases the first term in the constant part of the upper bound on the
convergence rate. Thus, the effect shown on the plots in Fig. 3 is logical
given the changes of $\alpha_{0}$ and $\gamma_{0}$.