# Federated Learning over Wireless Networks:
A Band-limited Coordinated Descent Approach

Junshan Zhang1, Na Li2, and Mehmet Dedeoglu3 1,3School of Electrical, Computer
and Energy Engineering, Arizona State University
2School of Engineering and Applied Sciences, Harvard University
<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

We consider a many-to-one wireless architecture for federated learning at the
network edge, where multiple edge devices collaboratively train a model using
local data. The unreliable nature of wireless connectivity, together with
constraints in computing resources at edge devices, dictates that the local
updates at edge devices should be carefully crafted and compressed to match
the wireless communication resources available and should work in concert with
the receiver. Thus motivated, we propose SGD-based bandlimited coordinate
descent algorithms for such settings. Specifically, for the wireless edge
employing over-the-air computing, a common subset of k-coordinates of the
gradient updates across edge devices are selected by the receiver in each
iteration, and then transmitted simultaneously over k sub-carriers, each
experiencing time-varying channel conditions. We characterize the impact of
communication error and compression, in terms of the resulting gradient bias
and mean squared error, on the convergence of the proposed algorithms. We then
study learning-driven communication error minimization via joint optimization
of power allocation and learning rates. Our findings reveal that optimal power
allocation across different sub-carriers should take into account both the
gradient values and channel conditions, thus generalizing the widely used
water-filling policy. We also develop sub-optimal distributed solutions
amenable to implementation.

##  I Introduction

In many edge networks, mobile and IoT devices collecting a huge amount of data
are often connected to each other or a central node wirelessly. The unreliable
nature of wireless connectivity, together with constraints in computing
resources at edge devices, puts forth a significant challenge for the
computation, communication and coordination required to learn an accurate
model at the network edge. In this paper, we consider a many-to-one wireless
architecture for distributed learning at the network edge, where the edge
devices collaboratively train a machine learning model, using local data, in a
distributed manner. This departs from conventional approaches which rely
heavily on cloud computing to handle high complexity processing tasks, where
one significant challenge is to meet the stringent low latency requirement.
Further, due to privacy concerns, it is highly desirable to derive local
learning model updates without sending data to the cloud. In such distributed
learning scenarios, the communication between the edge devices and the server
can become a bottleneck, in addition to the other challenges in achieving edge
intelligence.

In this paper, we consider a wireless edge network with $M$ devices and an
edge server, where a high-dimensional machine learning model is trained using
distributed learning. In such a setting with unreliable and rate-limited
communications, local updates at sender devices should be carefully crafted
and compressed to make full use of the wireless communication resources
available and should work in concert with the receiver (edge server) so as to
learn an accurate model. Notably, lossy wireless communications for edge
intelligence presents unique challenges and opportunities [1], subject to
bandwidth and power requirements, on top of the employed multiple access
techniques. Since it often suffices to compute a function of the sum of the
local updates for training the model, over-the-air computing is a favorable
alternative to the standard multiple-access communications for edge learning.
More specifically, over-the-air computation [2, 3] takes advantage of the
superposition property of wireless multiple-access channel via simultaneous
analog transmissions of the local messages, and then computes a function of
the messages at the receiver, scaling signal-to-noise ratio (SNR) well with an
increasing number of users. In a nutshell, when multiple edge devices
collaboratively train a model, it is plausible to employ distributed learning
over-the-air.

We seek to answer the following key questions: 1) What is the impact of the
wireless communication bandwidth/power on the accuracy and convergence of the
edge learning? 2) What coordinates in local gradient signals should be
communicated by each edge device to the receiver? 3) How should the
coordination be carried out so that multiple sender devices can work in
concert with the receiver? 4) What is the optimal way for the receiver to
process the received noisy gradient signals to be used for the stochastic
gradient descent algorithm? 5) How should each sender device carry out power
allocation across subcarriers to transmit its local updates? Intuitively, it
is sensible to allocate more power to a coordinate with larger gradient value
to speed up the convergence. Further, power allocation should also be channel-
aware.

Figure 1: A bandlimited coordinate descent algorithm for distributed learning
over wireless multi-access channel

To answer the above questions, we consider an integrated learning and
communication scheme where multiple edge devices send their local gradient
updates over multi-carrier communications to the receiver for learning. Let
$K$ denote the number of subcarriers for communications, where $K$ is
determined by the wireless bandwidth. First, $K$ dimensions of the gradient
updates are determined (by the receiver) to be transmitted. Multiple methods
can be used for selecting $K$ coordinates, e.g., selecting the top-$k$ (in
absolute value) coordinates of the sum of the gradients or randomized uniform
selection. This paper will focus on randomly uniform selection (we elaborate
further on this in Section V). During the subsequent communications, the
gradient updates are transmitted only in the $K$-selected dimensions via over-
the-air computing over $K$ corresponding sub-carriers, each experiencing time-
varying channel conditions and hence time-varying transmission errors. The
devices are subject to power constraints, giving rise to a key question on how
to allocate transmission power across dimension, at each edge device, based on
the gradient update values and channel conditions. Thus, we explore joint
optimization of the power allocation and the learning rate to obtain the best
estimate of the gradient updates and minimize the impact of the communication
error. We investigate a centralized solution to this problem as a benchmark,
and then devise sub-optimal distributed solutions amenable to practical
implementation. We note that we have also studied the impact of errors of
synchronization across devices in this setting (we omit the details due to
limited space).

The main contributions of this paper are summarized as follows:

  * •

We take a holistic approach to study federated learning algorithms over
wireless MAC channels, and the proposed bandlimited coordinated descent(BLCD)
algorithm is built on innovative integration of computing in the air, multi-
carrier communications, and wireless resource allocation.

  * •

We characterize the impact of communication error and compression, in terms of
its resulting gradient bias and mean squared error (MSE), on the convergence
performance of the proposed algorithms. Specifically, when the communication
error is unbiased, the BLCD algorithm would converge to a stationary point
under very mild conditions on the loss function. In the case the bias in the
communication error does exist, the iterates of the BLCD algorithm would
return to a contraction region centered around a scaled version of the bias
infinitely often.

  * •

To minimize the impact of the communication error, we study joint optimization
of power allocation at individual devices and learning rates at the receiver.
Observe that since there exists tradeoffs between bias and variance,
minimizing the MSE of the communication error does not necessarily amount to
minimizing the bias therein. Our findings reveal that optimal power allocation
across different sub-carriers should take into account both the gradient
values and channel conditions, thus generalizing the widely used water-filling
policy. We also develop sub-optimal distributed solutions amenable to
implementation. In particular, due to the power constraints at individual
devices, it is not always feasible to achieve unbiased estimators of the
gradient signal across the coordinates. To address this complication, we
develop a distributed algorithm which can drive the bias in the communication
error to (close to) zero under given power constraints and then reduce the
corresponding variance as much as possible.

##  II Related Work

Communication-efficient SGD algorithms are of great interest to reduce latency
caused by the transmission of the high dimensional gradient updates with
minimal performance loss. Such algorithms in the ML literature are based on
compression via quantization [4, 5, 6, 7], sparsification [8, 9, 10] and
federated learning [11] (or local updates [12]), where lossless communication
is assumed to be provided. At the wireless edge, physical-layer design and
communication loss should be taken into consideration for the adoption of the
communication-efficient algorithms.

Power allocation for over-the-air computation is investigated for different
scenarios in many other works [13, 14, 15, 16, 17] including MIMO, reduced
dimensional MIMO, standard many to one channel and different channel models.
In related works on ML over wireless channels, [18, 19, 20, 21, 22, 23, 24,
25] consider over-the-air transmissions for training of the ML model. The
authors in [21] propose sparsification of the updates with compressive sensing
for further bandwidth reduction, and recovered sum of the compressed sparse
gradients is used for the update. They also apply a similar framework for
federated learning and fading channels in [22]. [18] considers a broadband
aggregation for federated learning with opportunistic scheduling based on the
channel coefficients for a set of devices uniformly distributed over a ring.
Lastly, [25] optimize the gradient descent based learning over multiple access
fading channels. It is worth noting that the existing approaches for
distributed learning in wireless networks do not fully account for the
characteristics of lossy wireless channels. It is our hope that the proposed
BLCD algorithms can lead to an innovative architecture of distributed edge
learning over wireless networks that accounts for computation, power, spectrum
constraints and packet losses.

##  III Federated Learning over Wireless Multi-access Networks

###  III-A Distributed Edge Learning Model

Consider an edge computing environment with $M$ devices
$\mathcal{M}=\\{1,\ldots,M\\}$ and an edge server. As illustrated in Figure 1,
a high-dimensional ML model is trained at the server by using an SGD based
algorithm, where stochastic gradients are calculated at the devices with the
data points obtained by the devices and a (common) subset of the gradient
updates are transmitted through different subcarriers via over-the-air.

The general edge learning problem is as follows:

$\min_{w\in\mathbb{R}^{d}}f(w):=\frac{1}{M}\sum_{m=1}^{M}\mathbb{E}_{\xi_{m}}[l(w,\xi_{m})],$
(1)

in which $l(\cdot)$ is the loss function, and edge device $m$ has access to
inputs $\xi_{m}$. Such optimization is typically performed through empirical
risk minimization iteratively. In the sequel, we let $w_{t}$ denote the
parameter value of the ML model at communication round $t$, and at round $t$
edge device $m$ uses its local data $\xi_{m,t}$ to compute a stochastic
gradient $g^{m}_{t}(w_{t}):=\nabla l(w_{t},\xi_{m,t})$. Define
$g_{t}(w_{t})=\frac{1}{M}\sum_{m=1}^{M}g^{m}_{t}(w_{t})$. The standard vanilla
SGD algorithms is given as

$w_{t+1}=w_{t}-\gamma g_{t}(w_{t})$ (2)

with $\gamma$ being the learning rate. Nevertheless, different updates can be
employed for different SGD algorithms, and this study will focus on
communication-error-aware SGD algorithms.

###  III-B Bandlimited Coordinate Descent Algorithm

Due to the significant discrepancy between the wireless bandwidth constraint
and the high-dimensional nature of the gradient signals, we propose a sparse
variant of the SGD algorithm over wireless multiple-access channel, named as
bandlimited coordinate descent (BLCD), in which at each iteration only a
common set of $K$ coordinates, $I(t)\subset\\{1,\ldots,d\\}$ (with $K\ll d$),
of the gradients are selected to be transmitted through over-the-air computing
for the gradient updates. The details of coordinate selection for the BLCD
algorithm are relegated to Section VI. Worth noting is that due to the
unreliable nature of wireless connectivity, the communication is assumed to be
lossy, resulting in erroneous estimation of the updates at the receiver.
Moreover, gradient correction is performed by keeping the difference between
the update made at the receiver and the gradient value at the transmitter for
the subsequent rounds, as gradient correction dramatically improves the
convergence rate with sparse gradient updates [9].

For convenience, we first define the gradient sparsification operator as
follows.

######  Definition 1.

$C_{I}:\mathbb{R}^{d}\rightarrow\mathbb{R}^{d}$ for a set
$I\subseteq\\{1,\ldots,d\\}$ as follows: for every input $x\in\mathbb{R}^{d}$,
$\big{(}C_{I}(x)\big{)}_{j}$ is $(x)_{j}$ for $j\in I$ and $0$ otherwise.

Since this operator $C_{I}$ compress a $d$-dimensional vector to a
$k$-dimension one, we will also refer this operator as compression operator in
the rest of the paper.

Algorithm 1 Bandlimited Coordinate Descent Algorithm

1: Input: Sample batches $\xi_{m,t}$, model parameters $w_{1}$, initial
learning rate $\gamma$, sparsification operator $C_{t}(.)$, $\forall
m=1,\dots,M;\forall t=1,\dots,T.$

2: Initialize: $r_{t}^{m}:=0$.

3: for $t=1:T$ do

4: for $m=1:M$ do

5: $g_{t}^{m}(w_{t}):=\text{stochasticGradient}(f(w_{t},\xi_{m,t}))$

6: $u_{t}^{m}:=\gamma g_{t}^{m}(w_{t})+r_{t}^{m}$

7: $r_{t+1}^{m}:=u_{t}^{m}-C_{t}(u_{t}^{m})$

8: Compute power allocation coefficients $b_{km}^{*},\forall k=1,\dots,K$.

9: Transmit $\mathbf{b}^{*}\odot C_{t}(u_{t}^{m})$

10: end for

11: Compute gradient estimator $\hat{G}_{t}(w_{t})$

12: $w_{t+1}:=w_{t}-\hat{G}_{t}(w_{t})$.

13: Broadcast $w_{t+1}$ back to all transmitters.

14: end for

With a bit abuse of notation, we let $C_{t}$ denote $C_{I(t)}$ for convenience
in the following. Following [26], we incorporate the sparsification error made
in each iteration (by the compression operator $C_{t}$) into the next step to
alleviate the possible gradient bias therein and improve the convergence
possible. Specifically, as in [26], one plausible way for compression error
correction is to update the gradient correction term as follows:

$\displaystyle r_{t+1}^{m}$ $\displaystyle=u_{t}^{m}-C_{t}(u_{t}^{m}),$ (3)
$\displaystyle u_{t}^{m}$ $\displaystyle\triangleq\gamma
g^{m}_{t}(w_{t})+r_{t}^{m}$ (4)

which $r_{t+1}^{m}$ keeps the error in the sparsification operator that is in
the memory of user $m$ at around $t$, and $u_{t}^{m}$ is the scaled gradient
with correction at device $m$ where the scaling factor $\gamma$ is the
learning rate in equation (2). (We refer readers to [26] for more insights of
this error-feedback based compression SGD.) Due to the lossy nature of
wireless communications, there would be communication errors and the gradient
estimators at the receiver would be erroneous. In particular, the gradient
estimator at the receiver in the BLCD can be written as

$\hat{G}_{t}(w_{t})=\frac{1}{M}\sum_{m=1}^{M}C_{t}\left(u_{t}^{m}\right)+\epsilon_{t},$
(5)

where $\epsilon_{t}$ denotes the random communication error in round $t$. In a
nutshell, the bandlimited coordinate descent algorithm is outlined in
Algorithm 1.

Recall that $g_{t}(w_{t})=\frac{1}{M}\sum_{m=1}^{M}g^{m}_{t}(w_{t})$ and
define $r_{t}\triangleq\frac{1}{M}\sum_{m=1}^{M}r^{m}_{t}$. Thanks to the
common sparsification operator across devices, the update in the SGD algorithm
at communicatioon round $t$ is given by

$w_{t+1}=w_{t}-\big{[}C_{t}(\gamma g_{t}(w_{t})+r_{t})+\epsilon_{t}\big{]}.$
(6)

To quantify the impact of the communication error, we use the corresponding
communication-error free counterpart as the benchmark, defined as follows:

$\hat{w}_{t+1}=w_{t}-C_{t}(\gamma g_{t}(w_{t})+r_{t}).$ (7)

It is clear that $w_{t+1}=\hat{w}_{t+1}-\epsilon_{t}$.

For convenience, we define $\tilde{w}_{t}\triangleq{w}_{t}-r_{t}$. It can be
shown that $\tilde{w}_{t+1}=\tilde{w}_{t}-\gamma g_{t}(w_{t})-\epsilon_{t}$.
Intuitively, $w_{t+1}$ in (6) is a noisy version of the iterate
$\hat{w}_{t+1}$ in (7), which implies that $\tilde{w}_{t+1}$ is a noisy
version of the compression-error correction of $\hat{w}_{t+1}$ in (7), where
the “noisy perturbation” is incurred by the communication error.

###  III-C BLCD Coordinate Transmissions over Multi-Access Channel

Figure 2: A multi-access communication protocol for bandlimited coordinate
selection and transmission.

A key step in the BLCD algorithm is to achieve coordinate synchronization of
the transmissions among many edge devices. To this end, we introduce a
receiver-driven low-complexity multi-access communication protocol, as
illustrated in Fig. 2, with the function $C_{t}(x)$ denoting the compression
of $x$ at round $t$. Let $I(t)$ (of size $K$) denote the subset of coordinates
chosen for transmission by the receiver at round $t$. Observe that the updates
at the receiver are carried out only in the dimensions $I(t)$. Further, the
edge receiver can broadcast its updated iterate to participant devices, over
the reverse link. This task is quite simple, given the broadcast nature of
wireless channels. In the transmissions, each coordinate of the gradient
updates is mapped to a specific subcarrier and then transmitted through the
wireless MAC channel, and the coordinates transmitted by different devices
over the same subcarrier are received by the edge server in the form of an
aggregate sum. It is worth noting that the above protocol is also applicable
to the case when the SGD updates are carried out for multiple rounds at the
devices.

When there are many edge devices, over-the-air computation can be used to take
advantage of superposition property of wireless multiple-access channel via
simultaneous analog transmissions of the local updates. More specifically, at
round t, the received signal in subcarrier $k$ is given by:

$y_{k}(t)=\sum_{m=1}^{M}b_{km}(t)h_{km}(t)x_{km}(t)+n_{k}(t)$ (8)

where $b_{km}(t)$ is a power scaling factor, $h_{km}(t)$ is the channel gain,
and $x_{km}(t)$ is the message of user $m$ through the subcarrier $k$,
respectively, and $n_{k}(t)\sim\mathcal{N}(0,\sigma^{2})$ is the channel
noise.

To simplify notation, we omit $(t)$ when it is clear from the context in the
following. Specifically, the message $x_{km}=(C_{t}(u^{m}_{t}))_{l(k)}$, with
a one-to-one mapping $l(k)=(I(t))_{k}$, which indicates the $k$-th element of
$I(t)$, transmitted through the $k$-th subcarrier. The total power that a
device can use in the transmission is limited in practical systems. Without
loss of generality, we assume that there is a power constraint at each device,
given by $\sum_{k=1}^{K}\left\lvert b_{km}x_{km}\right\rvert^{2}\leq E_{m},\
\forall m\in\\{1,\ldots,M\\}$. Note that $b_{km}$ hinges heavily upon both
$\bm{h}_{m}=[h_{1m},\ldots,h_{Km}]^{\top}$ and
$\bm{x}_{m}=[x_{1m},\ldots,x_{Km}]^{\top}$, and a key next step is to optimize
$b_{km}(\bm{h}_{m},\bm{x}_{m})$. In each round, each device optimizes its
power allocation for transmitting the selected coordinates of its update
signal over the $K$ subcarriers, aiming to minimize the communication error so
as to achieve a good estimation of $G_{t}(w_{t})$ (or its scaled version) for
the gradient update, where

$G_{t}(w_{t})\triangleq\frac{1}{M}\sum_{m=1}^{M}C_{t}(u_{t}^{m}).$

From the learning perspective, based on $\\{y_{k}\\}_{k=1}^{K}$, it is of
paramount importance for the receiver to get a good estimate of
$G_{t}(w_{t})$. Since $n_{k}(t)$ is Gaussian noise, the optimal estimator is
in the form of

$\big{(}\widehat{G}_{t}(w_{t})\big{)}_{k}=\begin{cases}\alpha_{l(k)}y_{l(k)},&k\in
I(t)\\\ 0&\text{otherwise}\end{cases}\vspace{0.05in}$ (9)

where $\\{\alpha_{k}\\}_{k=1}^{K}$ are gradient estimator coefficients for
subcarriers. It follows that the communication error (i.e., the gradient
estimation error incurred by lossy communications) is given by

$\epsilon_{t}=\widehat{G}_{t}(w_{t})-G_{t}(w_{t}).$ (10)

We note that $\\{\alpha_{k}\\}_{k=1}^{K}$ are intimately related to the
learning rates for the $K$ coordinates, scaling the learning rate to be
$\\{\gamma\alpha_{k}\\}_{k=1}^{K}$. It is interesting to observe that the
learning rates in the proposed BLCD algorithm are essentially different across
the dimensions, due to the unreliable and dynamically changing channel
conditions across different subcarriers.

##  IV Impact of Communication Error and Compression on BLCD Algorithm

Recall that due to the common sparsification operator across devices, the
update in the SGD algorithm at communication round $t$ is given by

$w_{t+1}=w_{t}-\big{[}C_{t}(\gamma g_{t}(w_{t})+r_{t})+\epsilon_{t}\big{]}.$

Needless to say, the compression operator $C_{t}$ plays a critical role in
sparse transmissions. In this study, we impose the following standard
assumption on the compression rate of the operator.

######  Assumption 1.

For a set of the random compression operators $\\{C_{t}\\}_{t=1}^{T}$ and any
$x\in\mathbb{R}^{d}$, it holds

$\operatorname{\mathbb{E}}\left\lVert
x-C_{t}(x)\right\rVert^{2}\leq(1-\delta)\left\lVert x\right\rVert^{2}$ (11)

for some $\delta\in(0,1]$.

We impose the following standard assumptions on the non-convex objective
function $f(\cdot)$ and the corresponding stochastic gradients
$g^{m}_{t}(w_{t})$ computed with the data samples of device $m$ in round $t$.
(We assume that the data samples $\\{\xi_{m,t}\\}$ are i.i.d. across the
devices and time.)

######  Assumption 2.

(Smoothness) A function $f:\mathbb{R}^{d}\rightarrow\mathbb{R}$ is L-smooth if
for all ${x},{y}\in\mathbb{R}^{d}$, it holds

$\lvert f({y})-f({x})-\langle\nabla
f({x}),{y}-{x}\rangle\rvert\leq\frac{L}{2}\left\lVert{y}-{x}\right\rVert^{2}.$
(12)

######  Assumption 3.

For any $x\in\mathbb{R}^{d}$ and for any $m=1,\ldots,M$, a stochastic gradient
$g_{t}^{m}(x),\forall t$, satisfies

$\operatorname{\mathbb{E}}[g_{t}^{m}(x)]=\nabla f(x),\textrm{
}\operatorname{\mathbb{E}}\left\lVert g_{t}^{m}(x)\right\rVert^{2}\leq G^{2}$
(13)

where $G>0$ is a constant.

$\displaystyle\mathbb{E}_{t}[f(\tilde{w}_{t+1})]\leq$ $\displaystyle
f(\tilde{w}_{t})+\langle\nabla
f(\tilde{w}_{t}),\mathbb{E}_{t}[\tilde{w}_{t+1}-\tilde{w}_{t}]\rangle+\frac{L}{2}\mathbb{E}_{t}[\lVert\tilde{w}_{t+1}-\tilde{w}_{t}\rVert^{2}]$
$\displaystyle=f(\tilde{w}_{t})-\langle\nabla
f(\tilde{w}_{t}),\gamma\mathbb{E}_{t}[g_{t}(w_{t})]+\mathbb{E}_{t}[\epsilon_{t}]\rangle+\frac{L}{2}\mathbb{E}_{t}[\lVert\gamma
g_{t}(w_{t})\rVert^{2}]+\frac{L}{2}\mathbb{E}_{t}[\lVert\epsilon_{t}\rVert^{2}]+L\mathbb{E}_{t}[\langle\gamma
g_{t}(w_{t}),\epsilon_{t}\rangle]$
$\displaystyle=f(\tilde{w}_{t})-\langle\nabla
f({w}_{t}),\gamma\mathbb{E}_{t}[g_{t}(w_{t})]+\mathbb{E}_{t}[\epsilon_{t}]\rangle-\langle\nabla
f(\tilde{w}_{t})-\nabla
f({w}_{t}),\gamma\mathbb{E}_{t}[g_{t}(w_{t})]+\mathbb{E}_{t}[\epsilon_{t}]\rangle+\frac{L}{2}\mathbb{E}_{t}[\|\epsilon_{t}\rVert_{2}^{2}]+L\mathbb{E}_{t}[\langle\gamma
g_{t}(w_{t}),\epsilon_{t}\rangle]+\frac{L}{2}\mathbb{E}_{t}[\|\gamma
g_{t}(w_{t})\rVert^{2}$ $\displaystyle\leq f(\tilde{w}_{t})-\gamma\lVert\nabla
f(w_{t})\rVert_{2}^{2}-\langle\nabla
f(w_{t}),\mathbb{E}_{t}[\epsilon_{t}]\rangle+\frac{\rho}{2}\lVert\gamma\nabla
f(w_{t})+\mathbb{E}_{t}[\epsilon_{t}]\rVert_{2}^{2}+\frac{L^{2}}{2\rho}\mathbb{E}_{t}[\lVert
r_{t}\rVert_{2}^{2}]+\frac{L}{2}\mathbb{E}_{t}[\lVert\epsilon_{t}\rVert_{2}^{2}]+L\langle\nabla
f(w_{t}),\mathbb{E}_{t}[\epsilon_{t}]\rangle+\frac{L\gamma^{2}}{2}\mathbb{E}_{t}\lVert
g_{t}(w_{t})\rVert_{2}^{2}$ $\displaystyle\leq
f(\tilde{w}_{t})-\gamma\lVert\nabla f(w_{t})\rVert_{2}^{2}+(L-1)\lVert\nabla
f(w_{t})\rVert\lVert\mathbb{E}_{t}[\epsilon_{t}]\rVert+\frac{\rho}{2}\left(\gamma^{2}\lVert\nabla
f(w_{t})\rVert_{2}^{2}+\lVert\mathbb{E}_{t}[\epsilon_{t}]\rVert_{2}^{2}+2\gamma\langle\nabla
f(w_{t}),\mathbb{E}_{t}[\epsilon_{t}]\rangle\right)+\frac{L^{2}}{2\rho}\mathbb{E}_{t}[\lVert
r_{t}\rVert_{2}^{2}]+\frac{L}{2}\mathbb{E}_{t}[\lVert\epsilon_{t}\rVert_{2}^{2}]+\frac{L\gamma^{2}}{2}G^{2}$
$\displaystyle\leq f(\tilde{w}_{t})-\gamma\lVert\nabla
f(w_{t})\rVert_{2}^{2}+(L-1+2\gamma)\lVert\nabla
f(w_{t})\rVert\lVert\mathbb{E}_{t}[\epsilon_{t}]\rVert+\frac{\gamma^{2}\rho}{2}\lVert\nabla
f(w_{t})\rVert_{2}^{2}+\frac{L^{2}}{2\rho}\mathbb{E}_{t}[\lVert
r_{t}\rVert_{2}^{2}]+\lVert\mathbb{E}_{t}[\epsilon_{t}]\rVert_{2}^{2}+\frac{L}{2}\mathbb{E}_{t}[\lVert\epsilon_{t}\rVert_{2}^{2}]+\frac{L\gamma^{2}}{2}G^{2}$
$\displaystyle=f(\tilde{w}_{t})-\gamma\left[1-\frac{\rho}{2}\gamma\right]\lVert\nabla
f(w_{t})\rVert_{2}^{2}+(L-1+2\gamma)\lVert\nabla
f(w_{t})\rVert\lVert\mathbb{E}_{t}[\epsilon_{t}]\rVert+\frac{L^{2}}{2\rho}\mathbb{E}_{t}[\lVert
r_{t}\rVert_{2}^{2}]+\lVert\mathbb{E}_{t}[\epsilon_{t}]\rVert_{2}^{2}+\frac{L}{2}\mathbb{E}_{t}[\lVert\epsilon_{t}\rVert_{2}^{2}]+\frac{L\gamma^{2}}{2}G^{2}$
(14)

It follows directly from [26] that $\mathbb{E}[\lVert
r_{t}\rVert_{2}^{2}]\leq\frac{4(1-\delta)}{\delta^{2}}\gamma^{2}G^{2}.$ Recall
that $\tilde{w}_{t+1}=\tilde{w}_{t}-\gamma g_{t}(w_{t})-\epsilon_{t}$ and that
$\tilde{w}_{t+1}$ can be viewed as a noisy version of the compression-error
correction of $\hat{w}_{t+1}$ in (7), where the “noisy perturbation” is
incurred by the communication error. For convenience, let
$\mathbb{E}_{t}[\epsilon_{t}]$ denote the gradient bias incurred by the
communication error and $\mathbb{E}_{t}[\lVert\epsilon_{t}\rVert_{2}^{2}]$ be
the corresponding mean square error, where $\mathbb{E}_{t}$ is taken with
respect to channel noise.

Let $\eta=\frac{L-1+2\gamma}{\gamma(2-\rho\gamma)}$ with $0<\rho<2$. Let
$f^{*}$ denote the globally minimum value of $f$. We have the following main
result on the iterates in the BLCD algorithm.

######  Theorem 1.

Under Assumptions 1, 2 and 3, the iterates $\\{w_{t}\\}$ in the BLCD algorithm
satisfies that

$\displaystyle\frac{1}{T+1}\sum_{t=0}^{T}\left(\lVert\nabla
f(w_{t})\rVert_{2}-\eta\lVert\underbrace{\mathbb{E}_{t}[\epsilon_{t}]}_{\mbox{bias}}\rVert_{2}\right)^{2}$
$\displaystyle\hskip
7.22743pt\leq\frac{1}{T+1}\sum_{t=0}^{T}\left[\frac{L\eta}{L-1+2\gamma}\underbrace{\mathbb{E}_{t}[\lVert\epsilon_{t}\rVert_{2}^{2}]}_{\mbox{MSE}}+\left(1+\eta^{2}\right)\lVert\underbrace{\mathbb{E}_{t}[\epsilon_{t}]}_{\mbox{bias}}\rVert_{2}^{2}\right]$
$\displaystyle\hskip
7.22743pt+\frac{2}{T+1}\frac{f(w_{0})-f^{*}}{\gamma(2-\rho\gamma)}+\left(\frac{L}{\rho}\frac{2(1-\delta)}{\delta^{2}}+\frac{1}{2}\right)\frac{2L\gamma
G^{2}}{2-\rho\gamma}.$ (15)

###### Proof.

Due to the limited space, we outline only a few main steps for the proof.
Recall that $\tilde{w}_{t}={w}_{t}-r_{t}$. It can be shown that
$\tilde{w}_{t+1}=\tilde{w}_{t}-\gamma g_{t}(w_{t})-\epsilon_{t}$. As shown in
(14), using the properties of the iterates in the BLCD algorithm and the
smoothness of the objective function $f$, we can establish an upper bound on
$\mathbb{E}_{t}[f(\tilde{w}_{t+1})]$ in terms of $f(\tilde{w}_{t})$ the
corresponding gradient $\nabla f(w_{t})$, and the gradient bias and MSE due to
the communication error. Then, (15) can be obtained after some further
algebraic manipulation. ∎

Remarks. Based on Theorem 1, we have a few observations in order.

  * •

We first examine the four terms on the right hand side of (15): The first two
terms capture the impact on the gradient by the time average of the bias in
the communication error $\epsilon_{t}$ and that of the corresponding the mean
square, denoted as MSE; the two items would go to zero if the bias and the MSE
diminish; the third term is a scaled version of $f(w_{0})-f^{*}$ and would go
to zero as long as $\gamma=O(T^{-\beta})$ with $\beta<1$; and the fourth term
is proportional to $\gamma$ and would go to zero when $\gamma\rightarrow 0$.

  * •

If the right hand side of (15) diminishes as $T\rightarrow\infty$, the
iterates in the BLCD algorithm would “converge” to a neighborhood around
$\eta\lVert\mathbb{E}_{t}[\epsilon_{t}]\rVert_{2}$, which is a scaled version
of the bias in the communication error. For convenience, let
$\bar{\epsilon}=\limsup_{t}\lVert\mathbb{E}_{t}[\epsilon_{t}]\rVert_{2}$, and
define a contraction region as follows:

$A_{\gamma}=\left\\{w_{t}:\lVert\nabla
f(w_{t})\rVert_{2}\leq(\eta+\Delta)\bar{\epsilon}\right\\}.$

where $\Delta>0$ is an arbitrarily small positive number. It then follows that
the iterates in the BLCD algorithm would “converge” to a contraction region
given by $A_{\gamma}$, in the sense that the iterates return to $A_{\gamma}$
infinitely often. Note that $f$ is assumed to be any nonconvex smooth
function, and there can be many contraction regions, each corresponding to a
stationary point.

  * •

When the communication error is unbiased, the gradients would diminish to $0$
and hence the BLCD algorithm would converge to a stationary point. In the case
the bias in the communication error does exist, there exists intrinsic
tradeoff between the size of the contraction region and
$\eta\lVert\mathbb{E}_{t}[\epsilon_{t}]\rVert_{2}$. When the learning rate
$\gamma$ is small, the right hand side of (15) would small, but $\eta$ can be
large, and vice verse. It makes sense to choose a fixed learning rate that
would make $\eta$ small. In this way, the gradients in the BLCD algorithm
would “concentrate” around a (small) scaled version of the bias.

  * •

Finally, the impact of gradient sparsification is captured by $\delta$. For
instance, when (randomly) uniform selection is used, $\delta=\frac{k}{d}$. We
will elaborate on this in Section VI.

Further, we have the following corollary.

######  Corollary 1.

Under Assumptions 1, 2, and 3, we have that if
$\operatorname{\mathbb{E}}_{t}[\epsilon_{t}]=0$ and
$\gamma=\frac{1}{\sqrt{T+1}}$, the BLCD algorithm converges to a stationary
point and satisfies that

$\displaystyle\frac{1}{T+1}\sum_{t=0}^{T}\lVert\nabla f(w_{t})\rVert_{2}^{2}$
$\displaystyle\leq\frac{1}{2-\frac{\rho}{\sqrt{T+1}}}\left\\{\frac{2(f(w_{0})-f^{*})}{\sqrt{T+1}}+\frac{2LG^{2}}{\sqrt{T+1}}\left(\frac{L}{\rho}\frac{2(1-\delta)}{\delta^{2}}+\frac{1}{2}\right)\right.$
$\displaystyle\left.\hskip 36.135pt+\
\frac{L}{T+1}\sum_{t=0}^{T}\underbrace{\mathbb{E}_{t}[\lVert\epsilon_{t}\rVert_{2}^{2}]}_{\text{MSE}}\right\\}$
(16)

##  V Communication Error Minimization via Joint Optimization of Power
Allocation and Learning Rates

Theorem 1 reveals that the communication error has a significant impact on the
convergence behavior of the BLCD algorithm. In this section, we turn our
attention to minimizing the communication error (in term of MSE and bias) via
joint optimization of power allocation and learning rates.

Without loss of generality, we focus on iteration $t$ (with abuse of notation,
we omit $t$ in the notation for simplicity). Recall that the coordinate
updates in the BLCD algorithm, sent by different devices over the same
subcarrier, are received by the edge server as an aggregate sum, which is used
to estimate the gradient value in that specific dimension. We denote the power
coefficients and estimators as
$\bm{b}\triangleq[b_{11},b_{12},\ldots,b_{1M},b_{21},\ldots,b_{KM}]$ and
$\bm{\alpha}\triangleq[\vec{\alpha}_{1},\ldots,\vec{\alpha}_{K}]$. In each
round, each sender device optimizes its power allocation for transmitting the
selected coordinates of their updates over the $K$ subcarriers, aiming to
achieve the best convergence rate. We assume that the perfect channel state
information is available at the corresponding transmitter, i.e.,
$\bm{h}_{m}=[h_{1m},\ldots,h_{Km}]^{\top}$ is available at the sender $m$
only.

Based on (10), the mean squared error of the communication error in iteration
$t$ is given by

$\mathbb{E}_{t}[\lVert\epsilon_{t}\rVert_{2}^{2}]=\operatorname{\mathbb{E}}\bigg{[}\left\lVert\widehat{G}_{t}(w_{t})-G_{t}(w_{t})\right\rVert^{2}\bigg{]}$
(17)

where the expectation is taken over the channel noise. For convenience, we
denote $\mathbb{E}_{t}[\lVert\epsilon_{t}\rVert_{2}^{2}]$ as
$\mbox{{MSE}}_{1}$, and after some algebra, it can be rewritten as the sum of
the variance and the square of the bias:

$\displaystyle\mbox{{MSE}}_{1}(\bm{\alpha},\bm{b})=\sum_{k=1}^{K}\bigg{[}\underbrace{\sum_{m=1}^{M}\left(\alpha_{k}b_{km}h_{km}-\frac{1}{M}\right)x_{km}}_{\mbox{bias
in $k$th
coordinate}}\bigg{]}^{2}+\underbrace{\sum_{k=1}^{K}\sigma^{2}\alpha^{2}_{k}}_{\mbox{variance\
}}$ (18)

Recall that $\\{\alpha_{k}\\}_{k=1}^{K}$ are intimately related to the
learning rates for the $K$ coordinates, making the learning rate effectively
$\\{\gamma\alpha_{k}\\}_{k=1}^{K}$.

###  V-A Centralized Solutions to Minimizing MSE (Scheme 1)

In light of the above, we can cast the MSE minimization problem as a learning-
driven joint power allocation and learning rate problem, given by

$\displaystyle\textbf{P1:\ }\min_{\bm{\alpha},\bm{b}}\quad$
$\displaystyle\mbox{{MSE}}_{1}(\bm{\alpha},\bm{b})$ (19) s.t.
$\displaystyle\sum_{k=1}^{K}\left\lvert b_{km}x_{km}\right\rvert^{2}\leq
E_{m},\quad\forall m$ (20) $\displaystyle b_{km}\geq 0,\ \alpha_{k}\geq
0\quad\forall k,m$ (21)

which minimizes the MSE for every round.

The above formulated problem is non-convex because the objective function
involves the product of variables. Nevertheless, it is biconvex, i.e., for one
of the variables being fixed, the problem is convex for the other one. In
general, we can solve the above bi-convex optimization problem in the same
spirit as in the EM algorithm, by taking the following two steps, each
optimizing over a single variable, iteratively:

$\displaystyle\textbf{P1-a:}\ \min_{\bm{\alpha}}$
$\displaystyle\mbox{{MSE}}_{1}(\bm{\alpha},\bm{b})\quad\textrm{s.t.}\quad\alpha_{k}\geq
0,\quad\forall k$ $\displaystyle\textbf{P1-b:}\min_{\bm{b}}$
$\displaystyle\mbox{{MSE}}_{11}(\bm{\alpha},\bm{b})$ s.t.
$\displaystyle\sum_{k=1}^{K}\left\lvert b_{km}x_{km}\right\rvert^{2}\leq
E_{m}\ \ \forall m,\quad b_{km}\geq 0\ \ \forall k,m.$

Since (P1-a) is unconstrained,for given $\\{b_{km}\\}$, the optimal solution
to (P1-a) is given by

$\displaystyle\alpha^{*}_{k}\\!=\\!\max\left\\{\frac{\big{(}\sum_{m=1}^{M}x_{km}\big{)}\big{(}\sum_{m=1}^{M}b_{km}h_{km}x_{km}\big{)}}{M\big{[}\sigma^{2}+\big{(}\sum_{m=1}^{M}b_{km}h_{km}x_{km}\big{)}^{2}\big{]}},0\right\\}.$
(22)

Then, we can solve (P1-b) by optimizing $\bm{b}$ only. Solving the sub-
problems (P1-a) and (P1-b) iteratively leads to a local minimum, however, not
necessarily to the global solution.

Observe that the above solution requires the global knowledge of $x_{km}$’s
and $h_{km}$’s of all devices, which is difficult to implement in practice. We
will treat it as a benchmark only. Next, we turn our attention to developing
distributed sub-optimal solutions.

###  V-B Distributed Solutions towards Zero Bias and Variance Reduction
(Scheme 2)

As noted above, the centralized solution to (P1) requires the global knowledge
of $x_{km}$’s and $h_{km}$’s and hence is not amenable to implementation.
Further, minimizing the MSE of the communication error does not necessarily
amount to minimizing the bias therein since there exists tradeoffs between
bias and variance. Thus motivated, we next focus on devising distributed sub-
optimal solutions which can drive the bias in the communication error to
(close to) zero, and then reduce the corresponding variance as much as
possible.

Specifically, observe from (18) that the minimization of MSE cost does not
necessarily ensure $\hat{G}$ to be an unbiased estimator, due to the intrinsic
tradeoff between bias and variance. To this end, we take a sub-optimal
approach where the optimization problem is decomposed into two subproblems. In
the subproblem at the transmitters, each device $m$ utilizes its available
power and local gradient/channel information to compute a power allocation
policy in terms of $\\{b_{1m},b_{2m},\dots,b_{Km}\\}$. In the subproblem at
the receiver, the receiver finds the best possible $\alpha_{k}$ for all
$k=1,\dots,K$. Another complication is that due to the power constraints at
individual devices, it is not always feasible to achieve unbiased estimators
of the gradient signal across the coordinates. Nevertheless, for given power
constraints, one can achieved unbiased estimators of a scaled down version of
the coordinates of the gradient signal. In light of this, we formulate the
optimization problem at each device (transmitter) $m$ to ensure an unbiased
estimator of a scaled version $\zeta_{m}$ of the transmitted coordinates, as
follows:

Device m: $\displaystyle\underset{\\{b_{km}\\}_{k=1:K}}{\max}~{}\zeta_{m}$
(23) $\displaystyle\text{s.t.}\hskip 7.22743pt\sum_{k=1}^{K}b_{km}^{2}$
$\displaystyle x_{km}^{2}\leq E_{m},\ \ b_{km}\geq 0,$ (24)
$\displaystyle\zeta_{m}x_{km}$ $\displaystyle-b_{km}h_{km}x_{km}=0,$
$\displaystyle\forall k=1,\dots,K,$ (25)

where maximizing $\zeta_{m}$ amounts to maximizing the corresponding SNR (and
hence improving the gradient estimation accuracy). The first constraint in the
above is the power constraint, and the second constraint is imposed to ensure
that there is no bias of the same scaled version of the transmitted signals
across the dimensions for user $m$. The power allocation solution can be found
using Karush-Kuhn-Tucker (KKT) conditions as follows:

$\displaystyle\zeta_{m}^{*}=\sqrt{\frac{E_{m}}{\sum_{k=1}^{K}\frac{x_{km}^{2}}{h_{km}^{2}}}},~{}b_{km}^{*}=\frac{\zeta^{*}_{m}}{h_{km}},~{}\forall
k.$ (26)

Observe that using the obtained power allocation policy in (66), all $K$
transmitted coordinates for device $m$ have the same scaling factor
$\zeta_{m}$. Next, we will ensure zero bias by choosing the right
$\bm{\alpha}$ for gradient estimation at the receiver, which can be obtained
by solving the following optimization problem since all transmitted gradient
signals are superimposed via the over-the-air transmission:

$\displaystyle\mbox{\bf Receiver side:}\
\underset{\\{\alpha_{k}\\},}{\min}~{}$
$\displaystyle\sum_{k=1}^{K}{\nu}_{k}^{2}(\alpha_{k},\\{b_{km}^{*}\\})$ (27)
$\displaystyle\text{s.t.}~{}e_{k}(\alpha_{k},\\{b_{km}^{*}\\})=0,$
$\displaystyle~{}~{}~{}\alpha_{k}\geq 0,\forall k=1,\dots,K,$ (28)

where $e_{k}$ and $\nu_{k}^{2}$ denote the bias and variance components, given
as follows:

$\displaystyle e_{k}(\alpha_{k},\\{b_{km}^{*}\\})$
$\displaystyle=\alpha_{k}\left(\sum_{m=1}^{M}\zeta_{m}^{*}x_{km}\right)-\frac{1}{M}\sum_{m=1}^{M}x_{km},$
$\displaystyle\nu_{k}^{2}(\alpha_{k},\\{b_{km}^{*}\\})$
$\displaystyle=\alpha_{k}^{2}\sigma^{2},$ (29)

for all $k=1,\dots,K$. For given $\\{\zeta_{m}^{*}\\}$, it is easy to see that

$\displaystyle\alpha_{k}^{*}=\frac{\frac{1}{M}\sum_{m=1}^{M}x_{km}}{\sum_{m=1}^{M}\zeta_{m}^{*}x_{km}}\simeq\frac{1}{\sum_{m=1}^{M}\zeta_{m}^{*}},\
\ \forall k.$ (30)

We note that in the above, from an implementation point of view, since
$\\{x_{km}\\}$ is not available at the receiver, it is sensible to set
$\alpha_{k}^{\dagger}\simeq\frac{1}{\sum_{m=1}^{M}\zeta_{m}^{*}}$. Further,
$\\{\zeta_{m}^{*}\\}$ is not readily available at the receiver either.
Nevertheless, since there is only one parameter $\zeta_{m}^{*}$ from each
sender $m$, the sum $\sum_{m=1}^{M}\zeta_{m}^{*}$ can be sent over a control
channel to the receiver to compute $\alpha_{k}^{\dagger}$. It is worth noting
that in general the bias exists even if $E_{m}$ is the same for all senders.

Next, we take a closer look at the case when the number of subchannels $K$ is
large (which is often the case in practice). Suppose that $\\{x_{km}\\}$ are
i.i.d. across subchannels and users, and so are $\\{h_{km}\\}$. We can then
simplify $\zeta_{m}^{*}$ further. For ease of exposition, we denote
$\mathbb{E}[x_{km}^{2}]=\varphi^{2}+\bar{x}^{2}$ and
$\mathbb{E}\left[\frac{1}{h_{km}^{2}}\right]=\varpi^{2}$. When $K$ is large,
for every user $m$ we have that:

$\displaystyle\zeta_{m}^{*}=\frac{\sqrt{E_{m}}}{\sqrt{\sum_{k=1}^{K}\frac{x_{km}^{2}}{h_{km}^{2}}}}\underset{\begin{subarray}{c}\text{when
$K$}\\\ \text{is
large}\end{subarray}}{\Longrightarrow}\zeta_{m}^{*}\approx\frac{\sqrt{E_{m}}}{\sqrt{K(\varphi^{2}+\bar{x}^{2})\varpi^{2}}}$
(31)

As a result, the bias and variance for each dimension $k$ could be written as,

$\displaystyle e_{k}(\alpha_{k}^{*},\\{b_{km}^{*}\\})$
$\displaystyle=\sum_{m=1}^{M}\left[\frac{\sqrt{E_{m}}}{\sum_{m=1}^{M}\sqrt{E_{m}}}-\frac{1}{M}\right]x_{km},\forall
k.$ (32) $\displaystyle{\nu}_{k}^{2}$
$\displaystyle=\frac{K\varpi^{2}(\varphi^{2}+\bar{x}^{2})}{\left(\sum_{m=1}^{M}\sqrt{E_{m}}\right)^{2}}\sigma^{2},\forall
k.$ (33)

Observe that when $E_{m}$ is the same across the senders, the bias term
$\mathbb{E}_{t}[\epsilon_{t}]=\mathbf{0}$ in the above setting according to
(32).

###  V-C A User-centric Approach Using Single-User Solution (Scheme 3)

In this section, we consider a suboptimal user-centric approach, which
provides insight on the power allocation across the subcarriers from a single
device perspective. We formulate the single device (say user $m$) problem as

P2:
$\displaystyle\min_{\\{b_{km}\\},\\{\alpha_{k}\\}}~{}\sum_{k=1}^{K}\bigg{[}\big{(}\alpha_{k}b_{km}h_{km}-1\big{)}x_{km}\bigg{]}^{2}+\sigma^{2}\sum_{k=1}^{K}\alpha^{2}_{k}$
$\displaystyle\textrm{s.t.}~{}\sum_{k=1}^{K}\left\lvert
b_{km}x_{km}\right\rvert^{2}\leq E_{m};\ b_{k}\geq 0,\ \alpha_{k}\geq
0,\forall k.$

######  Theorem 2.

The optimal solution $\\{b_{km}^{*},\alpha_{k}^{*}\\}$ to (P2) is given by

$(b^{*}_{km})^{2}=\bigg{[}\sqrt{\frac{\sigma^{2}}{\lambda
x_{km}^{2}h_{km}^{2}}}-\frac{\sigma^{2}}{h_{km}^{2}x_{km}^{2}}\bigg{]}^{+},\forall
k,$ (34)
$\displaystyle\alpha^{*}_{k}=\frac{b^{*}_{km}h_{km}x^{2}_{k}}{\sigma^{2}+(b^{*}_{km})^{2}h^{2}_{km}x^{2}_{km}},~{}~{}~{}~{}\forall
k,$ (35)

where $\lambda_{m}$ is a key parameter determining the waterfilling level:

$\sum_{k=1}^{K}\bigg{[}\sqrt{\frac{1}{\lambda_{m}}}\sqrt{\frac{x_{km}^{2}\sigma^{2}}{h_{km}^{2}}}-\frac{\sigma^{2}}{h_{km}^{2}}\bigg{]}^{+}=E_{m}.$
(36)

The proof of Theorem 2 is omitted due to space limitation. Observe that
Theorem 2 reveals that the larger the gradient value (and the smaller channel
gain) in one subcarrier, the higher power the it should be allocated to in
general, and that $\\{x_{km}/h_{km)}\\}$ can be used to compute the water
level for applying the water filling policy.

Based on the above result, in the multi-user setting, each device can adopt
the above single-user power allocation solution as given in Theorem 2. This
solution can be applied individually without requiring any coordination
between devices.

Next, we take a closer look at the case when the number of subchannels $K$ is
large. Let $\bar{E}_{m}$ denote the average power constraint per subcarrier.
When $K$ is large, after some algebra, the optimization problem P2 can be
further approximated as follows:

$\displaystyle\textbf{P3: }\underset{{b}_{km}}{\min}$
$\displaystyle~{}\mathbb{E}\left[\frac{{{x}^{2}_{km}}\sigma^{2}}{{{b}^{2}_{km}}{{h}^{2}_{km}}{{x}^{2}_{km}}+\sigma^{2}}\right]$
s.t.
$\displaystyle~{}\mathbb{E}\left[{b}^{2}_{km}{{x}^{2}_{km}}\right]\leq\bar{E}_{m},~{}{b}_{km}\geq
0,$ (37)

where the expectation is taken with respect to $\\{{h}_{km}\\}$ and
$\\{{x}_{km}\\}$.

The solution for $k=1,\dots,K$ is obtained as follows:

$\displaystyle b_{km}^{*}=\sqrt{\left[\frac{\sigma\lvert
x_{km}\rvert^{-1}}{h_{km}\sqrt{\lambda_{m}}}-\frac{\sigma^{2}}{x_{km}^{2}h_{km}^{2}}\right]^{+}}$
(38)
$\displaystyle\lambda_{m}<\frac{h_{km}^{2}x_{km}^{2}}{\sigma_{k}^{2}}\Rightarrow
b^{*}_{km}>0$ (39)

We can compute the bias and the variance accordingly.

##  VI Coordinate Selection for Bandlimited Coordinate Descent Algorithms

The selection of which coordinates to operate on is crucial to the performance
of sparsified SGD algorithms. It is not hard to see that selecting the top-$k$
(in absolute value) coordinates of the sum of the gradients provides the best
performance. However, in practice it may not always be feasible to obtain
top-$k$ of the sum of the gradients, and in fact there are different solutions
for selecting $k$ dimensions with large absolute values; see e.g., [22, 27].
Note that each device individually transmitting top-$k$ coordinates of their
local gradients is not applicable to the scenario of over-the-air
communications considered here. Sequential device-to-device transmissions
provides an alternative approach [28], but these techniques are likely to
require more bandwidth with wireless connection.

Another approach that is considered is the use of compression and/or sketching
for the gradients to be transmitted. For instance, in [22], a system that
updates SGD via decompressing the compressed gradients transmitted through
over-the-air communication is examined. To the best of our knowledge, such
techniques do not come with rigorous convergence guarantees. A similar
approach is taken in [27], where the sketched gradients are transmitted
through an error-free medium and these are then used to obtain top-$k$
coordinates; the devices next simply transmit the selected coordinates.
Although such an approach can be taken with over-the-air computing since only
the summation of the sketched gradients is necessary; this requires the
transmission of $\mathcal{O}(k\log d)$ dimensions. To provide guarantees with
such an approach $\mathcal{O}(k\log d+k)$ up-link transmissions are needed.
Alternatively, uniformly selected $\mathcal{O}(k\log d+k)$ coordinates can be
transmitted with similar bandwidth and energy requirements. For the practical
learning models with non-sparse updates, uniform coordinate selection tend to
perform better. Moreover, the common $K$ dimensions can be selected uniformly
via synchronized pseudo-random number generators without any information
transfer. To summarize, uniform selection of the coordinates is more
attractive based on the energy, bandwidth and implementation considerations
compared to the methods aiming to recover top-$k$ coordinates; indeed, this is
the approach we adopt.

##  VII Experimental Results

In this section, we evaluate the accuracy and convergence performance of the
BLCD algorithm, when using one of the following three schemes for power
allocation and learning rate selection (aiming to minimize the impact of
communication error): 1) the bi-convex program based solution (Scheme 1), 2)
the distributed solution towards zero bias in Section V. (Scheme 2); 3) the
single-user solution (Scheme 3). We use the communication error free scheme as
the baseline to evaluate the performance degradation. We also consider the
naive scheme (Scheme 4) using equal power allocation for all dimensions, i.e.,
$b_{km}=\sqrt{E/{\sum_{k=1}^{K}x^{2}_{km}}}$.

In our first experiment, we consider a simple single layer neural network
trained on the MNIST dataset. The network consists of two 2-D convolutional
layers with filter size $5\times 5$ followed by a single fully connected layer
and it has 7840 parameters. $K=64$ dimensions are uniformly selected as the
support of the sparse gradient transmissions. For convenience, we define
$E_{avg}$ as the average sum of the energy (of all devices) per dimension
normalized by the channel noise variance, i.e.,
$E_{avg}=EM\operatorname{\mathbb{E}}[h_{km}^{2}]/K\sigma^{2}.$ Without loss of
generality, we take the variance of the channel noise as $\sigma^{2}=1$ and
$\\{h_{km}\\}$ are independent and identically distributed Rayleigh random
variables with mean $1$. The changes on $E_{avg}$ simply amount to different
SNR values. In Fig. 4, we take $K=64$, $M=8$, batch size $4$ to calculate each
gradient, and the learning rate $\gamma=0.01$. In the second experiment, we
expand the model to a more sophisticated $5$-layer neural network and an
$18$-layer ResNet [29] with $61706$ and $11175370$ parameters, respectively.
The $5$-layer network consists of two 2-D convolutional layers with filter
size $5\times 5$ followed by three fully connected layers. In all experiments,
we have used a learning rate of $0.01$. the local dataset of each worker is
randomly selected from the entire MNIST dataset. We use $10$ workers with
varying batch sizes and we utilize $K=1024$ sub-channels for sparse gradient
signal transmission.

It can be seen from Fig. 4 that in the presence of the communication error,
the centralized solution (Scheme 1) based on bi-convex programming converges
quickly and performs the best, and it can achieve accuracy close to the ideal
error-free scenario. Further, the distributed solution (Scheme 2) can
eventually approach the performance of Scheme 1, but the single-user solution
(Scheme 3) performs poorly, so does the naive scheme using equal power
allocation (Scheme 4). Clearly, there exists significant gap between its
resulting accuracy and that in the error-free case, and this is because the
bias in Scheme 3 is more significantly.

Figure 3: Testing accuracy over training iterations for $\alpha_{k}=1/8$,
$E_{avg}=0.1$ and a batch size of $4$. Training model consists of a single
layer neural network with 7840 differentiable parameters.

Figure 4: Testing accuracy over training iterations for $10$ workers and a
batch size of $256$. Training model consists of a $5$-layer deep neural
network with 61706 differentiable parameters.

Next, Figures 4, 6 and 6 depict the results in the second experiment using
much larger-scale deep neural networks. It can be observed from Figs. 4, 6 and
6 that the SNR can have significant impact on the final accuracy. As expected,
the convergence on the ResNet network is slower in comparison to other DNNs
due to the huge number of parameters and small batch size. Nevertheless, it is
clear that the learning accuracy improves significantly at high SNR. (The
solution of the distributed algorithm for $E_{avg}=10$ is omitted in Fig. 6,
since it is indistinguishably close to error-free solution.) It is interesting
to observe that when the SNR increases, the distributed solution (Scheme 2)
can achieve accuracy close to the ideal error-free case, but the single-user
solution (Scheme 3) would not. It is worth noting that due to the
computational complexity of bi-convex programming in this large-scale case,
Scheme 4 could be solved effectively (we did not present it here). Further,
the batch size at each worker can impact the convergence rate, but does not
impact the final accuracy.

Figure 5: Testing accuracy over training iterations for $10$ workers and a
batch size of $4$. Training model consists of a $5$-layer deep neural network
with 61706 differentiable parameters.

Figure 6: Testing accuracy over training iterations for $10$ devices and a
batch size of $4$. Training model consists of an $18$-layer ResNet network
with more than 11 million differentiable parameters.

##  VIII Conclusions

In this paper, we consider a many-to-one wireless architecture for distributed
learning at the network edge, where multiple edge devices collaboratively
train a machine learning model, using local data, through a wireless channel.
Observing the unreliable nature of wireless connectivity, we design an
integrated communication and learning scheme, where the local updates at edge
devices are carefully crafted and compressed to match the wireless
communication resources available. Specifically, we propose SGD-based
bandlimited coordinate descent algorithms employing over-the-air computing, in
which a subset of k-coordinates of the gradient updates across edge devices
are selected by the receiver in each iteration and then transmitted
simultaneously over k sub-carriers. We analyze the convergence of the
algorithms proposed, and characterize the effect of the communication error.
Further, we study joint optimization of power allocation and learning rates
therein to maximize the convergence rate. Our findings reveal that optimal
power allocation across different sub-carriers should take into account both
the gradient values and channel conditions. We then develop sub-optimal
solutions amenable to implementation and verify our findings through numerical
experiments.

## Acknowledgements

The authors thank Gautam Dasarathy for stimulating discussion in the early
stage of this work. This work is supported in part by NSF Grants CNS-2003081,
CNS-CNS-2003111, CPS-1739344 and ONR YIP N00014-19-1-2217.

## References

  * [1] G. Zhu, D. Liu, Y. Du, C. You, J. Zhang, and K. Huang, “Towards an intelligent edge: Wireless communication meets machine learning,” _arXiv preprint arXiv:1809.00343_ , 2018.
  * [2] M. Goldenbaum and S. Stanczak, “Robust analog function computation via wireless multiple-access channels,” _IEEE Transactions on Communications_ , vol. 61, no. 9, pp. 3863–3877, 2013.
  * [3] O. Abari, H. Rahul, and D. Katabi, “Over-the-air function computation in sensor networks,” _arXiv preprint arXiv:1612.02307_ , 2016.
  * [4] D. Alistarh, D. Grubic, J. Li, R. Tomioka, and M. Vojnovic, “Qsgd: Communication-efficient sgd via gradient quantization and encoding,” in _Advances in Neural Information Processing Systems_ , 2017, pp. 1709–1720.
  * [5] W. Wen, C. Xu, F. Yan, C. Wu, Y. Wang, Y. Chen, and H. Li, “Terngrad: Ternary gradients to reduce communication in distributed deep learning,” in _Advances in Neural Information Processing Systems_ , 2017, pp. 1509–1519.
  * [6] J. Bernstein, Y.-X. Wang, K. Azizzadenesheli, and A. Anandkumar, “Signsgd: Compressed optimisation for non-convex problems,” in _International Conference on Machine Learning_ , 2018, pp. 559–568.
  * [7] J. Wu, W. Huang, J. Huang, and T. Zhang, “Error compensated quantized sgd and its applications to large-scale distributed optimization,” in _International Conference on Machine Learning_ , 2018, pp. 5321–5329.
  * [8] A. F. Aji and K. Heafield, “Sparse communication for distributed gradient descent,” in _Proceedings of the 2017 Conference on Empirical Methods in Natural Language Processing_ , 2017, pp. 440–445.
  * [9] S. U. Stich, J.-B. Cordonnier, and M. Jaggi, “Sparsified sgd with memory,” in _Advances in Neural Information Processing Systems_ , 2018, pp. 4447–4458.
  * [10] D. Alistarh, T. Hoefler, M. Johansson, N. Konstantinov, S. Khirirat, and C. Renggli, “The convergence of sparsified gradient methods,” in _Advances in Neural Information Processing Systems_ , 2018, pp. 5973–5983.
  * [11] 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,” _arXiv preprint arXiv:1610.05492_ , 2016.
  * [12] S. U. Stich, “Local sgd converges fast and communicates little,” in _ICLR 2019 International Conference on Learning Representations_ , 2019.
  * [13] J. Dong, Y. Shi, and Z. Ding, “Blind over-the-air computation and data fusion via provable wirtinger flow,” _arXiv preprint arXiv:1811.04644_ , 2018.
  * [14] W. Liu and X. Zang, “Over-the-air computation systems: Optimization, analysis and scaling laws,” _arXiv preprint arXiv:1909.00329_ , 2019.
  * [15] D. Wen, G. Zhu, and K. Huang, “Reduced-dimension design of MIMO over-the-air computing for data aggregation in clustered iot networks,” _IEEE Transactions on Wireless Communications_ , vol. 18, no. 11, pp. 5255–5268, 2019\.
  * [16] G. Zhu and K. Huang, “MIMO over-the-air computation for high-mobility multi-modal sensing,” _IEEE Internet of Things Journal_ , 2018.
  * [17] X. Cao, G. Zhu, J. Xu, and K. Huang, “Optimal power control for over-the-air computation in fading channels,” _arXiv preprint arXiv:1906.06858_ , 2019\.
  * [18] G. Zhu, Y. Wang, and K. Huang, “Broadband analog aggregation for low-latency federated edge learning,” _IEEE Transactions on Wireless Communications_ , 2019.
  * [19] K. Yang, T. Jiang, Y. Shi, and Z. Ding, “Federated learning via over-the-air computation,” _arXiv preprint arXiv:1812.11750_ , 2018.
  * [20] Q. Zeng, Y. Du, K. K. Leung, and K. Huang, “Energy-efficient radio resource allocation for federated edge learning,” _arXiv preprint arXiv:1907.06040_ , 2019.
  * [21] M. M. Amiri and D. Gündüz, “Machine learning at the wireless edge: Distributed stochastic gradient descent over-the-air,” _arXiv preprint arXiv:1901.00844_ , 2019.
  * [22] ——, “Federated learning over wireless fading channels,” _arXiv preprint arXiv:1907.09769_ , 2019.
  * [23] ——, “Over-the-air machine learning at the wireless edge,” in _2019 IEEE 20th International Workshop on Signal Processing Advances in Wireless Communications (SPAWC)_. IEEE, 2019, pp. 1–5.
  * [24] J.-H. Ahn, O. Simeone, and J. Kang, “Wireless federated distillation for distributed edge learning with heterogeneous data,” in _2019 IEEE 30th Annual International Symposium on Personal, Indoor and Mobile Radio Communications (PIMRC)_. IEEE, 2019, pp. 1–6.
  * [25] T. Sery and K. Cohen, “On analog gradient descent learning over multiple access fading channels,” _arXiv preprint arXiv:1908.07463_ , 2019.
  * [26] S. P. Karimireddy, Q. Rebjock, S. U. Stich, and M. Jaggi, “Error feedback fixes signsgd and other gradient compression schemes,” _arXiv preprint arXiv:1901.09847_ , 2019.
  * [27] N. Ivkin, D. Rothchild, E. Ullah, V. Braverman, I. Stoica, and R. Arora, “Communication-efficient distributed sgd with sketching,” _arXiv preprint arXiv:1903.04488_ , 2019.
  * [28] S. Shi, Q. Wang, K. Zhao, Z. Tang, Y. Wang, X. Huang, and X. Chu, “A distributed synchronous sgd algorithm with global top-k sparsification for low bandwidth networks,” in _2019 IEEE 39th International Conference on Distributed Computing Systems (ICDCS)_. IEEE, 2019, pp. 2238–2247.
  * [29] K. He, X. Zhang, S. Ren, and J. Sun, “Deep residual learning for image recognition,” _arXiv preprint arXiv:1512.03385_ , 2015.

##  Appendix A Proofs

###  A-A Proof of Theorem 1

We here restate equations (6) and (7) as follows:

$\displaystyle{w}_{t+1}=$ $\displaystyle w_{t}-[C_{t}(\gamma
g_{t}(w_{t})+r_{t})+\epsilon_{t}]$ (40) $\displaystyle\hat{w}_{t+1}=$
$\displaystyle w_{t}-C_{t}(\gamma g_{t}(w_{t})+r_{t})$ (41)

It is clear that ${w}_{t+1}=\hat{w}_{t+1}-\epsilon_{t}$. For convenience, we
define $\tilde{w}_{t}={w}_{t}-r_{t}=\hat{w}_{t}-r_{t}-\epsilon_{t-1}$. It can
be shown that $\tilde{w}_{t+1}=\tilde{w}_{t}-\gamma
g_{t}(w_{t})-\epsilon_{t}$.

$\displaystyle\mathbb{E}_{t}[f(\tilde{w}_{t+1})]\leq$ $\displaystyle
f(\tilde{w}_{t})+<\nabla
f(\tilde{w}_{t}),\mathbb{E}_{t}[\tilde{w}_{t+1}-\tilde{w}_{t}]>+\frac{L}{2}\mathbb{E}_{t}[||\tilde{w}_{t+1}-\tilde{w}_{t}||^{2}]$
(42) $\displaystyle=$ $\displaystyle f(\tilde{w}_{t})-<\nabla
f(\tilde{w}_{t}),\gamma\mathbb{E}_{t}[g_{t}(w_{t})]+\mathbb{E}_{t}[\epsilon_{t}]>+\frac{L}{2}\mathbb{E}_{t}[||\gamma
g_{t}(w_{t})||^{2}]$ $\displaystyle\hskip
36.135pt+\frac{L}{2}\mathbb{E}_{t}[||\epsilon_{t}||^{2}]+L\mathbb{E}_{t}[<\gamma
g_{t}(w_{t}),\epsilon_{t}>]$ (43) $\displaystyle=$ $\displaystyle
f(\tilde{w}_{t})-<\nabla
f({w}_{t}),\gamma\mathbb{E}_{t}[g_{t}(w_{t})]+\mathbb{E}_{t}[\epsilon_{t}]>-<\nabla
f(\tilde{w}_{t})-\nabla
f({w}_{t}),\gamma\mathbb{E}_{t}[g_{t}(w_{t})]+\mathbb{E}_{t}[\epsilon_{t}]>$
$\displaystyle\hskip
36.135pt+\frac{L}{2}\mathbb{E}_{t}[||\epsilon_{t}||_{2}^{2}]+L\mathbb{E}_{t}[<\gamma
g_{t}(w_{t}),\epsilon_{t}>]+\frac{L}{2}\mathbb{E}_{t}[||\gamma
g_{t}(w_{t})||^{2}$ (44) $\displaystyle\leq$ $\displaystyle
f(\tilde{w}_{t})-\gamma\lVert\nabla f(w_{t})\rVert_{2}^{2}-\langle\nabla
f(w_{t}),\mathbb{E}_{t}[\epsilon_{t}]\rangle+\frac{\rho}{2}\lVert\gamma\nabla
f(w_{t})+\mathbb{E}_{t}[\epsilon_{t}]\rVert_{2}^{2}+\frac{L^{2}}{2\rho}\mathbb{E}_{t}[\lVert
r_{t}\rVert_{2}^{2}]$ $\displaystyle\hskip
36.135pt+\frac{L}{2}\mathbb{E}_{t}[\lVert\epsilon_{t}\rVert_{2}^{2}]+L\langle\nabla
f(w_{t}),\mathbb{E}_{t}[\epsilon_{t}]\rangle+\frac{L\gamma^{2}}{2}\mathbb{E}_{t}\lVert
g_{t}(w_{t})\rVert_{2}^{2}$ (45) $\displaystyle\leq$ $\displaystyle
f(\tilde{w}_{t})-\gamma\lVert\nabla f(w_{t})\rVert_{2}^{2}+(L-1)\lVert\nabla
f(w_{t})\rVert\lVert\mathbb{E}_{t}[\epsilon_{t}]\rVert+\frac{\rho}{2}\left(\gamma^{2}\lVert\nabla
f(w_{t})\rVert_{2}^{2}+\lVert\mathbb{E}_{t}[\epsilon_{t}]\rVert_{2}^{2}+2\gamma\langle\nabla
f(w_{t}),\mathbb{E}_{t}[\epsilon_{t}]\rangle\right)$ $\displaystyle\hskip
36.135pt+\frac{L^{2}}{2\rho}\mathbb{E}_{t}[\lVert
r_{t}\rVert_{2}^{2}]+\frac{L}{2}\mathbb{E}_{t}[\lVert\epsilon_{t}\rVert_{2}^{2}]+\frac{L\gamma^{2}}{2}G^{2}$
(46) $\displaystyle\leq$ $\displaystyle f(\tilde{w}_{t})-\gamma\lVert\nabla
f(w_{t})\rVert_{2}^{2}+(L-1+2\gamma)\lVert\nabla
f(w_{t})\rVert\lVert\mathbb{E}_{t}[\epsilon_{t}]\rVert+\frac{\gamma^{2}\rho}{2}\lVert\nabla
f(w_{t})\rVert_{2}^{2}$ $\displaystyle\hskip
36.135pt+\frac{L^{2}}{2\rho}\mathbb{E}_{t}[\lVert
r_{t}\rVert_{2}^{2}]+\lVert\mathbb{E}_{t}[\epsilon_{t}]\rVert_{2}^{2}+\frac{L}{2}\mathbb{E}_{t}[\lVert\epsilon_{t}\rVert_{2}^{2}]+\frac{L\gamma^{2}}{2}G^{2}$
(47) $\displaystyle=$ $\displaystyle
f(\tilde{w}_{t})-\gamma\left[1-\frac{\rho}{2}\gamma\right]\lVert\nabla
f(w_{t})\rVert_{2}^{2}+(L-1+2\gamma)\lVert\nabla
f(w_{t})\rVert\lVert\mathbb{E}_{t}[\epsilon_{t}]\rVert$ $\displaystyle\hskip
36.135pt+\frac{L^{2}}{2\rho}\mathbb{E}_{t}[\lVert
r_{t}\rVert_{2}^{2}]+\lVert\mathbb{E}_{t}[\epsilon_{t}]\rVert_{2}^{2}+\frac{L}{2}\mathbb{E}_{t}[\lVert\epsilon_{t}\rVert_{2}^{2}]+\frac{L\gamma^{2}}{2}G^{2}$
(48)

Based on [26], we have that

$\displaystyle\mathbb{E}_{t}[\lVert
r_{t}\rVert_{2}^{2}]\leq\frac{4(1-\delta)}{\delta^{2}}\gamma^{2}G^{2}.$ (49)

It follows that

$\displaystyle\frac{1}{T+1}\sum_{t=0}^{T}\left\\{\gamma(1-\frac{\rho}{2}\gamma)\lVert\nabla
f(w_{t})\rVert^{2}_{2}-(L-1+2\gamma)\lVert\mathbb{E}_{t}[\epsilon_{t}]\rVert\lVert\nabla
f(w_{t})\rVert\right\\}$ $\displaystyle\hskip
36.135pt\leq\frac{1}{T+1}[f(w_{0})-f^{*}]+\frac{L^{2}}{\rho}\frac{2(1-\delta)}{\delta^{2}}\gamma^{2}G^{2}+\frac{L\gamma^{2}}{2}G^{2}+\frac{1}{T+1}\sum_{t=0}^{T}\left[\lVert\mathbb{E}_{t}[\epsilon_{t}]\rVert_{2}^{2}+\frac{L}{2}\mathbb{E}_{t}[\lVert\epsilon_{t}\rVert_{2}^{2}]\right]$
(50)

Through some further algebraic manipulation, we have that

$\displaystyle\frac{1}{T+1}\sum_{t=0}^{T}\left(\lVert\nabla
f(w_{t})\rVert_{2}-\frac{L-1+2\gamma}{\gamma(2-\rho\gamma)}\lVert\mathbb{E}_{t}[\epsilon_{t}]\rVert_{2}\right)^{2}\leq\frac{2}{T+1}\frac{f(w_{0})-f^{*}}{\gamma(2-\rho\gamma)}+\left(\frac{L}{\rho}\frac{2(1-\delta)}{\delta^{2}}+\frac{1}{2}\right)\frac{2L\gamma
G^{2}}{2-\rho\gamma}$ $\displaystyle\hskip
36.135pt+\frac{1}{T+1}\sum_{t=0}^{T}\left[\frac{L}{\gamma(2-\rho\gamma)}\mathbb{E}_{t}[\lVert\epsilon_{t}\rVert_{2}^{2}]+\left(1+\frac{(L-1+2\gamma)^{2}}{\gamma^{2}(2-\rho\gamma)^{2}}\right)\lVert\mathbb{E}_{t}[\epsilon_{t}]\rVert_{2}^{2}\right]$
(51)

For convenience, let $\eta=\frac{L-1+2\gamma}{\gamma(2-\rho\gamma)}$, and
define a contraction region as follows:

$C_{\beta}=\left\\{\lVert\nabla
f(w_{t})\rVert_{2}\geq(\eta+\Delta)\lVert\mathbb{E}_{t}[\epsilon_{t}]\rVert_{2}\right\\}.$

It follows from (51) that iterates in the BLCD algorithm returns to the
contraction region infinitely often with probability one. Further, when
setting $\gamma=\frac{1}{\sqrt{T+1}}$, we have that

$\displaystyle\frac{1}{T+1}\sum_{t=0}^{T}\left(\lVert\nabla
f(w_{t})\rVert_{2}-\frac{L-1+2\gamma}{\gamma(2-\rho\gamma)}\lVert\mathbb{E}_{t}[\epsilon_{t}]\rVert_{2}\right)^{2}\leq\frac{f(w_{0})-f^{*}}{(1-\frac{1}{2}\rho\gamma)}+\frac{1}{\sqrt{T+1}}\left(\frac{L}{\rho}\frac{2(1-\delta)}{\delta^{2}}+\frac{1}{2}\right)\frac{2LG^{2}}{2-\rho\gamma}$
$\displaystyle\hskip
36.135pt+\frac{L}{(2-\rho\gamma)}\frac{1}{T+1}\sum_{t=0}^{T}\mathbb{E}_{t}[\lVert\epsilon_{t}\rVert_{2}^{2}]+\left(1+\frac{(L-1+2\gamma)^{2}}{\gamma^{2}(2-\rho\gamma)^{2}}\right)\frac{1}{T+1}\sum_{t=0}^{T}\lVert\mathbb{E}_{t}[\epsilon_{t}]\rVert_{2}^{2}$
(52)

###  A-B Solution of Problem (P1-a)

Since the problem (P1-a) is convex, the Lagrangian function is given as:

$L_{1a}(\bm{\alpha},\bm{\lambda})=\text{MSE}_{1}(\bm{b},\bm{\alpha})+\sum_{k=1}^{K}\lambda_{k}\alpha_{k}.$
(53)

Then, the Karush-Kuhn-Tucker (KKT) conditions are given as follows:

$\begin{split}\frac{\partial
L_{1a}(\bm{\alpha},\bm{\lambda})}{\partial\alpha_{k}}=&2\bigg{(}\sum_{m=1}^{M}b_{km}h_{km}x_{km}\bigg{)}\bigg{(}\sum_{m=1}^{M}\big{(}\alpha_{k}b_{km}h_{km}-\frac{1}{M}\big{)}x_{km}\bigg{)}+2\sigma^{2}\alpha_{k}+\lambda_{k}=0,\end{split}$
(54) $\lambda_{k}\geq
0,\quad\lambda_{k}\alpha^{*}_{k}=0,\quad\alpha^{*}_{k}\geq 0.$ (55)

It follows that

$\begin{split}\alpha^{*}_{k}&=\max\left\\{\frac{\big{(}\sum_{m=1}^{M}x_{km}\big{)}\big{(}\sum_{m=1}^{M}b_{km}h_{km}x_{km}\big{)}}{M\big{[}\sigma^{2}+\big{(}\sum_{m=1}^{M}b_{km}h_{km}x_{km}\big{)}^{2}\big{]}},0\right\\}\triangleq\max\left\\{\frac{\bar{x}_{k}\beta_{k}}{\sigma^{2}+\beta_{k}^{2}}\right\\},\end{split}$
(56)

where the auxiliary variables are defined as
$\beta_{k}=\sum_{m=1}^{M}b_{km}h_{km}x_{km}$ and
$\bar{x}_{k}=\sum_{m=1}^{M}x_{km}/M$.

###  A-C Proof of Theorem 2

###### Proof.

Observing that the problem (P2) is defined only in terms of $b_{k}^{2}$, we
define the auxiliary variables $\tilde{b}_{k}=b_{k}^{2}$,
$\tilde{h}_{k}=h_{k}^{2}/\sigma^{2}$ and $\tilde{x}_{k}=1/x_{k}^{2}$, and re-
formulate (P2) as:

$\displaystyle\textbf{P2-1:\ }\min_{\tilde{\bm{b}}}$
$\displaystyle\sum_{k=1}^{K}(\tilde{b}_{k}\tilde{h_{k}}+\tilde{x}_{k})^{-1}$
s.t. $\displaystyle\sum_{k=1}^{K}\frac{\tilde{b}_{k}}{\tilde{x}_{k}}\leq E$
$\displaystyle\tilde{b}_{k}\geq 0,\quad\forall k,$

which is convex and can be solved in closed form. Then, we have the Lagrangian
as

$L_{22}(\tilde{\bm{b}},\bm{\lambda},\bm{\mu})=\sum_{k=1}^{K}(\tilde{b}_{k}\tilde{h_{k}}+\tilde{x}_{k})^{-1}+\lambda(\frac{\tilde{b}_{k}}{\tilde{x}_{k}}-E)-\sum_{k=1}^{K}\mu_{k}\tilde{b}_{k},$
(57)

which leads to the following KKT conditions:

$\frac{\partial
L_{12}(\tilde{\bm{b}},\bm{\lambda},\bm{\mu})}{\partial\tilde{b}_{k}}=-\tilde{h}_{k}(\tilde{b}_{k}\tilde{h}_{k}+\tilde{x}_{k})^{-2}+\frac{\lambda}{\tilde{x}_{k}}-\mu_{k}=0$
$\sum_{k=1}^{K}\frac{\tilde{b}^{*}_{k}}{\tilde{x}_{k}}\leq
E,\quad\tilde{b}^{*}_{k}\geq 0$ $\lambda\geq 0,\quad\mu_{k}\geq 0$
$\lambda\bigg{(}E-\sum_{k=1}^{K}\frac{\tilde{b}^{*}_{k}}{\tilde{x}_{k}}\bigg{)}=0,\quad\mu_{k}\tilde{b}^{*}_{k}=0.$

For $\mu_{k}$=0, and $\lambda>0$, we have that

$\tilde{b}_{k}^{*}=\max\bigg{\\{}\frac{\sqrt{\frac{\tilde{h}_{k}\tilde{x}_{k}}{\lambda}}-\tilde{x}_{k}}{\tilde{h}_{k}},0\bigg{\\}}=\bigg{[}\frac{\sqrt{\frac{\tilde{h}_{k}\tilde{x}_{k}}{\lambda}}-\tilde{x}_{k}}{\tilde{h}_{k}}\bigg{]}^{+}$
(58)

with $E=\sum_{k=1}^{K}\frac{\tilde{b}^{*}_{k}}{\tilde{x}_{k}}$. By combining
(58) and $E=\sum_{k=1}^{K}\frac{\tilde{b}^{*}_{k}}{\tilde{x}_{k}}$, we obtain
the following result:

$E=\sum_{k=1}^{K}\bigg{[}\frac{\sqrt{\tilde{h}_{k}\tilde{x}_{k}}\lambda^{\prime}-\tilde{x}_{k}}{\tilde{h}_{k}\tilde{x}_{k}}\bigg{]}^{+}$
(59)

for $\lambda^{\prime}=\sqrt{1/\lambda}$. (59) can be solved by using the
water-filling algorithm, where the solution can be found by increasing
$\lambda^{\prime}$ until the equality is satisfied. The optimal
$\lambda^{\prime}$ can be plugged int ${\tilde{b}^{*}_{k}}$ to yield
$b^{*}_{k}=\sqrt{\tilde{b}^{*}_{k}}$ as a solution to (P2). ∎

###  A-D Distributed Solutions towards Zero Bias and Variance Reduction
(Scheme 2)

The over-the-air gradient estimation requires a more comprehensive estimator
design. To this end, a generalized optimization problem is defined for
computing the optimal estimator. We define the MSE cost for the communication
error, $\epsilon_{t}$, in terms of the received signal
$\mathbf{y}=[y_{1},y_{2},\dots,y_{K}]$ as

$\displaystyle\mathbb{E}[(\hat{{G}}-{G})^{2}]$
$\displaystyle=\mathbb{E}[(\bm{\alpha}\odot\mathbf{y}-\mathbf{G})^{2}]=\frac{1}{K}\sum_{k=1}^{K}\left[\alpha_{k}\left(\sum_{m=1}^{M}b_{km}h_{km}x_{km}+n_{k}\right)-G_{k}\right]^{2}$
(60)
$\displaystyle=\frac{1}{K}\sum_{k=1}^{K}{\underbrace{\left(\mathbb{E}\left[\alpha_{k}\left(\sum_{m=1}^{M}b_{km}h_{km}x_{km}+n_{k}\right)\right]-\mathbb{E}\left[\frac{1}{M}\sum_{m=1}^{M}x_{km}\right]\right)}_{e_{k}}}^{2}$
$\displaystyle\hskip
36.135pt+\underbrace{\mathbb{E}\left[\left(\alpha_{k}\left(\sum_{m=1}^{M}b_{km}h_{km}x_{km}+n_{k}\right)-\mathbb{E}\left[\alpha_{k}\left(\sum_{m=1}^{M}b_{km}h_{km}x_{km}+n_{k}\right)\right]\right)^{2}\right]}_{\nu_{k}^{2}},$
(61)

where the estimator is
$\hat{G}=\alpha_{k}\left(\sum_{m=1}^{M}b_{km}h_{km}x_{km}+n_{k}\right)$ for
$k=1,2,\dots,K$. In (61), $\alpha_{k},b_{km},h_{km}$ and $n_{k}$ respectively
denote the correction factor for recovering the $k$th dimension of the true
gradient, the power allocation for the $k$th dimension of the local gradient
$x_{km}$ in the $m$th transmitter, the channel fading coefficient of the $k$th
sub-channel between the $m$th transmitter and the receiver, and the thermal
additive noise for the $k$th sub-channel. Further, $\nu_{k}^{2}$ and $e_{k}$
denote the estimator variance and bias, respectively. As apparent in (61), the
minimization of the MSE cost does not ensure $\hat{G}_{k}$ to be an unbiased
estimator of $G_{k}$. To resolve this issue, we formulate the unbiased
optimization problem as

$\displaystyle\underset{\\{\alpha_{k}\\},\\{b_{km}\\}}{\operatorname*{argmin}}~{}$
$\displaystyle\sum_{k=1}^{K}\nu_{k}^{2}(\alpha_{k},\\{b_{km}\\})$ (62) subject
to $\displaystyle
e_{k}(\alpha_{k},\\{b_{km}\\})=0,~{}~{}\sum_{k=1}^{K}b_{km}^{2}x_{km}^{2}\leq
E_{m},~{}~{}b_{km}\geq 0,~{}~{}\alpha_{k}\geq 0,$ $\displaystyle\forall
k=1,\dots,K;~{}~{}\forall m=1,\dots,M,$ (63)

where $E_{m}$ denotes the power budget of the $m$th transmitter. We note that
$E_{m}$, $\\{x_{1m},x_{2m},\dots,x_{Km}\\}$ and
$\\{b_{1m},b_{2m},\dots,b_{Km}\\}$ are only available to the $m$th transmitter
and $\alpha_{k},\forall k,$ are only available to the receiver. The
optimization problem defined in (60)-(61) can then be decomposed into two
stages. In the first stage, each transmitter $m$ utilizes all of its available
power and the local gradient information to compute the optimal power
allocation $\\{b_{1m},b_{2m},\dots,b_{Km}\\}$. In the second stage, the
receiver solves a consecutive optimization problem for finding the optimal
$\alpha_{k}$ for all $k=1,\dots,K$. The optimization problem at each
transmitter is formulated as

$\displaystyle\underset{\\{b_{km}\\}_{k=1:K}}{\operatorname*{argmax}}~{}$
$\displaystyle\zeta_{m}$ (64) subject to
$\displaystyle\mathbb{E}\left[\left(\zeta_{m}x_{km}-b_{km}h_{km}x_{km}\right)^{2}\right]=0,~{}~{}~{}\sum_{k=1}^{K}b_{km}^{2}x_{km}^{2}\leq
E_{m},~{}~{}~{}b_{km}\geq 0,$ $\displaystyle\forall k=1,\dots,K.$ (65)

The first constraint in (65) ensures that there is no additive bias in the
transmitted signal (i.e., the bias can be removed by a multiplicative factor),
while the second constraint is the power constraint. The first constraint can
be restated in a simpler form as $\zeta_{m}=b_{km}h_{km}$, and then the
solution can simply be obtained via the KKT conditions as

Lagrangian:
$\displaystyle\mathcal{L}(\\{b_{km}\\},\\{\lambda_{k}\\},\vartheta,\\{\beta_{k}\\})=\zeta_{m}-\sum_{k=1}^{K}\lambda_{k}(b_{km}h_{km}-\zeta_{m})-\vartheta\left(\sum_{k=1}^{K}b_{km}^{2}x_{km}^{2}-E_{m}\right)+\sum_{k=1}^{K}\beta_{k}b_{km}$
Stationarity: $\displaystyle\frac{\partial\mathcal{L}}{\partial
b_{km}}=0-\lambda_{k}h_{km}-2\vartheta x_{km}^{2}b_{km}+\beta_{k}=0\rightarrow
b_{km}=\frac{\beta_{k}-\lambda_{k}h_{km}}{2\vartheta x_{km}^{2}}$ Primal
Feasibility: $\displaystyle
b_{km}h_{km}=\zeta_{m},~{}~{}\sum_{k=1}^{K}b_{km}^{2}x_{km}^{2}\leq
E_{m},~{}~{}b_{km}\geq 0,~{}~{}\forall k=1\dots,K,$ Dual Feasibility:
$\displaystyle\vartheta\geq 0,~{}~{}\beta_{k}\geq 0,~{}~{}\forall
k=1,\dots,K,$ Comp. Slackness:
$\displaystyle\vartheta\left(\sum_{k=1}^{K}b_{km}^{2}x_{km}^{2}-E_{m}\right)=0,~{}~{}\sum_{k=1}^{K}\beta_{k}b_{km}=0.$

Then, the corresponding solution is given by

$\displaystyle\zeta_{m}^{*}=\sqrt{\frac{E_{m}}{\sum_{k=1}^{K}\frac{x_{km}^{2}}{h_{km}^{2}}}},$
$\displaystyle b_{km}^{*}=\frac{\zeta_{m}}{h_{km}},~{}~{}~{}\forall
k=1,\dots,K.$ (66)

(66) illustrates that the $m$th transmitter utilizes all of its power budget
to amplify its transmitted local gradient signal. Then, the corresponding
received signal is equivalent to the $\zeta_{m}^{*}$ times of the local
gradient signal, inducing a multiplicative bias. Yet, this bias can be removed
by multiplying the received signal with $\bm{\alpha}$ in the receiver.
However, the received signal at the receiver is a superposition of all
transmitted signals because of the over-the-air transmission. Therefore, a
single vector of optimal $\bm{\alpha}$ must be computed for removing the bias
and minimizing the estimator variance. To this end, in the second stage, the
receiver solves the following optimization problem

$\displaystyle\underset{\\{\alpha_{k}\\},}{\operatorname*{argmin}}~{}$
$\displaystyle\sum_{k=1}^{K}\nu_{k}^{2}(\alpha_{k},\\{b_{km}^{*}\\})$ (67)
subject to $\displaystyle
e_{k}(\alpha_{k},\\{b_{km}^{*}\\})=0,~{}~{}\alpha_{k}\geq 0,$
$\displaystyle\forall k=1,\dots,K.$ (68)

Next, we derive $e_{k}(\alpha_{k},\\{b_{km}^{*}\\})$ and
$\nu_{k}^{2}(\alpha_{k},\\{b_{km}^{*}\\})$. For ease of exposition, the
cumbersome steps are omitted here:

$\displaystyle e_{k}(\alpha_{k},\\{b_{km}^{*}\\})$
$\displaystyle=\mathbb{E}\left[\alpha_{k}y_{k}-\frac{1}{M}\sum_{m=1}^{M}x_{km}\right]=\mathbb{E}\left[\alpha_{k}\left(\sum_{m=1}^{M}h_{km}b_{km}^{*}x_{km}+n_{k}\right)-\frac{1}{M}\sum_{m=1}^{M}x_{km}\right]$
$\displaystyle=\alpha_{k}\left(\sum_{m=1}^{M}h_{km}b_{km}^{*}x_{km}\right)-\frac{1}{M}\sum_{m=1}^{M}x_{km}=\alpha_{k}\left(\sum_{m=1}^{M}\zeta_{m}^{*}x_{km}\right)-\frac{1}{M}\sum_{m=1}^{M}x_{km}$
$\displaystyle\nu_{k}^{2}(\alpha_{k},\\{b_{km}^{*}\\})$
$\displaystyle=\mathbb{E}\left[\left(\alpha_{k}\left(\sum_{m=1}^{M}h_{km}b_{km}^{*}x_{km}+n_{k}\right)-\mathbb{E}\left[\alpha_{k}\left(\sum_{m=1}^{M}h_{km}b_{km}^{*}x_{km}+n_{k}\right)\right]\right)^{2}\right]$
(69)
$\displaystyle=\mathbb{E}\left[\left(\alpha_{k}\left(\sum_{m=1}^{M}h_{km}b_{km}^{*}x_{km}-h_{km}b_{km}^{*}x_{km}+n_{k}\right)\right)^{2}\right]=\alpha_{k}^{2}\sigma^{2},$
(70)

where the expectation is taken with respect to $n_{k}$ for the realizations of
all $x_{km}$ and $\sigma^{2}=\mathbb{E}[n_{k}^{2}]$. By solving the KKT
conditions, we obtain the following solution

$\displaystyle\alpha_{k}^{*}$
$\displaystyle=\frac{\frac{1}{M}\sum_{m=1}^{M}x_{km}}{\sum_{m=1}^{M}h_{km}b_{km}^{*}x_{km}}=\frac{\frac{1}{M}\sum_{m=1}^{M}x_{km}}{\sum_{m=1}^{M}\zeta_{m}^{*}x_{km}}.$
(71)

From the implementation point of view, it is sensible to set
$\alpha_{k}^{\dagger}\simeq\frac{1}{\sum_{m=1}^{M}\zeta_{m}^{*}}$ since
$\\{x_{km}\\}$ is not available at the receiver. We herein also notice that
$\zeta_{m}^{*}$, for all $m=1,\dots,M$, are not available at the receiver as
well. Luckily, $\alpha_{k}^{\dagger}$ is a function of
$\sum_{m=1}^{M}\zeta_{m}^{*}$, and hence a subchannel could be allocated for
the over-the-air transmission of $\sum_{m=1}^{M}\zeta_{m}^{*}$. Subsequently,
it follows that the bias is given by

$\displaystyle e_{k}(\alpha_{k}^{*},\\{b_{km}^{*}\\})$
$\displaystyle=\mathbb{E}\left[\frac{1}{\sum_{m=1}^{M}\zeta_{m}^{*}}y_{k}-\frac{1}{M}\sum_{m=1}^{M}x_{km}\right]=\mathbb{E}\left[\frac{1}{\sum_{m=1}^{M}\zeta_{m}^{*}}\left(\sum_{m=1}^{M}h_{km}b_{km}^{*}x_{km}+n_{k}\right)-\frac{1}{M}\sum_{m=1}^{M}x_{km}\right]$
$\displaystyle=\mathbb{E}$
$\displaystyle\left[\frac{1}{\sum_{m=1}^{M}\zeta_{m}^{*}}\left(\sum_{m=1}^{M}\zeta_{m}^{*}x_{km}+n_{k}\right)-\frac{1}{M}\sum_{m=1}^{M}x_{km}\right]=\mathbb{E}\left[\frac{\sum_{m=1}^{M}\zeta_{m}^{*}x_{km}+n_{k}}{\sum_{m=1}^{M}\zeta_{m}^{*}}-\frac{1}{M}\sum_{m=1}^{M}x_{km}\right].$
(72)

Assuming that the distributions of $\\{x_{km}\\}$, for all
$k=1,\dots,K;m=1,\dots,M$, are identical across the subchannels and users, and
so are $\\{h_{km}\\}$, $\zeta_{m}^{*}$ can be simplified. For ease of
exposition, we denote $\mathbb{E}[x_{km}^{2}]=\varphi^{2}+\bar{x}^{2}$ and
$\mathbb{E}\left[\frac{1}{h_{km}^{2}}\right]=\varpi^{2}$. When the number of
subchannels, $K$, is large, we have that

$\displaystyle\zeta_{m}^{*}$
$\displaystyle=\frac{\sqrt{E_{m}}}{\sqrt{\sum_{k=1}^{K}\frac{x_{km}^{2}}{h_{km}^{2}}}}\underset{\text{when
$K$ is
large}}{\Longrightarrow}\zeta_{m}^{*}=\frac{\sqrt{E_{m}}}{\sqrt{K\mathbb{E}[x_{km}^{2}]\mathbb{E}\left[\frac{1}{h_{km}^{2}}\right]}}=\frac{\sqrt{E_{m}}}{\sqrt{K(\varphi^{2}+\bar{x}^{2})\varpi^{2}}}$
(73) $\displaystyle e_{k}(\alpha_{k}^{*},\\{b_{km}^{*}\\})$
$\displaystyle=\mathbb{E}\left[\frac{\sum_{m=1}^{M}\zeta_{m}^{*}x_{km}+n_{k}}{\sum_{m=1}^{M}\zeta_{m}^{*}}-\frac{1}{M}\sum_{m=1}^{M}x_{km}\right]=\sum_{m=1}^{M}\left[\frac{\sqrt{E_{m}}}{\sum_{m=1}^{M}\sqrt{E_{m}}}-\frac{1}{M}\right]x_{km}.$
(74)

Then, the norm $\lVert\mathbb{E}_{t}[\epsilon_{t}]\rVert_{2}^{2}$ in Theorem 1
can be expressed in terms of $e_{k},\forall k,$ as
$\lVert\mathbb{E}_{t}[\epsilon_{t}]\rVert_{2}^{2}=\sum_{k=1}^{K}\mathbb{E}_{t}[\epsilon_{t}]_{k}^{2}=\sum_{k=1}^{K}e_{k}^{2}$.
Similarly, the variance $\nu_{k}^{2}$ can also be computed as

$\displaystyle\nu_{k}^{2}\left(\alpha_{k}^{*},\\{b_{km}^{*}\\}\right)$
$\displaystyle={\alpha_{k}^{*}}^{2}\sigma^{2}=\left(\frac{1}{\sum_{m=1}^{M}\sqrt{\frac{E_{m}}{K\varpi(\varphi^{2}+\bar{x}^{2})}}}\right)^{2}\sigma^{2}=\frac{K\varpi^{2}(\varphi^{2}+\bar{x}^{2})}{\left(\sum_{m=1}^{M}\sqrt{E_{m}}\right)^{2}}\sigma^{2}.$
(75)

Finally, the MSE cost in Theorem 1 can be written as
$\mathbb{E}_{t}[\lVert\epsilon_{t}\rVert_{2}^{2}]=\sum_{k=1}^{K}({\nu}_{k}^{2}+e_{k}^{2})$.

###  A-E Alternative Formulations and Baselines

A Receiver Centric Approach. In what follows, we take a receiver-centric
approach by selecting a fixed estimator, $\alpha_{k}=\frac{1}{Mp},\forall k$,
at the receiver. Given the fixed estimator, the MSE objective function is set
as

$\begin{split}\mbox{{MSE}}(1/(Mp),\bm{b})&=\sum_{k=1}^{K}\sum_{m=1}^{M}\big{(}\frac{b_{km}h_{km}}{Mp}-\frac{1}{M}\big{)}^{2}x_{km}^{2}\\\
&+\sum_{k=1}^{K}\sum_{m=1}^{M}\sum_{\begin{subarray}{c}m^{\prime}=1\\\
m^{\prime}\neq
m\end{subarray}}^{M}\big{(}\frac{b_{km}h_{km}}{Mp}-\frac{1}{M}\big{)}\times\big{(}\frac{b_{km^{\prime}}h_{km^{\prime}}}{Mp}-\frac{1}{M}\big{)}x_{km}x_{km^{\prime}}+\sigma^{2}\frac{K}{M^{2}p^{2}}.\end{split}$

We note that the first term can be decomposed for different devices and be
solved in a distributed manner. The second term is coupled across different
users, and the third term is only expressed in terms of the variable $p$. If
there were no power constraints, the solution would be $b_{km}h_{km}=p$ for a
very large $p$. With this intuition, it is sensible to assign most (if not
all) of the multipliers $(\frac{b_{km}h_{km}}{Mp}-\frac{1}{M})$ zero. That is
to say, $p$ should be selected so that it fits most users to provide enough
power for the transmission of all the dimensions, while every user
individually solves the power allocation problem using the first term for a
given value of $p$, i.e.,

$\displaystyle\textbf{P2-1:\ }\min_{\\{b_{km}\\}}$
$\displaystyle\sum_{k=1}^{K}\big{(}\frac{b_{km}h_{km}}{Mp}-\frac{1}{M}\big{)}^{2}x_{km}^{2}$
s.t. $\displaystyle\sum_{k=1}^{K}b_{km}^{2}x_{km}^{2}\leq E.$

The optimal solution $\\{b_{km}^{*}\\}$ to the problem (P2-1) can be obtained
by finding $\lambda^{*}\geq 0$ satisfying the following two inequalities:

$\sum_{k=1}^{K}(b_{km}^{*})^{2}x_{km}^{2}\leq E,\quad
b^{*}_{km}=\frac{x_{km}^{2}p}{h_{km}(x_{km}^{2}-M^{2}p^{2}\lambda^{*}_{m})}.$
(76)

When $\lambda^{*}=0$, then the solution is $b_{km}^{*}=p/h_{km}$, the
objective function is minimized to 0 and the power constraint is satisfied
with strict inequality. On contrary, when $\lambda^{*}>0$, then the power
constraint is satisfied with equality and the optimal solution $b_{km}^{*}$
deviates from $p/h_{km}$. Then, the optimal solution $\\{b_{km}^{*}\\}$ to the
problem (P2-1) can be obtained by finding $\lambda^{*}>0$ satisfying the
following two equalities:

$\sum_{k=1}^{K}(b_{km}^{*})^{2}x_{km}^{2}=E,\quad
b^{*}_{km}=\frac{x_{km}^{2}p}{h_{km}(x_{km}^{2}-M^{2}p^{2}\lambda^{*}_{m})}.$
(77)

An equal power allocation approach. For comparison, we also consider the equal
power allocation scheme, in which equal power is allocated to each dimension,
i.e., we set $b_{km}=b_{m}$. Enforcing the power constraint of the devices,
e.g., $\sum_{k=1}^{K}b_{km}^{2}x_{km}^{2}=E$, leads to the equal power
solution that can be written as
$b_{m}^{*}=\sqrt{E/{\sum_{k=1}^{K}x^{2}_{km}}}$. Therefore, each device
applies the power $b_{m}^{*}$ in each dimension, taking the advantage of the
distribution of the data which is independent and identical.

An Alternative Formulation to P3.

When $K$ is large, the channel coefficients $\\{h_{km}\\}$ and the local
gradients $\\{x_{km}\\}$ are independent and identically distributed (i.i.d.)
across all users and subchannels. Hence, the optimization problem P3 can be
further simplified via the law of large numbers as

$\displaystyle\textbf{P4: }\underset{\tilde{b}}{\min}$
$\displaystyle~{}\mathbb{E}\left[\frac{1}{\tilde{b}_{km}\tilde{h}_{km}+\tilde{x}_{km}}\right]$
s.t.
$\displaystyle~{}\mathbb{E}[\tilde{b}_{km}/\tilde{x}_{km}]\leq\bar{E}_{m},$
$\displaystyle\tilde{b}_{km}\geq 0.$ (78)

We here notice that the general law of large numbers only applies if all
$\tilde{b}_{km}$ are identically distributed. By further noting that
$\tilde{b}_{km}$ depends on $\tilde{x}_{km}$ and $\tilde{h}_{km}$, we conclude
that the objective function for any $k$ is identically distributed in (A-E) .
Then, (A-E) simplifies as

$\displaystyle\textbf{P4-1: }\underset{\tilde{b}_{km}}{\min}$
$\displaystyle~{}\mathbb{E}\left[\frac{1}{\tilde{b}_{km}\tilde{h}_{km}+\tilde{x}_{km}}\right]$
s.t.
$\displaystyle~{}\mathbb{E}\left[\frac{\tilde{b}_{km}}{\tilde{x}_{km}}\right]\leq\bar{E}_{m},$
$\displaystyle\tilde{b}_{km}\geq 0.$ (79)

Subsequently, the KKT conditions are written as

Lagrangian:
$\displaystyle\mathcal{L}=\mathbb{E}\left[\frac{1}{\tilde{b}_{km}\tilde{h}_{km}+\tilde{x}_{km}}\right]+\lambda_{m}\left(\mathbb{E}\left[\frac{\tilde{b}_{km}}{\tilde{x}_{km}}\right]-\bar{E}_{m}\right)-\beta\tilde{b}_{km}$
(80) Stationarity:
$\displaystyle\frac{\partial\mathcal{L}}{\partial\tilde{b}_{km}}=\mathbb{E}\left[\frac{\lambda_{m}}{\tilde{x}_{km}}-\frac{\tilde{h}_{km}}{(\tilde{b}_{km}\tilde{h}_{km}+\tilde{x}_{km})^{2}}\right]=0\Leftarrow\frac{\lambda_{m}^{*}}{\tilde{x}_{km}}=\frac{\tilde{h}_{km}}{(\tilde{b}_{km}\tilde{h}_{km}+\tilde{x}_{km})^{2}}$
(81) Primal Feasibility:
$\displaystyle\mathbb{E}\left[\frac{\tilde{b}_{km}}{\tilde{x}_{km}}\right]\leq\bar{E}_{m},~{}~{}~{}\tilde{b}_{km}\geq
0,$ (82) Dual Feasibility: $\displaystyle\lambda_{m}\geq 0,~{}~{}~{}\beta\geq
0,$ (83) Comp. Slackness:
$\displaystyle\lambda_{m}\left(\mathbb{E}\left[\frac{\tilde{b}_{km}}{\tilde{x}_{km}}\right]-\bar{E}_{m}\right)=0,~{}~{}~{}-\beta\tilde{b}_{km}=0.$
(84)

Then, the solution is computed as follows

$\displaystyle\tilde{b}_{km}=$
$\displaystyle\left[\sqrt{\frac{\tilde{x}_{km}}{\lambda_{m}^{*}\tilde{h}_{km}}}-\frac{\tilde{x}_{km}}{\tilde{h}_{km}}\right]^{+}\Rightarrow
b_{km}=\sqrt{\left[\frac{\sigma}{h_{km}|x_{km}|\sqrt{\lambda_{m}^{*}}}-\frac{\sigma^{2}}{x_{km}^{2}h_{km}^{2}}\right]^{+}}$
(85) $\displaystyle\lambda_{m}^{*}<$
$\displaystyle\frac{\tilde{h}_{km}}{\tilde{x}_{km}}\Rightarrow\tilde{b}_{km}>0\text{
and }\lambda_{m}^{*}<\frac{h_{km}^{2}x_{km}^{2}}{\sigma^{2}}\Rightarrow
b_{km}>0$ (86)
$\displaystyle\mathbb{E}\left[\frac{\tilde{b}_{km}}{\tilde{x}_{km}}\right]=$
$\displaystyle\int_{0}^{\infty}\int_{0}^{\infty}\frac{\tilde{b}_{km}}{\tilde{x}_{km}}p(\tilde{h}_{km})p(\tilde{x}_{km})d\tilde{h}_{km}d\tilde{x}_{km}$
(87) $\displaystyle=$
$\displaystyle\int_{0}^{\infty}\int_{0}^{\lambda_{m}\tilde{x}_{km}}\left(\sqrt{\frac{1}{\lambda_{m}^{*}\tilde{h}_{km}\tilde{x}_{km}}}-\frac{1}{\tilde{h}_{km}}\right)p(\tilde{h}_{km})p(\tilde{x}_{km})d\tilde{h}_{km}d\tilde{x}_{km}=\bar{E}_{m}$
(88) $\displaystyle\mathbb{E}\left[{b}_{km}^{2}{x}_{km}^{2}\right]=$
$\displaystyle\int_{-\infty}^{\infty}\int_{0}^{\infty}\left[\frac{\sigma}{h_{km}|x_{km}|\sqrt{\lambda_{m}^{*}}}-\frac{\sigma^{2}}{x_{km}^{2}h_{km}^{2}}\right]^{+}x_{km}^{2}p(h_{km})p(x_{km})dh_{km}dx_{km}$
(89) $\displaystyle=$
$\displaystyle\int_{-\infty}^{\infty}\int_{\sqrt{\frac{\sigma^{2}\lambda_{m}^{*}}{x_{km}^{2}}}}^{\infty}\left(\frac{|x_{km}|\sigma}{h_{km}\sqrt{\lambda_{m}^{*}}}-\frac{\sigma^{2}}{h_{km}^{2}}\right)p(h_{km})p(x_{km})dh_{km}dx_{km}=\bar{E}_{m}.$
(90)

By using the results above, the bias term, $e_{k}$, is derived as

$\displaystyle e_{k}(\alpha_{k}^{*},\\{b_{km}\\})=$
$\displaystyle\mathbb{E}\left[\alpha_{k}^{*}y_{k}\right]-\frac{1}{M}\sum_{m=1}^{M}x_{km}=\mathbb{E}\left[\sum_{m=1}^{M}\left(\alpha_{k}^{*}b_{km}h_{km}-\frac{1}{M}\right)x_{km}+\alpha_{k}^{*}n_{k}\right]$
(91) $\displaystyle=$
$\displaystyle\mathbb{E}\left[\sum_{m=1}^{M}\left(\frac{1}{M}\sqrt{\left[\frac{\sigma}{|x_{km}|h_{km}\sqrt{\lambda_{m}^{*}}}-\frac{\sigma^{2}}{x_{km}^{2}h_{km}^{2}}\right]^{+}}h_{km}-\frac{1}{M}\right)x_{km}+\frac{n_{k}}{M}\right]$
(92) $\displaystyle=$
$\displaystyle\frac{1}{M}\sum_{m\in\mathcal{S}^{\dagger}}\left(\sqrt{\frac{\sigma
h_{km}}{|x_{km}|\sqrt{\lambda_{m}^{*}}}-\frac{\sigma^{2}}{x_{km}^{2}}}-1\right)x_{km},$
(93)

where $\mathcal{S}^{\dagger}$ represents the set of transmitters for which
(86) is satisfied. To compute the MSE cost, we compute the variance
${\nu}^{2}(\alpha_{k}^{*},\\{b_{km}\\})$ by taking the expectation with
respect to $x_{km}$

$\displaystyle{\nu}_{k}^{2}(\alpha_{k}^{*},\\{b_{km}\\})=$
$\displaystyle\mathbb{E}\left[\left(\sum_{m=1}^{M}\alpha_{k}^{*}b_{km}h_{km}x_{km}+\alpha_{k}^{*}n_{k}-\mathbb{E}\left[\sum_{m=1}^{M}\alpha_{k}^{*}b_{km}h_{km}x_{km}\right]\right)^{2}\right]$
(94)
$\displaystyle=\mathbb{E}\left[\left(\sum_{m=1}^{M}\frac{1}{M}\left\\{\sqrt{\left[\frac{\sigma
h_{km}}{\lvert
x_{km}\rvert\sqrt{\lambda_{m}^{*}}}-\frac{\sigma^{2}}{x_{km}^{2}}\right]^{+}}x_{km}-\sqrt{\left[\frac{\sigma
h_{km}}{\lvert
x_{km}\rvert\sqrt{\lambda_{m}^{*}}}-\frac{\sigma^{2}}{x_{km}^{2}}\right]^{+}}x_{km}\right\\}+\frac{n_{k}}{M}\right)^{2}\right]=\frac{\sigma^{2}}{M^{2}}.$
(95)

Consequently, the cost
$\mathbb{E}_{t}\left[\lVert\epsilon_{t}\rVert_{2}^{2}\right]$ in Theorem 1 is
computed as
$\mathbb{E}_{t}[\lVert\epsilon_{t}\rVert_{2}^{2}]=\sum_{k=1}^{K}(e_{k}^{2}+{{\nu}}_{k}^{2})$.