text:
\section{Introduction}
The availability of massive data sets for machine learning has led to the development of distributed computing platforms that can improve the computational efficiency of training tasks by partitioning a global data set among multiple workers \cite{stoica2017berkeley, li2014scaling,abadi2016tensorflow,dhar2019device,zhang2019deep,celik2018wireless}. Decentralization may come at the cost of a possible degradation in accuracy due to the reliance on local optimization steps \cite{ben2019demystifying}. Moreover, when workers are connected to the central server via capacity-limited channels, communication-related impairments can cause additional performance degradation \cite{zhu2019broadband,amiri2020machine,yang2020federated,lan2020capacity,ren2020scheduling, yang2019scheduling,liu2020privacy, cao2020optimized,zhang2020gradient,celik2018wireless}.

As we will review,
existing works on wireless distributed learning focus on a standard \emph{frequentist} formulation  of the learning problem, whereby the goal is to identify a single vector of model parameters by minimizing (a function of) the training loss  \cite{amiri2020machine, liu2020privacy, cao2020optimized, zhu2020one, Xing2020Federated}. Much research has been specifically carried out on studying the advantages of \emph{uncoded transmission} to enable ``over-the-air'' computing \cite{amiri2020machine, zhu2019broadband, liu2020privacy, cao2020optimized}. Frequentist learning is known to be limited in its ability to quantify epistemic uncertainty, i.e., uncertainty in the model parameter space, yielding overconfident decisions \cite{guo2017calibration,lakshminarayanan2016simple}.  In contrast, \emph{Bayesian} learning provides a principled way to obtain predictive models that account for epistemic uncertainty. This is done by optimizing not over individual model parameter vectors, but over distributions in the model parameter space \cite{barber2012bayesian,simeone2017brief}. Other advantages of Bayesian learning over frequentist learning include the capacity to incorporate prior knowledge, to guide active learning, and to enable continual learning \cite{barber2012bayesian}. Practical implementations of Bayesian learning involve either approximating the posterior distribution over the model parameters -- an approach known as \emph{variational inference (VI)}, or producing approximate samples from the posterior -- an approach known as \emph{Monte Carlo (MC) sampling} \cite{angelino2016patterns}.

This paper considers, for the first time, distributed Bayesian learning in wireless networks. We specifically focus on a wireless data center setting as illustrated in Fig. \ref{Fig: sys} \cite{celik2018wireless}, and on the standard class of ``embarrassingly parallel'' \emph{consensus Monte Carlo (CMC)} protocols \cite{scott2016bayes, rabinovich2015variational}. In wireless data centers, the global data set is available at the central server, which distributes subsets of data points to the workers. Wireless data centers have been introduced as an alternative to conventional wired topologies to reduce cost and to cope with unbalanced traffic \cite{celik2018wireless}. In CMC protocols, first, each worker produces samples from a local posterior based on the locally available data; then, the server carries out a one-shot aggregation of the received messages in order to approximate samples from the global posterior distribution.

As a key contribution, we introduce the novel idea of \emph{channel-driven MC sampling}. We specifically demonstrate that, under uncoded wireless transmission, channel noise can be utilized as an integral component of MC sampling. As a result, if properly accounted for, channel noise may not cause harm to the performance of MC-based distributed Bayesian learning. This is in stark contrast to wireless distributed learning based on frequentist principles, for which channel noise is generally detrimental to the learning performance. {\color{black}An exception to this statement was demonstrated in \cite{sery2021over}, which argues that the presence of noise can be useful when training using non-convex objectives.}

\subsection{Wireless Distributed Learning}
Most work on wireless distributed learning focuses on the ``federated'' setting, in which data originates at the workers and it remains unknown to the server. In contrast, in a wireless data center system, the server has access to the global data set, and it distributes data to the workers to benefit from computational parallelism. In both cases, distributed learning protocols can be embarrassingly parallel, that is, one-shot, or iterative. In the former case, which we focus on, the workers carry out local computations and communicate only once with the server, which aggregates the results of the local computations \cite{neiswanger2013asymptotically, scott2016bayes,rabinovich2015variational}; while, in the latter, local computing and communication steps are alternated. To counteract hostile channel fading and additive channel noise conditions in wireless distributed learning, various strategies have been proposed including optimal power control \cite{liu2020privacy, cao2020optimized,zhang2020gradient} and over-the-air computing \cite{zhu2019broadband, amiri2020machine, yang2020federated}.  Over-the-Air computing utilizes non-orthogonal multi-access (NOMA) and uncoded transmission to enable the aggregation of signals transmitted by the workers ``over-the-air", leveraging the waveform-superposition property of a multi-access channel. For the wireless data center setting, existing works address issues such as optimizing the number of workers to balance the tradeoff between communication latency and computation time \cite{song2020wireless}, and designing data partition strategy to cope with straggling workers \cite{lyu2019optimal}.

\subsection{Distributed Bayesian Learning}

To the best of our knowledge, as reviewed above, prior work on distributed wireless learning -- in both federated and wireless data center settings -- is based on frequentist learning. As mentioned, since exact Bayesian learning is generally intractable, existing Bayesian learning schemes are based on approximate methods leveraging either VI or MC sampling. VI assumes a parametric family for the posterior distribution of the model parameters, and its performance is generally limited by the bias due to the variational approximation \cite{angelino2016patterns}. In contrast, MC sampling aims at obtaining samples from the posterior distribution, and is asymptotically exact in the number of samples. Given their reliance on parametric optimization, distributed VI schemes can be implemented in a manner similar to existing frequentist protocols in federated learning, including local stochastic gradient descent (SGD) optimization and global aggregation \cite{corinzia2019variational,kassab2020federated}.

Distributed MC sampling can also be implemented either in a one-shot fashion or an iterative manner. Most common is the former class of protocols, which are known as CMC \cite{scott2016bayes}. The key problem in the design of CMC schemes is how to aggregate the local samples provided by the workers \cite{neiswanger2013asymptotically}. In this regard, by approximating the posterior distribution as Gaussian, one can obtain optimal weights in closed form, which can be estimated using the local samples \cite{scott2016bayes}. For general distributions, reference \cite{rabinovich2015variational} proposed to design the aggregation function via a VI formulation. Less well studied are iterative protocols, which include distributed stochastic gradient Markov Chain MC strategies \cite{ahn2014distributed}. These protocols are inherently sequential, making parallel implementations challenging.

All prior works on CMC assume ideal communication.
Nevertheless, the existence of  channel noise in wireless data centers can severely distort the distribution of the received samples. This observation motivates us to seek optimized aggregation strategies that account for the impairments caused by wireless transmission. We refer to the resulting setting as \emph{wireless CMC}.

\subsection{Contributions, Organization, and Notations}

In this paper, we study for the first time the design of wireless CMC protocols for distributed one-shot Bayesian learning in a wireless data center setting. As illustrated in Fig. \ref{Fig: sys}, the global data set is partitioned into disjoint subsets for the distributed workers to enable the implementation of CMC. After local sampling, each worker uploads the produced local samples  to the server via uncoded transmission via either orthogonal or non-orthogonal multi-accesses protocols.  At the server, the received signals are aggregated by using an optimized function to produce approximated samples from the global posterior. The main findings and contributions of the paper can be summarized as follows.

\noindent $\bullet$ {\bf Wireless Gaussian Consensus Monte Carlo (WGCMC): } To start, we study as in \cite{scott2016bayes} the baseline case in which the local posteriors follow Gaussian distributions, and we introduce a wireless CMC protocol that is asymptotically optimal with respect to the number of produced samples.  The solution automatically adapts to the signal-to-noise ratio (SNR) level, and, under Gaussian local posteriors, it is demonstrated to leverage channel noise as a component of the sample generation process,  instead of a nuisance to be mitigated. We consider orthogonal multiple access (OMA), as well as NOMA  with over-the-air computing to enhance spectral efficiency.

\noindent $\bullet$ {\bf Wireless Variational Consensus Monte Carlo (WVCMC):}  We then target  the general case with arbitrary  local posteriors. We propose a variational design methodology that aims at reducing the Kullback-Liebler (KL) divergence of the distribution of the aggregated samples with respect to the true global posterior.  The proposed principled solution approach is based on SGD steps at the server.

\noindent $\bullet$ {\bf Experiments:}
 We provide extensive numerical results on both synthetic and real data sets to demonstrate the advantages of channel-driven CMC.

{\bf Organization:}  The remainder of the paper is organized as follows.
Section~\ref{sec: system model} introduces the system model.
Section~\ref{sec: GCMC} introduces the necessary background on CMC.
Section~\ref{sec: tx rx} presents the general form of signal design for wireless CMC, while WGCMC is introduced in Section \ref{sec: WGCMC} and WVCMC in Section~\ref{sec: WVCMC}.
Section~\ref{sec: sim} provides numerical results, followed by conclusions in Section~\ref{sec: conclusions}.

{\bf Notations:} The following notations are used throughout the paper.  $\otimes$ denotes Kronecker product. ${\bf 1}_l$ is all-ones $l \times 1$ vector. $\delta[\cdot]$ is the Dirac delta function. ${\mathbf{A}}^{-1/2}={\mathbf{U}} {{\mathbf{m}} \Lambda}^{-1/2} {\mathbf{U}}^{\sf T}$ is the operation for any positive semidefinite matrix ${\mathbf{A}}={\mathbf{U}} {{\mathbf{m}} \Lambda} {\mathbf{U}}^{\sf T}$, where ${\mathbf{U}}$ is comprised by singular vectors and   $ {{\mathbf{m}} \Lambda} $ represents the non-zero diagonalized singular values.  $[{\mathbf{A}}]^+$ sets the negative singular values of ${\mathbf{A}}$ to 0, and ${\mathcal{I}}(\cdot)$ is the indicator function.

\section{System Model} \label{sec: system model}
\begin{figure}[t]
\centering
\includegraphics[width=8 cm]{./System.pdf}
\caption{Bayesian learning via consensus Monte Carlo (CMC), i.e., embarrassingly parallel distributed MC sampling, in a wireless data center.}
\label{Fig: sys}
\end{figure}

\subsection{Setting}
As shown in Fig. \ref{Fig: sys}, we consider a wireless data center setting \cite{celik2018wireless} including a server and $K$ workers connected to it via a shared wireless channel.  The global dataset $\mathcal{Z}$ at the server consists of $N$ data samples $\mathcal{Z}= \{{\mathbf{z}}_n\}_{n=1}^N$.  For supervised learning applications, each data sample ${\mathbf{z}}_n=({\mathbf{u}}_n, v_n)$ is in turn partitioned into a covariate vector ${\mathbf{u}}_n$ and a label $v_n$; while, for unsupervised learning applications such as generative modeling, it consists of a single vector ${\mathbf{z}}_n$ (see, e.g., \cite{simeone2017brief}). The goal of the system is to carry out Bayesian learning via Monte Carlo (MC) sampling. This entails generating samples that are (approximately) distributed according to posterior distribution $p({{\mathbf{m}} \vartheta} | \mathcal{Z} )$ of model parameter ${{\mathbf{m}} \vartheta} \in \mathbb{R}^d$ given data set $\mathcal{Z}$.
Samples from the posterior can be used to obtain ensemble predictive distributions that approximate the optimal Bayesian predictor (see, e.g., \cite{barber2012bayesian,murphy2012machine, simeone2017brief}).

{\color{black}Denote as $p({\mathbf{z}} | {{\mathbf{m}} \vartheta} )$ the likelihood of model parameter vector ${\mathbf{m}} \vartheta$ for data point ${\mathbf{z}}$, and as $p({{\mathbf{m}} \vartheta})$  the prior distribution of the model parameter vector.  Under the standard assumption of conditionally independent and identical (i.i.d.) data points, the global posterior distribution can be written as
\begin{align}\label{eq: global posterior}
\text{(Global posterior)}  \quad p({{\mathbf{m}} \vartheta} | \mathcal{Z} ) \propto p({{\mathbf{m}} \vartheta}) \prod_{n=1}^N p({\mathbf{z}}_n | {{\mathbf{m}} \vartheta} ).
\end{align}}
 Note that, for supervised learning problems, the likelihood can be written as $p({\mathbf{z}}|{{\mathbf{m}} \vartheta})=p({\mathbf{u}},v|{{\mathbf{m}} \vartheta})$ for generative models and as $p({\mathbf{z}}|{{\mathbf{m}} \vartheta})=p(v|{\mathbf{u}},{\mathbf{m}} \vartheta)$ for discriminative models \cite{simeone2017brief}.
We make the key assumption that exact sampling from the posterior $p({{\mathbf{m}} \vartheta} | \mathcal{Z} )$ is computationally intensive due to the size of the data set $\mathcal{Z}$. Possible solutions include  adaptive subsampling procedures \cite{bardenet2014adaptive} and efficient stochastic gradient-based algorithms such as stochastic gradient Langevin dynamics (SGLD) \cite{welling2011bayesian}, which require limited access to data at the server at any given time. This paper considers an alternative, and potentially more effective approach that leverages the presence of the $K$ distributed workers in the wireless data center system in Fig.~\ref{Fig: sys}.

We specifically propose, and study, the wireless implementation of distributed MC sampling solutions,  whereby workers perform parallel MC sampling on distinct subsets of the data \cite{angelino2016patterns}.
We focus on the class of embarrassingly parallel strategies, known as consensus Monte Carlo (CMC), in which communication from workers to server occurs only once, after local processing at the workers is completed.
In CMC, the global dataset $\mathcal{Z}$ is partitioned into $K$ non-overlapping subsets $\mathcal{Z}_1, \cdots, \mathcal{Z}_K$, each of size $N_k$, and each data set $\mathcal{Z}_k$ is allocated to the $k$-th worker. We have $\sum_{k=1}^K N_k = N$.  {\color{black}Data delivery is offline, and hence one can consider data to be available without distortion at the workers.}
 It is assumed that the $k$-th worker has sufficient computational resources to generate samples from the subposterior $p_k({{\mathbf{m}} \vartheta}|\mathcal{Z}_k)$. The subposterior depends on the allocated subset $\mathcal{Z}_k$ and on the underweighted prior $p({{\mathbf{m}} \vartheta}) ^{1/K}$ as
\begin{align}\label{eq: subposterior}
\text{(Subposterior)} \quad p_k({{\mathbf{m}} \vartheta}|\mathcal{Z}_k) \propto   p({{\mathbf{m}} \vartheta}) ^{1/K} p(\mathcal{Z}_k|{{\mathbf{m}} \vartheta}) ,
\end{align}
where we have defined the local likelihood function as $p(\mathcal{Z}_k|{{\mathbf{m}} \vartheta})=\prod_{{\mathbf{z}}\in\mathcal{Z}_k}p({\mathbf{z}}|{{\mathbf{m}}\vartheta})$. {\color{black}This assumption is justified by the fact that the complexity of MC sampling scales up with the size of the underlying data set \cite{scott2016bayes, rabinovich2015variational, angelino2016patterns} (see Appendix \ref{app: sampling complexity} for some additional discussion on this point).}
Note that the global posterior \eqref{eq: global posterior} can be expressed as  the product of subposteriors \eqref{eq: subposterior}, upon normalization,~as
\begin{align}
 p({{\mathbf{m}} \vartheta}|\mathcal{Z}) \propto    \prod_{k=1}^Kp_k({{\mathbf{m}} \vartheta}|\mathcal{Z}_k).  \label{eq: fact global posterior}
\end{align}

Each worker $k=1,\dots, K$ generates $S$ samples ${{\mathbf{m}} \vartheta}_k^{(s)}\sim  p_k({{\mathbf{m}} \vartheta}|\mathcal{Z}_k)$ for $s=1,\dots,S$ and transmits them one-by-one, as detailed in the next section, over the wireless channel to the server. Based on the received signals on the wireless channel, the server aims at obtaining samples ${{\mathbf{m}} \vartheta}^{(s)}$ for $s=1,\dots, S$ that are approximately drawn from the global posterior $p({{\mathbf{m}} \vartheta}|\mathcal{Z})$ by applying an aggregation function $F(\cdot)$. We note that the aggregation function $F(\cdot)$ may depend on the data set $\mathcal{Z}$, which is available at the server.

In existing CMC methods \cite{scott2016bayes, rabinovich2015variational}, the samples $\{{{\mathbf{m}} \vartheta}_k^{(s)}\}_{k=1}^K$ are assumed to be received noiselessly to the server, and the aggregation function is chosen as the linear combination
\begin{align}
{{\mathbf{m}} \vartheta}^{(s)}=F({{\mathbf{m}} \vartheta}_1^{(s)},\cdots,{{\mathbf{m}} \vartheta}_K^{(s)})=\sum_{k=1}^K {\mathbf{W}}_k {{\mathbf{m}} \vartheta}_k^{(s)} \label{eq: agg func}
\end{align}
for some $d \times d$ weight matrices ${\mathbf{W}}_1,\cdots, {\mathbf{W}}_k$. {\color{black}CMC does not require  feedback from the edge server or other communications beyond the transmission of the local posterior samples from workers to the server.}  The weight matrices are optimized based on the samples $\{\{{{\mathbf{m}} \vartheta}^{(s)}_k\}_{k=1}^K\}_{s=1}^S$ and on data set $\mathcal{Z}$ \cite{scott2016bayes, rabinovich2015variational}. CMC is successful when the distribution $q_F({{\mathbf{m}} \vartheta})$ of the samples ${{\mathbf{m}} \vartheta}^{(s)}$ computed at the server is a good approximation of the global posterior $p({{\mathbf{m}} \vartheta}|\mathcal{Z})$.

{\color{black}In this paper, we will consider for the first time the case in which the samples $\{\{{{\mathbf{m}} \vartheta}_k^{(s)}\}_{k=1}^K\}_{s=1}^S$  have to be communicated to the server over a wireless channel, and demonstrate the novel use of channel noise for Monte Carlo sampling. The proposed channel-driven sampling is enabled by the use of analog communication for the transmission of the samples $\{{{\mathbf{m}} \vartheta}^{(s)}_k\}_{k=1}^K$, whereby channel noise is directly added to the samples produced and transmitted by the workers. The CMC aggregation function \eqref{eq: agg func} can be implemented using both orthogonal and non-orthogonal multiple access schemes, which have distinct advantages and disadvantages.  With orthogonal multiple access (OMA) transmission, the local samples $\{{{\mathbf{m}} \vartheta}^{(s)}_k\}_{k=1}^K$ are received separately. Therefore, channel noise acts as a distinct  source of randomness for each sample, and arbitrary weight matrices $\{{\mathbf{W}}_k\}_{k=1}^K$ can be applied at the server. In contrast, under non-orthogonal multiple access (NOMA) transmission, AirComp can increase communication efficiency, but at the costs of sharing the channel noise across samples and of limiting the choice of the weight matrices $\{{\mathbf{W}}_k\}_{k=1}^K$ (see Sec. \ref{sec: tx rx}).
}

\subsection{Communication Model}\label{sec: comm model}
Following recent work on communication for machine learning \cite{amiri2020blind, zhu2020one, yang2020federated, zhu2019broadband,  zhu2020toward}, we assume that workers communicate with the server via a  general vector wireless communication model that may account for multi-carrier and/or single/multi-antenna setups. Accordingly, each worker transmits $m_t$  symbols, over time, space, and/or frequency domains, and the server receives $m_r$ symbols over a channel described by a real matrix and additive noise.  We focus on a real channel model with no loss of generality, since complex vector channels can be described in terms of real vector channel of twice the size of the original model \cite{liu2020privacy}. We will specifically consider two general classes of communication protocols, namely orthogonal multiple access (OMA) and non-orthogonal multiple access (NOMA).

\subsubsection{Orthogonal multiple access (OMA)} Under OMA,  workers transmit over orthogonal communication resources, and the $m_r \times 1$ received signal from the $k$-th worker is
\begin{align} \label{eq: rev signal}
{\mathbf{y}}_k={\mathbf{H}}_k{\mathbf{x}}_k+{\mathbf{n}}_k,
\end{align}
where ${\mathbf{H}}_k\in \mathbb{R}^{m_r\times m_t}$ is the channel matrix between worker $k$ and server; ${\mathbf{x}}_k\in \mathbb{R}^{m_t}$ is the signal transmitted by worker $k$; and ${\mathbf{n}}_k\in\mathbb{R}^{m_r}$ is channel noise, with entries being independent and identically distributed (i.i.d.) $\mathcal{N}(0, N_0)$ random variables (rvs).  Specifically, in each communication block, a given worker $k$ is scheduled to transmit a sample ${{\mathbf{m}} \vartheta}_k^{(s)}$. We focus on analog, i.e., uncoded transmission strategies, whereby ${{\mathbf{m}} \vartheta}_k^{(s)}$ is mapped through a linear transformation to ${\mathbf{x}}_k$, as we will detail in Sec. \ref{sec: tx rx}.

\subsubsection{Non-orthogonal multiple access (NOMA)} Under NOMA, all workers transmit simultaneously, and the $m_r \times 1$ received signal given by the noisy superposition
\begin{align} \label{eq: rev signal noma}
{\mathbf{y}}=\sum_{k=1}^K{\mathbf{H}}_k{\mathbf{x}}_k+{\mathbf{n}},
\end{align}
where ${\mathbf{H}}_k\in \mathbb{R}^{m_r\times m_t}$ is the channel matrix between worker $k$ and server; ${\mathbf{x}}_k\in \mathbb{R}^{m_t}$ is the signal transmitted by worker $k$; and ${\mathbf{n}}\in\mathbb{R}^{m_r}$ is channel noise, with entries being i.i.d. $\mathcal{N}(0, N_0)$ random rvs. As for OMA, each signal ${\mathbf{x}}_k$ encodes a sample ${{\mathbf{m}} \vartheta}_k^{(s)}$ produced by worker $k$ through a linear transformation.

In both cases, the channel matrices generally change from one communication block, given by \eqref{eq: rev signal} or \eqref{eq: rev signal noma}, to another. {\color{black}The channel matrix ${\mathbf{H}}_k$ and its statistics are assumed to be known at  worker $k$. This can be achieved via the broadcasting of pilots by the server in time-division duplex (TDD) systems, or via uplink pilot transmission and feedback from the server for  frequency-division duplex (FDD) systems.}  An important aspect of the model is that, in order to receive information from all devices, NOMA can operate with a fraction $1/K$ of the resources of OMA, since the server obtains a superposition of all $K$ signals in a single communication block. This property  has been successfully leveraged by AirComp strategies based on analog transmission \cite{zhu2018mimo, amiri2020machine, yang2020federated, zhu2020over, liu2020privacy}.

For both OMA and NOMA, the transmit per-block power constraint is set as the inequality
\begin{align}\label{eq: power constraint}
\mathsf{E}[\|{\mathbf{x}}_k\|^2]\leq P,
\end{align}
where the expectation is taken with respect to the distribution of the channels and of the encoded samples. The average over the channel implies consideration of a long-term power constraint.
{\color{black}This is sensible for the common setting in which the workers transmit multiple samples $S>1$ in order to enable posterior averaging at the server, e.g., for ensemble prediction (see Sec. \ref{sec: sim}). }
A short-term power constraint can be accounted for in a similar way.

Throughout the paper, we assume the conditions $m_t\geq m_r\geq d$, where we recall that $d$ is the dimension of the model vector ${\mathbf{m}} \vartheta$. {\color{black} As is common in the Bayesian learning literature, the value of $d$ is generally small in prior CMC works \cite{rabinovich2015variational,scott2016bayes}. For larger models, i.e., for larger $d$, the transmission of one sample would require multiple coherence blocks and/or use sufficiently large antenna arrays. This setting can be modelled by introducing block diagonal channel matrices, including the channel matrices for each block on the block diagonal.} The integration of analog compression techniques (see, e.g., \cite{abdi2020analog, abdi2020quantized}) with the wireless CMC strategies to be developed here can alleviate this problem and is left for future work.
Finally, we assume that the channel matrices $\{{\mathbf{H}}_k\}_{k=1}^K$ have full rank with probability one. This is a common assumption in many studies on multi-antenna wireless communication (see, e.g., \cite{jafar2011interference}).

\section{Background: \\ Gaussian Consensus Monte Carlo}\label{sec: GCMC}
In this section, we review the Gaussian CMC (GCMC) scheme introduced in \cite{scott2016bayes} that assumes noiseless, ideal communication. Accordingly, GCMC assumes that the local samples $\{{{\mathbf{m}} \vartheta_k^{(s)}}\}_{k=1}^K$, for $s=1,\dots,S$ are available at server. We will argue in Sec. \ref{sec: WGCMC} that a direct application of this scheme in the wireless setting under study is strictly suboptimal, even in the asymptotic regime of large number of samples, $S$, motivating the introduction of alternative schemes.

GCMC is based on the observation that, if each local subposterior $p_k({{\mathbf{m}} \vartheta}|\mathcal{Z}_k)$ is a  Gaussian $\mathcal{N}({{\mathbf{m}} \mu}_k,{\mathbf{C}}_k)$ distribution with invertible covariance matrix ${\mathbf{C}}_k$, the global posterior is the Gaussian $\mathcal{N}({{\mathbf{m}} \mu},{\mathbf{C}})$ distribution, where ${{\mathbf{m}} \mu}={\mathbf{C}} (\sum_{k=1}^K {\mathbf{C}}_k^{-1}{{\mathbf{m}} \mu}_k) $ and ${\mathbf{C}}=(\sum_{k=1}^K{\mathbf{C}}_k^{-1})^{-1}$ . If follows that, given local samples ${{\mathbf{m}} \vartheta}_k \sim \mathcal{N}({{\mathbf{m}} \mu}_k,{\mathbf{C}}_k)$, for $k=1,\dots,  K$, a global sample ${{\mathbf{m}} \vartheta}\sim \mathcal{N}({{\mathbf{m}} \mu},{\mathbf{C}})$ can be obtained by applying the linear combination \eqref{eq: agg func} with weight matrices
\begin{align}
{\mathbf{W}}_k=\Big(\sum_{k=1}^K{\mathbf{C}}_k^{-1}\Big)^{-1} {\mathbf{C}}_k^{-1}. \label{eq: w GCMC}
\end{align}

Since the matrices ${\mathbf{C}}_1,\cdots, {\mathbf{C}}_K$ are not known, GCMC estimates them by using the $S$ samples $\{{{\mathbf{m}} \vartheta}_k^{(s)}\}_{k=1}^K$ with $s=1,\dots,S$ received from the workers. In particular, given the $S$ samples $\{{{\mathbf{m}} \vartheta}_k^{(s)}\}_{s=1}^S$ received form worker $k$, covariance ${\mathbf{C}}_k$ is  estimated  as \begin{align}\label{eq: GCMC}
\widehat{{\mathbf{C}}}_k = \frac{1}{S-1} \sum_{s=1}^{S} \left({{\mathbf{m}} \vartheta}_k^{(s)} - \widehat{{\mathbf{m}} \mu}_k \right) \left({{\mathbf{m}} \vartheta}_k^{(s)} - \widehat{{\mathbf{m}} \mu}_k \right)^{\sf T},
\end{align}
with $\widehat{{\mathbf{m}} \mu}_k= \frac{1}{S} \sum_{s=1}^{S} {{\mathbf{m}} \vartheta}_k^{(s)} $. Then, the function \eqref{eq: agg func} is applied using the estimate $\widehat{{\mathbf{C}}}_k$ \eqref{eq: GCMC} in lieu of ${\mathbf{C}}_k$ in the weights~\eqref{eq: w GCMC}. Accordingly, GCMC computes each global sample ${{\mathbf{m}} \vartheta}^{(s)}$ as
\begin{align}\label{eq: est GCMC result}
\text{(GCMC)} \quad {{\mathbf{m}} \vartheta}^{(s)}= \sum_{k=1}^K  \underbrace{  \Big(  \sum_{k'=1}^K\widehat{{\mathbf{C}}}_{k'}^{-1}\Big)^{-1} \widehat{{\mathbf{C}}}_k^{-1}} _{{\mathbf{W}}_k} {{\mathbf{m}} \vartheta}_k^{(s)} .
\end{align}
Note that, as a result, the function \eqref{eq: est GCMC result} applied by the server is non-linear in the samples $\{\{{{\mathbf{m}} \vartheta}^{(s)}\}_{k=1}^K\}_{s=1}^S$.  It is also noted that GCMC does not require the use of the data set $\mathcal{Z}$ at the server. {\color{black}The matrix inversion in \eqref{eq: est GCMC result} incurs prohibitive computational complexity for high-dimensional models. One possible solution, as proposed in \cite{rabinovich2015variational,scott2016bayes}, is to restrict the weight matrices ${\mathbf{W}}_k$ to be diagonal, with each element corresponding to the inverse of the marginal variance of the samples.}

{\color{black}While GCMC is defined under the assumption of Gaussian subposteriors, reference \cite{scott2016bayes} argues that the choice \eqref{eq: w GCMC}--\eqref{eq: GCMC} for the combining matrices $\{{\mathbf{W}}_k\}_{k=1}^K$ is suitable for other settings as well. This is theoretically well justified for large data sets, i.e., when the data set sizes $\{N_k\}_{k=1}^K$ are large as compared to the model dimension, since the subposteriors in this case can be well approximated by Gaussian distributions according to Bernstein-von Mises theorem \cite{johnstone2010high}. In this paper, we will adopt throughout a linear aggregation as in \eqref{eq: agg func}, although the methods to be introduced in the next sections apply to more general parametric aggregation functions.}

\section{Wireless Consensus Monte Carlo: \\Signal Design}\label{sec: tx rx}
Unlike prior works on CMC \cite{rabinovich2015variational,scott2016bayes}, including the baseline GCMC scheme reviewed in the previous section, in the following sections, we propose solutions that account for the  wireless communication channel between workers and server.
 A key novel idea introduced by this work is that the noise added by the wireless channel
 may not cause any performance degradation as compared to an ideal scenario with noiseless transmission: Noise can contribute to MC sampling, and is not necessarily a nuisance to be mitigated. In this section,  we describe the general form of the considered uncoded transmission strategies for CMC. The next two sections propose specific wireless CMC protocols. Throughout, we focus on analog, or uncoded, transmission strategies with the aims of leveraging channel noise for sampling and enabling AirComp solutions via NOMA. Comparison with digital schemes follow along similar lines as in \cite{amiri2020machine, Xing2020Federated}, and is left as future work.

\subsection{Signal Design in OMA}
In OMA, in each communication block,  a single sample ${{\mathbf{m}} \vartheta}_k^{(s)} $ is transmitted by the scheduled device $k$ to the server. Dropping temporarily the superscript $s$ to simplify the notation, we assume uncoded transmission of the form
\begin{align}\label{eq: tx signal}
\text{(Transmitted Signal)}\quad {\mathbf{x}}_k&= {\mathbf{U}}_k  {{\mathbf{m}} \vartheta}_k,
\end{align}
where ${\mathbf{U}}_k$ is an $\mathbb{R}^{m_t\times d}$ precoding matrix. We assume that ${\mathbf{U}}_k$ factorizes into a zero-forcing pre-equalization matrix ${\mathbf{H}}_k^\dagger={\mathbf{H}}_k^{\sf T}({\mathbf{H}}_k {\mathbf{H}}_k^{\sf T} )^{-1}$ and a $m_r\times d$ encoding matrix as ${\mathbf{U}}_k={\mathbf{H}}_k^\dagger {\mathbf{E}}_k$. The use of  zero-forcing  is aligned with the majority of studies on uncoded transmission for machine learning \cite{zhu2018mimo,zhu2020over,wen2019reduced,liu2020privacy}, although generalizations are possible. The encoding matrix ${\mathbf{E}}_k$ is fixed, and not subject to optimization. For instance, if $m_r$ is a multiple of $d$, e.g., $m_r=ld$ for some integer $l$, one can set ${\mathbf{E}}_k\propto {\bf 1}_l \otimes {\mathbf{I}}_d$  to implement a form of analog repetition coding. Optimization of analog coding strategies  is left for future work. The encoding matrix ${\mathbf{E}}_k$ is subject to the power constraint \eqref{eq: power constraint}, which imposes the inequality
\begin{align}
\mathsf{E}[\|{\mathbf{x}}_k\|^2]=\mathsf{E}[\mathrm{tr}({\mathbf{U}}_k {{\mathbf{m}}  \vartheta}_k {{\mathbf{m}} \vartheta}_k^{\sf T}{\mathbf{U}}_k ^{\sf T})]\leq P. \label{eq: origin pc}
\end{align}
Since the distribution of the samples $\{{{\mathbf{m}} \vartheta}_k\}_{k=1}^K$ is unknown, we will specifically impose \eqref{eq: origin pc} by averaging over the available $S$ samples as
\begin{align}
\frac{1} {S}\sum_{s=1}^S\mathrm{tr}\left( \mathsf{E}\big[ \big({\mathbf{H}}_k {\mathbf{H}}_k^{\sf T})^{-1}\big]{\mathbf{E}}_k {{\mathbf{m}} \vartheta}_k^{(s)}  ({{\mathbf{m}} \vartheta}_k^{(s)})^{\sf T} {\mathbf{E}}_k^{\sf T}  \right) \leq P . \label{eq: approx pc}
\end{align}
The inequality \eqref{eq: approx pc} can be thought of an energy constraint across the transmission of the $S$ samples for each worker $k$.
For the example of analog repetition coding, constraint \eqref{eq: approx pc} is satisfied with ${\mathbf{E}}_k=\sqrt{{P_k}} {\bf 1}_l \otimes {\mathbf{I}}_d $,  where $P_k= PS/\sum_{s=1}^S\mathrm{tr}\big( \mathsf{E}\big[ \big({\mathbf{H}}_k {\mathbf{H}}_k^{\sf T})^{-1}\big]({\bf 1}_l \otimes {\mathbf{I}}_d ) {{\mathbf{m}} \vartheta}_k^{(s)}  ({{\mathbf{m}} \vartheta}_k^{(s)})^{\sf T} ({\bf 1}_l \otimes {\mathbf{I}}_d )^{\sf T}  \big) $.

 At the receiver side, for OMA, we propose to aggregate the received signals, ${\mathbf{y}}_1,\cdots, {\mathbf{y}}_K$ in a manner akin to  \eqref{eq: agg func} as
\begin{align}\label{eq: est theta}
\text{(Aggregation OMA)}\quad {{\mathbf{m}} \vartheta}=F({\mathbf{y}}_1,\cdots,{\mathbf{y}}_K)=\sum_{k=1}^K {\mathbf{W}}_k{\mathbf{y}}_k,
\end{align}
where ${\mathbf{W}}_1, \cdots, {\mathbf{W}} _K$ are  $d\times m_r$ weight matrices.  Note that the aggregation in \eqref{eq: est theta} is applied separately for each set of received signals $\{{{\mathbf{y}}}_k^{(s)}\}_{k=1}^K$. The weights $\{{\mathbf{W}}_k\}_{k=1}^K$ are to be optimized based on the received signals $\{{\mathbf{y}}_k^{(s)}\}_{k=1}^K$ for all samples $s=1,\dots,S$, as well as, possibly, based on the data set $\mathcal{Z}$.

\subsection{Signal Design in NOMA}
Under NOMA, in each communication block, all workers transmit one of their local samples simultaneously by using uncoded transmission as in \eqref{eq: tx signal}. As in OMA, matrix ${\mathbf{U}}_k$ is factorized as the product ${\mathbf{U}}_k={\mathbf{H}}_k^\dagger {\mathbf{E}}_k$. Furthermore, to enable the implementation of AirComp in a manner similar to \cite{zhu2018mimo,wen2019reduced,liu2020privacy}, we assume that all the transmitted signals are aligned by using the same encoding matrix, which denoted as ${\mathbf{E}}_k={\mathbf{E}}$ for all $k=1,\dots, K$. This matrix is subject to power constraint  \eqref{eq: power constraint} for all workers $k=1,\dots, K$ as
\begin{align} \label{eq: power constraint NOMA}
\frac{1} {S}\sum_{s=1}^S\mathrm{tr}\left( \mathsf{E}\big[ \big({\mathbf{H}}_k {\mathbf{H}}_k^{\sf T})^{-1}\big]{\mathbf{E}} {{\mathbf{m}} \vartheta}_k^{(s)}  ({{\mathbf{m}} \vartheta}_k^{(s)})^{\sf T} {\mathbf{E}}^{\sf T}  \right) \leq P .
\end{align}
 For the mentioned example of  repetition coding, the encoding matrix should be set as ${\mathbf{E}}=\sqrt{{\min_kP_k}} {\bf 1}_l \otimes {\mathbf{I}}_d $ in order to satisfy the constraint in \eqref{eq: power constraint NOMA}, where the value of $P_k $ is same as for~OMA.

The received signal \eqref{eq: rev signal noma} in NOMA carries out aggregation ``on air". Applying a single weight matrix ${\mathbf{W}}\in\mathbb{R}^{d\times m_r}$ yields the global sample
\begin{align}
\text{(Aggregation NOMA)} \quad {{\mathbf{m}} \vartheta}=F({\mathbf{y}}) = {\mathbf{W}} {\mathbf{y}}.  \label{eq: est theta NOMA}
\end{align}
The weight matrix ${\mathbf{W}}$ is again applied separately for each received signal ${\mathbf{y}}^{(s)}$, and is optimized by using all received signals $\{{\mathbf{y}}^{(s)}\}_{s=1}^S$, as well as, possibly, data set $\mathcal{Z}$. {\color{black}Note that, while in OMA one can potentially apply any aggregation function to the received samples $\{\{{\mathbf{y}}_k^{(s)}\}_{s=1}^S\}_{k=1}^K$, over-the-air computing via NOMA limits the aggregation to the linear superposition \eqref{eq: rev signal noma}
 carried out by the channel.}

To unify notation, we will write ${\mathbf{W}}=\{{\mathbf{W}}_k\}_{k=1}^K$ for the collection of weight matrices for OMA.

\section{Wireless Gaussian Consensus Monte Carlo}\label{sec: WGCMC}

In this section, we present a first wireless CMC protocol that is inspired by GCMC, and referred to as wireless Gaussian consensus Monte Carlo (WGCMC). Unlike GCMC, WGCMC can be proved to be asymptotically exact with respect to  the number of samples, $S$, under the described wireless channel model  when the  subposteriors are Gaussian. The limitations of WGCMC will motivate the introduction of a general optimization procedure for the decoding matrices in the next section.

To proceed, as in GCMC (see Sec. \ref{sec: GCMC}), we assume that each local subposterior $p_k({{\mathbf{m}} \vartheta}|\mathcal{Z}_k)$ is a  Gaussian $\mathcal{N}( {\bf 0} ,{\mathbf{C}}_k)$ distributions,  so that the global posterior is the Gaussian $\mathcal{N}({\bf 0} ,{\mathbf{C}})$ distribution, where ${\mathbf{C}}=(\sum_{k=1}^K{\mathbf{C}}_k^{-1})^{-1}$. Furthermore, in order to simplify the notation, we set $m_t=m_r=d$, and we choose the trivial encoding matrices ${\mathbf{E}}_k=\sqrt{{P_k}}{\mathbf{I}}_d $ for OMA and  ${\mathbf{E}}=\sqrt{\min_k{P_k}}{\mathbf{I}}_d $ for NOMA, where  $P_k= {PS}/\sum_{s=1}^S\mathrm{tr}\big( \mathsf{E}\big[ \big({\mathbf{H}}_k {\mathbf{H}}_k^{\sf T})^{-1}\big]{{\mathbf{m}} \vartheta}_k^{(s)}  ({{\mathbf{m}} \vartheta}_k^{(s)})^{\sf T}   \big) $. The generalization to any set of encoding matrices follows directly from the discussion below.

\subsection{Wireless GCMC: OMA}
Let us consider OMA first. For any weight matrices ${\mathbf{W}}$, the distribution of the output of the aggregated function \eqref{eq: est theta}  is given as
\begin{align}\label{eq: sample wgcmc oma}
{{\mathbf{m}} \vartheta}=\sum_{k=1}^K{ {\mathbf{W}}_k{\mathbf{y}}_k} \sim \mathcal{N}\left({\bf 0},\sum_{k=1}^{K} {\mathbf{W}}_k\big({P_k}{\mathbf{C}}_k+{N_0} {\mathbf{I}}\big) {\mathbf{W}}_k^{\sf T}  \right).
\end{align}
Accordingly,  if we know the covariance matrices $\{{\mathbf{C}}_k\}_{k=1}^K$, the power of encoded sample $P_k$, and noise variance $N_0$, optimal weight matrices can be obtained by solving the following equality
\begin{align}\label{eq: distr equality}
\sum_{k=1}^{K} {\mathbf{W}}_k\big({P_k}{\mathbf{C}}_k+{N_0} {\mathbf{I}}\big) {\mathbf{W}}_k^{\sf T}  =\Big(\sum_{k=1}^K{\mathbf{C}}_k^{-1}\Big)^{-1}.
\end{align}
In fact, this condition ensures that the distribution of the global sample ${{\mathbf{m}} \vartheta}=\sum_{k=1}^K {\mathbf{W}}_k {\mathbf{y}}_k $ is exactly the desired posterior, i.e., ${{\mathbf{m}} \vartheta}=\sum_{k=1}^K {\mathbf{W}}_k {\mathbf{y}}_k\sim\mathcal{N}({\bf 0} ,{\mathbf{C}})$.
A  solution to this system of equation is given by the matrices
\begin{align}
{\mathbf{W}}_k=\Big(\sum_{k=1}^K{\mathbf{C}}_k^{-1}\Big)^{-1} {\mathbf{C}}_k^{-1/2}\big({P_k}{\mathbf{C}}_k+{N_0} {\mathbf{I}}\big)^{-1/2}.  \label{eq: asymp opt}
\end{align}
Solution \eqref{eq: asymp opt} is not computable since the matrices ${\mathbf{C}}_1,\cdots, {\mathbf{C}}_K$ are unknown. They can be estimated by using the noisy received samples $\{{\mathbf{y}}_1^{(s)},\cdots, {\mathbf{y}}_K^{(s)}\}_{s=1}^{S}$ as
\begin{align}\label{est: ck}
\widehat{{\mathbf{C}}}_k= \! \frac{1}{P_k}\left[ \!\frac{1}{ (S-1)}\! \sum_{s=1}^{S} \left({\mathbf{y}}_k^{(s)}\!-\! \widehat{{\mathbf{m}} \mu}_{{\mathbf{y}}_k}\! \right)\! \left({\mathbf{y}}_k^{(s)} - \widehat{{\mathbf{m}} \mu}_{{\mathbf{y}}_k} \right)^{\sf T} \!\!-\!{N_0} {\mathbf{I}}\right]^+,
\end{align}
where $\widehat{{\mathbf{m}} \mu}_{{\mathbf{y}}_k}= \frac{1}{S}\sum_{s=1}^S {\mathbf{y}}_k^{(s)}$.  Note that the estimate $\widehat{{\mathbf{C}}}_k$ in \eqref{est: ck} is positive semidefinite thanks to the use of operation $[\ \cdot \ ]^+$. Plugging \eqref{est: ck} into \eqref{eq: asymp opt} yields WGCMC scheme
\begin{align}\label{eq: est WGCMC result}
&\text{(WGCMC for OMA)} \nonumber\\
& {{\mathbf{m}} \vartheta}^{(s)}= \sum_{k=1}^K  \underbrace{  \Big(  \sum_{k'=1}^K\widehat{{\mathbf{C}}}_{k'}^{-1}\Big)^{-1} \widehat{{\mathbf{C}}}_k^{-1/2} (P_k\widehat{{\mathbf{C}}}_k+N_0 {\mathbf{I}})^{-1/2}} _{{\mathbf{W}}_k} {{\mathbf{y}}}_k^{(s)}
\end{align}

The estimated covariance $\widehat{{\mathbf{C}}}_k $ in \eqref{est: ck} approaches  the true covariance ${\mathbf{C}}_k$ as  $S \rightarrow \infty$. Under this condition, the distribution of the global sample ${\mathbf{m}} \vartheta$ in \eqref{eq: sample wgcmc oma} equals the desired global posterior. In contrast,  the conventional GCMC scheme when applied directly to the received signals, i.e., by replacing scaled received signals $\{{\mathbf{y}}_k^{(s)}/\sqrt{P_k}\}_{s=1}^S$ for $\{{{\mathbf{m}} \vartheta}_k^{(s)}\}_{s=1}^S$ in \eqref{eq: est GCMC result}, is not asymptotically optimal as $S\rightarrow \infty$ since it disregards the presence of noise on the channel.

\subsection{Wireless Gaussian CMC: NOMA} \label{sec: WGCMC noma}
To develop an asymptotically optimal solution for NOMA, we further assume that local subposteriors are identical Gaussian $\mathcal{N}(0, {\mathbf{C}}_0)$ distributions, so that the global posterior is the Gaussian $\mathcal{N}({\bf 0} , \frac{1}{K}{\mathbf{C}}_0)$ distribution.   This highly simplified scenario with homogeneous posteriors is considered here to illustrate the potential benefits, but also the possible limitations, of NOMA.

The distribution of the output of the aggregated function \eqref{eq: est theta}  is given as
\begin{align}
{ {\mathbf{W}}{\mathbf{y}}} \sim \mathcal{N}\left({\bf 0},{\mathbf{W}}  \big(  \min_k P_k K {\mathbf{C}}_0+{N_0} {\mathbf{I}}\big) {\mathbf{W}}^{\sf T}  \right).
\end{align}
In a manner similar to OMA, the optimal weight matrix can be obtained by solving the following equality
\begin{align}\label{eq: distr equality GCMC NOMA}
{\mathbf{W}}  \big(  \min_k P_k K {\mathbf{C}}_0+{N_0} {\mathbf{I}}\big) {\mathbf{W}}^{\sf T}  = \frac{1}{K}{\mathbf{C}}_0,
\end{align}
and one of solutions is given as
\begin{align}
{\mathbf{W}}=\frac{1}{\sqrt{K}}{\mathbf{C}}_0^{1/2}\big(  K \min_k P_k  {\mathbf{C}}_0+{N_0} {\mathbf{I}}\big) ^{-1/2}. \label{eq: asymp opt NOMA}
\end{align}
The unknown matrix ${\mathbf{C}}_0$ can be estimated by using the noisy received samples $\{{\mathbf{y}}^{(s)}\}_{s=1}^{S}$ as
\begin{align}\label{est: ck noma}
\widehat{{\mathbf{C}}}_0\!= \!\frac{1}{K\min_kP_k}\!\left[ \! \frac{1}{ (S-1)} \! \sum_{s=1}^{S} \! \left( {\mathbf{y}}^{(s)} \!- \! \widehat{{\mathbf{m}} \mu}_{{\mathbf{y}}} \right) \!\! \left({\mathbf{y}}^{(s)}\! - \! \widehat{{\mathbf{m}} \mu}_{{\mathbf{y}}}\!\right)^{\sf T} \!\!-\!{N_0} {\mathbf{I}}\right]^+\!\!\!\!,
\end{align}
where $\widehat{{\mathbf{m}} \mu}_{{\mathbf{y}}}= \frac{1}{S}\sum_{s=1}^S {\mathbf{y}}^{(s)}$. Plugging \eqref{est: ck noma} into \eqref{eq: distr equality GCMC NOMA} yields the WGCMC scheme under NOMA,~i.e.,
\begin{align}\label{eq: est WGCMC result NOMA}
&\text{(WGCMC for NOMA)}\nonumber\\
&{{\mathbf{m}} \vartheta}^{(s)}= \underbrace{\frac{1}{\sqrt{K}}\widehat{{\mathbf{C}}}_0^{1/2}\big(  K \min_k P_k  \widehat{{\mathbf{C}}}_0+{N_0} {\mathbf{I}}\big) ^{-1/2}} _{{\mathbf{W}}} {\mathbf{y}}^{(s)}.
\end{align}

As for OMA, the sample ${{\mathbf{m}} \vartheta}={\mathbf{W}}{\mathbf{y}}$ generated by of WGCMC under NOMA in \eqref{eq: asymp opt NOMA} follows the desired global posterior distribution when  $S \rightarrow \infty$. The potential advantage of NOMA versus OMA is that, with the same number of communication blocks, $T$, the number of received samples is $S=T$ in NOMA while it equals $S=T/K$ in OMA. The large number of samples is generally advantageous for the downstream applications that use the samples, e.g., for ensemble prediction. Furthermore, a larger $S$ leads to a more accurate estimation of the relevant covariance matrices. On the flip side, optimality was only argued here for a simplified homogeneous setting, and one may expect a degraded performance in heterogeneous settings. A strategy that addresses the case of general subposterior is introduced next.

\section{Wireless Variational Consensus Monte Carlo} \label{sec: WVCMC}
The WGCMC schemes proposed in the previous subsection can be applied to any class of subposteriors in a manner similar to GCMC. However, apart form the special case of Gaussian subposteriors, this application is expected to be successful only in the regime of large data sets due to the Bernstein-von Mises Theorem.  In general, the WGCMC scheme is suboptimal. In this section, we introduce an optimization strategy for the design of the weight matrices ${\mathbf{W}}$ for OMA and NOMA that is able to adapt to arbitrary subposteriors and channels.

\subsection{Design Principle}
Inspired by \cite{rabinovich2015variational},
we adopt a variational approach to formulate the problem of ensuring that the distribution $q_F({{\mathbf{m}} \vartheta})$ of the aggregated samples \eqref{eq: est theta} for OMA and \eqref{eq: est theta NOMA} for NOMA is close to the true posterior $p({{\mathbf{m}} \vartheta}|\mathcal{Z})$.  We recall that reference \cite{rabinovich2015variational} assumed noiseless channels and hence focused on optimizing the aggregation function \eqref{eq: agg func} instead.
This condition can be expressed as the approximate equality
\begin{subequations} \label{eq: approx global posterior}
\begin{align}
 &\text{(Mean Field Approximation for OMA)} \nonumber\\
 &  p({{\mathbf{m}} \vartheta}|\mathcal{Z}) \approx  q_F({{\mathbf{m}} \vartheta})\nonumber\\
 &= \!\!\int\!\!\! \bigg( \! \prod_{k=1}^K  p({\mathbf{y}}_k)\!\! \bigg) \delta\left[{{\mathbf{m}} \vartheta}  = F({\mathbf{y}}_1,\cdots,{\mathbf{y}}_K)\right] d {\mathbf{y}}_1, \cdots d {\mathbf{y}}_K  \label{eq: mf oma y}\\
 &= \!\!\int\!\!\! \bigg( \prod_{k=1}^K p_k({{\mathbf{m}} \vartheta}_k |\mathcal{Z}_k) p({\mathbf{n}}_k) \bigg)\delta\left[{{\mathbf{m}} \vartheta}  = F({\mathbf{y}}_1,\cdots,{\mathbf{y}}_K)\right]    d {{\mathbf{m}} \vartheta}_1, d{\mathbf{n}}_1, \nonumber\\
 & \qquad \qquad  \qquad \qquad \qquad \qquad \quad \qquad \cdots  ,d{{\mathbf{m}} \vartheta}_K, d {\mathbf{n}}_K  \!  \label{eq: mf oma theta} \\
 &= \mathsf{E}_{{{\mathbf{m}} \vartheta}_1,\cdots,{{\mathbf{m}} \vartheta}_K, {\mathbf{n}}_1,\cdots, {\mathbf{n}}_K}\Big[\delta\left[{{\mathbf{m}} \vartheta}  = F({\mathbf{y}}_1,\cdots,{\mathbf{y}}_K)\right] \Big],
  \end{align}
  \end{subequations}
  for OMA, where $p({\mathbf{n}}_k)$ is the $\mathcal{N}(0,N_0 {\mathbf{I}}_{m_r})$ probability density function (pdf) and  $\delta[\cdot]$ is the Dirac delta function; and
\begin{subequations}\label{eq: approx global posterior NOMA}
\begin{align}
 &\text{(Mean Field Approximation for NOMA)}    \nonumber\\
 & p({{\mathbf{m}} \vartheta}|\mathcal{Z}) \approx   q_F({{\mathbf{m}} \vartheta})\nonumber\\
  &=\int  p({\mathbf{y}}) \delta\left[{{\mathbf{m}} \vartheta}  = F({\mathbf{y}})\right] d {\mathbf{y}} \label{eq: mf noma y} \\
&= \int  \bigg(p({\mathbf{n}}) \prod_{k=1}^K p_k({{\mathbf{m}} \vartheta}_k |\mathcal{Z}_k) \bigg) \delta\left[{{\mathbf{m}} \vartheta}  = F({\mathbf{y}})\right] d {\mathbf{n}},  d {{\mathbf{m}} \vartheta}_1,\cdots  ,d{{\mathbf{m}} \vartheta}_K   \qquad \quad \qquad \qquad \label{eq: mf noma theta} \\
&=\mathsf{E}_{{{\mathbf{m}} \vartheta}_1,\cdots,{{\mathbf{m}} \vartheta}_K, {\mathbf{n}}} \Big[\delta\left[{{\mathbf{m}} \vartheta}  = F({\mathbf{y}})\right] \Big],
\end{align}
\end{subequations}
for NOMA, where $p({\mathbf{n}})$ is the $\mathcal{N}(0,N_0 {\mathbf{I}}_{m_r})$ pdf. The distribution $q_F({{\mathbf{m}} \vartheta})$ in \eqref{eq: mf oma theta} and \eqref{eq: mf noma theta} of the obtained global sample is computed by averaging over the distribution of the local samples $({{\mathbf{m}} \vartheta}_1,\cdots,{{\mathbf{m}} \vartheta}_K)$ and of the channel noise. Note that, for OMA, as per \eqref{eq: mf oma y}, the received signal $\{{\mathbf{y}}_k\}_{k=1} ^K$ are independent due to the independence of local samples and noise samples.

In order to enforce the approximate equalities \eqref{eq: approx global posterior} through the optimization of the weight matrices ${\mathbf{W}}$, we propose to tackle the problem of minimizing the distance of the distributions $q_F({{\mathbf{m}} \vartheta})$ and $  p({{\mathbf{m}} \vartheta}|\mathcal{Z})$ as measured by Kullback-Leibler (KL) divergence
$${\rm KL}\big(q_F({{\mathbf{m}} \vartheta}) \| p({{\mathbf{m}} \vartheta}|\mathcal{Z})\big)=\mathsf{E}_{q_F({{\mathbf{m}} \vartheta})}\left[\log \left(\frac{q_F({{\mathbf{m}} \vartheta}) }{ p({{\mathbf{m}} \vartheta}|\mathcal{Z})}\right)\right].$$
This is the standard objective function used in variational inference (see, e.g., \cite{simeone2017brief, rabinovich2015variational}). This problem is equivalent to the minimization of the variational free energy \cite{josefree}
   \begin{align}\label{eq: max ELBO}
 \min_{{\mathbf{W}}}\Big\{ \mathcal{L}({\mathbf{W}})& = -\mathsf{E}_{{{\mathbf{m}} \vartheta}\sim q_{F}({\mathbf{m}} \vartheta)}\big[\log {p({{\mathbf{m}} \vartheta},\mathcal{Z})}\big] -{\sf H}\big[ q_F({\mathbf{m}} \vartheta)\big] \Big\},
 \end{align}
which is also known as  negative evidence lower bound (ELBO) \cite{josefree}. Note that the variational free energy  \eqref{eq: max ELBO} does not depend on the unknown posterior $p({{\mathbf{m}} \vartheta}|\mathcal{Z})$ but on the joint $p({{\mathbf{m}} \vartheta},\mathcal{Z})=p({{\mathbf{m}} \vartheta})p(\mathcal{Z}|{{\mathbf{m}} \vartheta})=p({{\mathbf{m}} \vartheta})\prod_{k=1}^Kp(\mathcal{Z}_k|{{\mathbf{m}} \vartheta})$.

In the following, we will develop algorithms that tackle problem \eqref{eq: max ELBO}  for OMA by using the received signals $\{\{{\mathbf{y}}_k^{(s)}\}_{k=1}^K\}_{s=1}^S$ and for NOMA by using the received signals $\{{\mathbf{y}}^{(s)}\}_{s=1}^S$. In both cases, the server can also use the data set $\mathcal{Z}$. The proposed approaches are based on stochastic gradient descent (SGD) applied at the server. Specifically, the server applies SGD steps based on the gradient $\nabla_{{\mathbf{m}} \vartheta} \log p({\mathbf{z}}_n|{{\mathbf{m}} \vartheta})$ for source data point ${\mathbf{z}}_n\in \mathcal{Z}$, which can be readily obtained at the server based on the global data set $\mathcal{Z}$. In fact, computing gradients $\nabla_{{\mathbf{m}} \vartheta} \log p({\mathbf{z}}_n|{{\mathbf{m}} \vartheta})$ is much less computationally intensive than computing the posterior $p({{\mathbf{m}} \vartheta}|\mathcal{Z})$: The latter requires integrating over the space of the model parameters, while the former can be carried out via efficient algorithms such as backpropagation.  For comparison, in Sec \ref{sec: sim}, we will consider the performance of  stochastic gradient Langevin dynamics (SGLD), an MCMC scheme that also requires the computation of the same stochastic gradients.

\subsection{Wireless Variational CMC: OMA}

\subsubsection{Approximation of the Objective Function}
For OMA, the global sample ${{\mathbf{m}} \vartheta}\sim q_F({{\mathbf{m}} \vartheta})$ is given as a function of the received signals in \eqref{eq: est theta}. Therefore, the objective in \eqref{eq: max ELBO} can be rewritten as the expectation over the received signals  $\{{\mathbf{y}}_k\}_{k=1}^K$ as
   \begin{align} \label{eq: min KL repar}
  \mathcal{L}({\mathbf{W}})&= -\mathsf{E}_{{\mathbf{y}}_1,\cdots,{\mathbf{y}}_K}\bigg[\log {p\Big( \sum_{k=1}^K {\mathbf{W}}_k {\mathbf{y}}_k,\mathcal{Z}\Big)}\bigg] -{\sf H}\big[ q_F({\mathbf{m}} \vartheta)\big].
 \end{align}
 This step implements the  ``reparameterization trick" \cite{simeone2017brief}, moving the optimization variables ${\mathbf{W}}$ from the averaging distribution to the argument  of the expectation.
The expectation can now be estimated by averaging over the samples $\{\{{\mathbf{y}}_k^{(s)}\}_{k=1}^K\}_{s=1}^S$ as
\begin{align}
 \widehat{ \mathcal{L}}({\mathbf{W}})& \!= \!-\frac{1}{S}\sum_{s=1}^{S}\Big[\log {p\Big( \sum_{k=1}^K {\mathbf{W}}_k {\mathbf{y}}_k^{(s)},\mathcal{Z}\Big)}\Big]  -  {\sf H}\big[ q_F({\mathbf{m}} \vartheta)\big].  \!\label{eq: obj sum}
\end{align}
The  direct minimization of \eqref{eq: obj sum}  is still intractable due to the presence of differential entropy of the distribution $q_F({{\mathbf{m}} \vartheta})$. This issue is addressed via a tractable lower bound, as derived in the following lemma.
 \begin{lemma}[Lower Bound on Differential Entropy for OMA]\label{lemma: LB entropy}\emph{With given encoding matrices $\{{\mathbf{E}}_k\}_{k=1}^K$,  the differential  entropy ${\sf H}[ q_F({\mathbf{m}} \vartheta)]$ of the distribution $q_F({{\mathbf{m}} \vartheta})$ of the global sample \eqref{eq: est theta} under  OMA is lower bounded~as
 \begin{align}\label{eq: LB entropy}
{\sf H}[ q_F({\mathbf{m}} \vartheta)] &\geq \frac{d}{2}\log\left(2K\sqrt{2\pi e N_0}\right)+   \frac{1}{2K}\sum_{k=1}^K  \log |\det ({\mathbf{W}}_k{\mathbf{E}}_k)| \nonumber\\
& \qquad \qquad  \qquad+{\sf H}[p_k({{\mathbf{m}} \vartheta}_k|\mathcal{Z}_k) ] +\frac{1}{2}\log \det \left({\mathbf{W}}_k {\mathbf{W}}_k^{\sf T }\right)  \nonumber\\
&:=\widetilde{\sf H}[ q_F({\mathbf{m}} \vartheta)],
\end{align}
where ${\sf H}[p_k({{\mathbf{m}} \vartheta}_k|\mathcal{Z}_k) ]$ is the differential entropy of the subposterior  $p_k({{\mathbf{m}} \vartheta}_k|\mathcal{Z}_k)$.}
\end{lemma}
\noindent{\emph{Proof:} } The proof follows the steps of \cite[Supplement A Proofs]{rabinovich2015variational} and is detailed in Appendix~\ref{proof entropy}.

 Using Lemma~\ref{lemma: LB entropy} in  \eqref{eq: obj sum} yields the following upper bound on the empirical free energy  \eqref{eq: obj sum} for~OMA
 \begin{align}
\widehat{\mathcal{L}}({\mathbf{W}})\leq  \widetilde{\mathcal{L}}({\mathbf{W}}) := -\frac{1}{S}\sum_{s=1}^{S}\Big[\log {p\Big( \sum_{k=1}^K  {{\mathbf{W}}_k {\mathbf{y}}_k^{(s)}},\mathcal{Z}\Big)}\Big] \nonumber\\
 -\widetilde{\sf H}\big[ q_F({\mathbf{m}} \vartheta)\big]. \label{obj OMA}
\end{align}

\subsubsection{Stochastic Gradient Descent Algorithm}
 The proposed strategy, referred to as wireless variational consensus Monte Carlo (WVCMC), aims at minimizing the upper bound \eqref{obj OMA} via SGD. Given a mini-batch $\mathcal{B} \in \mathcal{Z}$ of $N_b\leq N$ samples, stochastic gradient of $\widetilde{L}({\mathbf{W}})$ with respect to weight matrix ${\mathbf{W}}_k$ is given  as
\begin{align}
&\!\!\!\widehat{\nabla}_{{\mathbf{W}}_k} \widetilde{\mathcal{L}}({\mathbf{W}}) \nonumber\\
&\!\!\!=\!-\frac{1}{S}\!\!\sum_{s=1}^{S} \!\!\Big[\!\nabla_{{{\mathbf{m}} \vartheta}^{(s)}}\! \log p({{\mathbf{m}} \vartheta}^{(s)})\!+\! \frac{N}{N_b}\!\!\sum_{{\mathbf{z}}_n\in \mathcal{B}}\!\! \nabla_{{{\mathbf{m}} \vartheta}^{(s)}} \!\!\log p({\mathbf{z}}_n|{{\mathbf{m}} \vartheta}^{(s)}\!)\!  \Big] \!({\mathbf{y}}_k^{(s)})\!^{\sf T} \nonumber\\
& \qquad \qquad \qquad \qquad -\frac{1}{2K} \big[ ({\mathbf{W}}_k{\mathbf{E}}_k)^{-\sf T}  {\mathbf{E}}_k^{\sf  T}+({\mathbf{W}}_k^\dagger)^{\sf T}  \big] ,\label{eq: gradient Dk}
 \end{align}
where we recall that we have the global sample ${{\mathbf{m}} \vartheta}^{(s)}=\sum_{k=1}^K {\mathbf{W}}_k {\mathbf{y}}_k^{(s)}$; and ${\mathbf{W}}_k^\dagger $ is the pseudoinverse of ${\mathbf{W}}_k$. This yields SGD algorithm based on \eqref{eq: gradient Dk} as summarized in  Algorithm~\ref{algorithm: offline optimization for Dk}.  Note that, unlike GCMC and WGCMC, WVCMC requires the use of data set $\mathcal{Z}$ at server.
 \begin{algorithm}
{\bf Input:} Data set $\mathcal{Z}$, received signals   $\{{\mathbf{y}}_1^{(s)},\cdots, {\mathbf{y}}_K^{(s)}\}_{s=1}^{S}$, encoding matrices$\{{\mathbf{E}}_k\}_{k=1}^K$, step size $\eta$, number of iterations $t_m$.

{\bf Initialize:}   ${{\mathbf{W}}}^0$

{\bf For each iteration:} $t=1,\dots,t_m$

\quad  {\bf Estimate samples of global posterior using ${{\mathbf{W}}}^{t-1}$:}

\qquad       ${{\mathbf{m}} \vartheta}^{(s)}=\sum_{k=1}^K{ {\mathbf{W}}_k^{t-1}{\mathbf{y}}_k^{(s)}} $ for $s=1, \dots, S$

\quad  {\bf Update ${{\mathbf{W}}}^{t}$:}

\qquad    Select a mini-batch of data $\mathcal{B}=\{{\mathbf{z}}_n\}_{n=1}^{N_b} \in \mathcal{Z}$

\qquad     $ {{\mathbf{W}}}_k^{t}={{\mathbf{W}}}_k^{t-1}-\eta \widehat{\nabla}_{{\mathbf{W}}_k} \widetilde{\mathcal{L}}({\mathbf{W}}^{t-1})  $

\qquad   \quad       with $\widehat{\nabla}_{{\mathbf{W}}_k} \widetilde{\mathcal{L}}({\mathbf{W}}^{t-1})$ in \eqref{eq: gradient Dk}

{\bf Output:} Global samples ${{\mathbf{m}} \vartheta}^{(s)}=\sum_{k=1}^K{ {\mathbf{W}}_k^{t_m}{\mathbf{y}}_k^{(s)}}$  for $s=1,\dots,S$.

\caption{WVCMC for OMA}
\label{algorithm: offline optimization for Dk}
\end{algorithm}

{\color{black}In \cite{rabinovich2015variational},  matrices $\{{\mathbf{W}}\}_{k=1}^K$ were assumed to be positive semidefinite, resulting in a convex program when the  joint distributions $p({{\mathbf{m}} \vartheta}, \mathcal{Z})$  are log-concave in ${\mathbf{m}} \vartheta$. In  the more general setting considered here, the problem is not convex.  The iteration complexity of Algorithm 1 to guarantee convergence to stationary points following directly from standard results on SGD, and it depends on the curvature properties of function $\log p({{\mathbf{m}} \vartheta}, \mathcal{Z})$ \cite[Theorem 4.8]{bottou2018optimization}.}

\subsection{Wireless Variational CMC: NOMA}
\subsubsection{Approximation of the Objective Function}
In a manner similar to OMA, we first approximate the objective in \eqref{eq: max ELBO}  by averaging over the received signal $\{{\mathbf{y}}^{(s)}\}_{s=1}^S$. To this end, we recall that the global samples ${{\mathbf{m}} \vartheta}^{(s)}$ are obtained via \eqref{eq: est theta NOMA}, and hence we can write the empirical free energy \eqref{eq: max ELBO} as
\begin{align}
\widehat{\mathcal{L}}({\mathbf{W}})&\!=\! -\frac{1}{S}\sum_{s=1}^{S}\Big[\log {p\Big({\mathbf{W}}{\mathbf{y}}^{(s)},\mathcal{Z}\Big)}\Big] -{\sf H}\big[ q_F({\mathbf{m}} \vartheta)\big].   \label{eq: obj sum noma}
 \end{align}
 Again, the direct optimization of \eqref{eq: obj sum noma}  is  intractable given the presence of differential entropy ${\sf H}[q_F({\mathbf{m}} \vartheta)]$, which is addressed by the tractable lower bound derived below.
 \begin{lemma}[Lower Bound on Differential Entropy for NOMA]\label{lemma: LB entropy noma}\emph{With given encoding matrices $\{{\mathbf{E}}_k\}_{k=1}^K$,  the differential  entropy ${\sf H}[ q_F({\mathbf{m}} \vartheta)]$ of the distribution $q_F({{\mathbf{m}} \vartheta})$ of the global sample \eqref{eq: est theta NOMA} under  NOMA is lower bounded~as
 \begin{align}\label{eq: LB entropy NOMA}
{\sf H}[ q_F({\mathbf{m}} \vartheta)] &\geq \frac{d}{2}\log\left[(K+1)\left({2 \pi e N_0}\right)^{\frac{1}{K+1}}\right]  \nonumber\\
&\quad +\frac{1}{K+1}\bigg(K\log |\det ({\mathbf{W}}{\mathbf{E}})|+ \frac{1}{2} \log \det({\mathbf{W}} {\mathbf{W}}^{\sf T}) \nonumber\\
&\qquad \qquad\qquad\qquad \qquad \qquad+ \sum_{k=1}^K   {\sf H}[p_k({{\mathbf{m}} \vartheta}_k|\mathcal{Z}_k) ] \bigg)   \nonumber\\
&:=\widetilde{\sf H}[ q_F({\mathbf{m}} \vartheta)].
\end{align}}
\end{lemma}

 Using Lemma~\ref{lemma: LB entropy noma} in  \eqref{eq: obj sum noma} yields the following upper bound on the empirical free energy for NOMA as
 \begin{align}
\widehat{\mathcal{L}}({\mathbf{W}})&\leq  \widetilde{\mathcal{L}}({\mathbf{W}}) := -\frac{1}{S}\sum_{s=1}^{S}\Big[\log {p\Big(  {{\mathbf{W}} {\mathbf{y}}^{(s)}},\mathcal{Z}\Big)}\Big]  -\widetilde{\sf H}\big[ q_F({\mathbf{m}} \vartheta)\big]. \label{obj NOMA}
\end{align}

\subsubsection{Stochastic Gradient Descent Algorithm}
The proposed WVCMC algorithm for NOMA aims at   minimizing the upper bound  \eqref{obj NOMA}  via SGD.
Given a mini-batch $\mathcal{B} \in \mathcal{Z}$ of $N_b\leq N$ samples, stochastic gradient of $\widetilde{L}({\mathbf{W}})$ with respect to weight matrix ${\mathbf{W}}$ is given  as
\begin{align}
&\!\!\widehat{\nabla}_{{\mathbf{W}}} \widetilde{\mathcal{L}}({\mathbf{W}}) \nonumber\\
&\!\!=\!\!-\frac{1}{S}\!\!\sum_{s=1}^{S}\!\!  \Big[\!\nabla_{{{\mathbf{m}} \vartheta}^{(s)}}\! \log p({{\mathbf{m}} \vartheta}^{(s)})\!+\!\frac{N}{N_b}\!\!\sum_{{\mathbf{z}}_n\in \mathcal{B}}\!\! \nabla_{{{\mathbf{m}} \vartheta}^{(s)}} \!\log \! p({\mathbf{z}}_n|{{\mathbf{m}} \vartheta}^{(s)}) \! \Big]\!  ({\mathbf{y}}^{(s)})\!^{\sf T}  \nonumber \\
&\qquad \qquad \qquad   -\frac{1}{K+1}\Big[ K ({\mathbf{W}}{\mathbf{E}})^{-\sf T}  {\mathbf{E}}^{\sf  T}+({\mathbf{W}}^\dagger)^{\sf T}   \Big],\label{eq: gradient Dk noma}
 \end{align}
where we have  ${{\mathbf{m}} \vartheta}^{(s)}= {\mathbf{W}}{\mathbf{y}}^{(s)}$. The resulting algorithm is summarized in  Algorithm~\ref{algorithm: offline optimization for Dk NOMA}.
 \begin{algorithm}
{\bf Input:} Data set $\mathcal{Z}$, received signals   $\{{\mathbf{y}}^{(s)}\}_{s=1}^{S}$, encoding matrices ${\mathbf{E}}$, step size $\eta$, number of iterations $t_m$.

{\bf Initialize:}   ${{\mathbf{W}}}^0$

{\bf For each iteration:} $t=1,\cdots,t_m$

\quad  {\bf Estimate samples of global posterior using  ${{\mathbf{W}}}^{t-1}$:}

\qquad       ${{\mathbf{m}} \vartheta}^{(s)}={ {\mathbf{W}}^{t-1}{\mathbf{y}}^{(s)}} $ for $s=1, \dots, S$

\quad  {\bf Update ${{\mathbf{W}}}^{t}$:}

\qquad    Select a mini-batch of data $\mathcal{B}=\{{\mathbf{z}}_n\}_{n=1}^{N_b} \in \mathcal{Z}$

\qquad     $ {{\mathbf{W}}}^{t}={{\mathbf{W}}}^{t-1}-\eta \widehat{\nabla}_{{\mathbf{W}}} \widetilde{\mathcal{L}}({\mathbf{W}}^{t-1})  $

\qquad \quad with $\widehat{\nabla}_{{\mathbf{W}}} \widetilde{\mathcal{L}}({\mathbf{W}}^{t-1})$ in \eqref{eq: gradient Dk noma}

{\bf Output:} Global samples ${{\mathbf{m}} \vartheta}^{(s)}={ {\mathbf{W}}^{t_m}{\mathbf{y}}^{(s)}}$, for $s=1,\dots,S$.

\caption{WVCMC for NOMA}
\label{algorithm: offline optimization for Dk NOMA}
\end{algorithm}

\section{Experiments}\label{sec: sim}

In this section, we evaluate the performance of the proposed channel-driven CMC schemes for two standard problems: 1) CMC with Gaussian subposteriors \cite{scott2016bayes}; and 2) CMC for Bayesian probit regression on  synthetic data set and  MNIST data set. To this end, we first introduce the performance metrics and reference schemes that are adopted throughout this section.

\subsubsection{Performance Metrics}
As in prior works on CMC \cite{rabinovich2015variational}, the performance of a CMC algorithm is evaluated by measuring the relative difference between two sample averages of a test function $f({{\mathbf{m}} \vartheta})$, namely the reference average obtained under the true global posterior \eqref{eq: global posterior} and the average computed using samples produced by the CMC algorithm.  Given a test function $f(\cdot)$, denoting ${{\mathbf{m}} \vartheta}_g^{(s)}\sim p({{\mathbf{m}} \vartheta}|\mathcal{Z})$ for $s=1,\dots,S_g$ as a set of samples from  the distribution of global posterior \eqref{eq: global posterior},  the test error is defined as
\begin{align}
{\sf err}(f(\cdot))=\frac{\big|\frac{1}{S}\sum_{s=1}^{S} f( {{\mathbf{m}} \vartheta}^{(s)}) -  \frac{1}{S_g}\sum_{s=1}^{S_g} f( {{\mathbf{m}} \vartheta}_g^{(s)}) \big|}{\big| \frac{1}{S_g}\sum_{s=1}^{S_g} f( {{\mathbf{m}} \vartheta}_g^{(s)}) \big|},  \label{eq: eval func}
\end{align}
where $\{{{\mathbf{m}} \vartheta}^{(s)}\}_{s=1}^S$ are the samples produced by the CMC algorithm under evaluation.
Specifically, as in  \cite{rabinovich2015variational}, we consider multiple test functions, with each function given by one entry of the outer product matrix ${{\mathbf{m}} \vartheta}{{\mathbf{m}} \vartheta}^{\sf T}$. Accordingly, we have the test functions $f_{ij}({\mathbf{m}} \vartheta)= {{\mathbf{m}} \vartheta}[i]  {{\mathbf{m}} \vartheta}[j]$ for $i,j=1,\dots, d$, where $ {{\mathbf{m}} \vartheta}[i]$ is the $i$-th entry in vector ${{\mathbf{m}} \vartheta}$. We then average the per-test function error \eqref{eq: eval func} to obtain
\begin{align}\label{eq: err2}
{\sf err}=\frac{1}{d^2}\sum_{i=1}^d \sum_{j=1}^d  \frac{\big|\frac{1}{S}\sum_{s=1}^{S} f_{ij}({{\mathbf{m}} \vartheta}^{^{(s)}}) -  \frac{1}{S_g}\sum_{s=1}^{S_g}f_{ij}({{\mathbf{m}} \vartheta}_g^{^{(s)}})  \big|}{\big|\frac{1}{S_g}\sum_{s=1}^{S_g} f_{ij}({{\mathbf{m}} \vartheta}_g^{^{(s)}})\big|}.
\end{align}
Accordingly,  the error \eqref{eq: err2} quantifies the accuracy of the empirical correlation matrix of the posterior distribution that can be obtained from the samples $\{{{\mathbf{m}} \vartheta}^{(s)}\}_{s=1}^S$. Similar results are obtained for the estimate of other moments of the posterior such as the mean vector.

{\color{black}
We will also consider the inference performance of the ensemble predictive distribution
\begin{align}
\widehat{p}_{\{{{\mathbf{m}} \vartheta}^{(s)}\}_{s=1}^S}(v  |{\mathbf{u}}, \mathcal{Z})=\frac{1}{S}\sum_{s=1}^S p(v|{{\mathbf{m}} \vartheta}^{(s)},{\mathbf{u}},\mathcal{Z}) \label{eq: pred dist1}
\end{align}
obtained at the server given the samples $\{{{\mathbf{m}} \vartheta}^{(s)}\}_{s=1}^S$ for any test input ${\mathbf{u}}$. Specifically, the quality of predictive distribution \eqref{eq: pred dist1}  is measured by the Kullback-Leibler (KL) divergence between \eqref{eq: pred dist1}  and the ideal ensemble distribution $\widehat{p}_{\{{{\mathbf{m}} \vartheta}_g^{(s)}\}_{s=1}^S}(v  |{\mathbf{u}}, \mathcal{Z}) $ computed by using the samples of true posterior distribution, i.e.,
\begin{align}\label{eq: KL divergence}
\frac{1}{N_t}\sum_{i=1}^{N_t}{\rm KL}\Big(\widehat{p}_{\{{{\mathbf{m}} \vartheta}^{(s)}\}_{s=1}^S}(v  |{\mathbf{u}}_i, \mathcal{Z}) \Big \| \widehat{p}_{\{{{\mathbf{m}} \vartheta}_g^{(s)}\}_{s=1}^S}(v  |{\mathbf{u}}_i, \mathcal{Z}) \Big),
\end{align}
where $\{{\mathbf{u}}_i\}_{i=1}^{N_t}$ is a given test set.
}

\subsubsection{Reference Schemes}
We evaluate the performance of:  (\emph i) GCMC in  \eqref{eq: est GCMC result} with decoded received signals $\{{\mathbf{E}}^\dagger  {\mathbf{y}}_k^{(s)} \}_{s=1}^S$  in lieu of $\{{{\mathbf{m}} \vartheta}_k^{(s)}\}_{s=1}^S$; (\emph{ii}) WGCMC in  \eqref {eq: est WGCMC result} for OMA and in \eqref{eq: est WGCMC result NOMA} for NOMA; and (\emph{iii}) WVCMC in Algorithm 1 and Algorithm~2.  Furthermore, in order to evaluate the advantages of leveraging the $K$ workers in the wireless data center, we also implement a state-of-the-art scheme for MC sampling at the server, namely SGLD \cite{welling2011bayesian}.  SGLD is selected because it can be run with the same order of computational complexity at the server of the proposed WVCMC schemes. In fact, as detailed next, SGLD relies on the computation of stochastic gradients $\nabla \log p({{\mathbf{m}} \vartheta}, {\mathbf{z}}_n)$ for randomly selected examples ${\mathbf{z}}_n\in \mathcal{Z}$. Therefore, we can compare SGLD and WVCMC under the assumption of (approximately) equal computation load at the server.

At each iteration $t$, SGLD updates the model parameter vector by using a subset of data in the global data set $\mathcal{Z}$ as
 \begin{align}
\text{(SGLD)} \quad  {{\mathbf{m}} \vartheta}^{(t+1)}\!=\!{{\mathbf{m}} \vartheta}^{(t)}\!\!+\!\!\frac{\eta^{(t)}}{2} &\Big[ \!\frac{N}{N_b}\sum_{{\mathbf{z}}_n\in \mathcal{B}} \nabla_{{{\mathbf{m}} \vartheta}^{(t)}} \log p({\mathbf{z}}_n|{{\mathbf{m}} \vartheta}^{(t)})  \nonumber\\
&\!\!+\!\nabla_{{{\mathbf{m}} \vartheta}^{(t)}}\! \log p({{\mathbf{m}} \vartheta}^{(t)})\!\Big]\!\!+\!{\mathbf{n}}^{(t)}\!,\! \label{eq: update SGLD}
 \end{align}
 where the learning rate is set as $\eta^{(t)} =\alpha(\beta+t)^{-\gamma}$ for some $\beta>0$ and $\gamma \in (0.5,1]$, and ${\mathbf{n}}^{(t)} \sim\mathcal{N}(0, \eta^{(t)} {\mathbf{I}})$ is injected noise. SGLD is known to produce, in the limit of a large number of iterations, samples from the desired posterior. After removing the first $t_b$ samples, which we consider to constitute the burn-in phase, the remaining $\{{{\mathbf{m}} \vartheta}^{(t)}\}_{t=t_b+1}^{t_m}$ samples are used for evaluation in \eqref{eq: err2}.  We emphasize that SGLD does not make use of the presence of the $K$ workers. Distributed versions of SGLD also exist. These can be useful as alternatives to the CMC schemes studied in this paper when the workers do not have sufficient computing power to obtain samples from the local subposteriors -- a condition that we have assumed. We refer to Sec. \ref{sec: conclusions} for more discussion on this point. As for the  computation at the server, consider $C$ computed gradients $\nabla_{{{\mathbf{m}} \vartheta}} \log p({{\mathbf{m}} \vartheta},{\mathbf{z}}_n)$, SGLD requires $C/N_b$ iterations while  WVCMC requires $C/N_bS$ iterations since WVCMC process $S$ samples $\{{{\mathbf{m}} \vartheta}^{(s)}\}_{s=1}^S$ in parallel.

\subsection{Gaussian Subposteriors}

We start with the baseline example studied in Sec. \ref{sec: GCMC} and Sec. \ref{sec: WGCMC}, in which the local subposterior for each worker $k$ follows a Gaussian distribution $\mathcal{N}(0,{\mathbf{C}}_k)$. The model dimension is $d=5$, and each covariance matrix ${\mathbf{C}}_k$ is a symmetric Toeplitz matrix with first column given as  $[ 1, \rho_k,    \rho_k^2,   \rho_k^3,  \rho_k^4]$ for some $\rho_k \in [-1,1] $. We consider $K = 10$ workers, and set $\rho_k=(k-1)/K$ for $k$-th worker.  For reference, we also study a homogeneous setting in which the covariances are equal to  ${\mathbf{C}}_0=K( \sum_{k=1}^K {\mathbf{C}}_k^{-1})^{-1}$ for all workers.

To start, we consider frequency-flat channels with ${\mathbf{H}}_k={\mathbf{I}}_d$, and the $\mathsf{SNR}=P/(dN_0)$ is set as $5$ dB by default. Unless stated otherwise, the number of communication blocks is equal to $T=2000$; the number of per-output-sample iterations  is $t_m=30$ for NOMA and $t_m=300$ for OMA, which ensures that the same total number of iterations, namely $300$, are carried out for both OMA and NOMA;   the learning rate is $\eta=10^{-3}$ for NOMA and $\eta=5 \times 10^{-3}$ for OMA, which have been separately optimized through a grid search;  and the initial point used for WVCMC is ${\mathbf{W}}^0=1/K{\mathbf{I}}_d$ for NOMA, while in OMA it is set as GCMC solution. This initialization given by the GCMC solutions was found to be advantageous for OMA, but it is not applicable to NOMA.  Note that, in this simple example, the global posterior covariance matrix is  ${\mathbf{C}}=( \sum_{k=1}^K {\mathbf{C}}_k^{-1})^{-1}$, and hence we use this ground-truth matrix in \eqref{eq: err2} in lieu of the estimate $\frac{1}{S_g}\sum_{s=1}^{S_g}{{\mathbf{m}} \vartheta}_g^{{(s)}}( {{\mathbf{m}} \vartheta}_g^{{(s)}})^{\sf T} $. We use full mini-batches, i.e., we implement gradient descent, such that the gradients $\nabla_{{\mathbf{m}} \vartheta} p({{\mathbf{m}} \vartheta},\mathcal{Z})$ in  \eqref{eq: gradient Dk}  and \eqref{eq: gradient Dk noma} are  computed as $\nabla_{{\mathbf{m}} \vartheta} p({{\mathbf{m}} \vartheta},\mathcal{Z})=\nabla_{{\mathbf{m}} \vartheta} p({{\mathbf{m}}{\vartheta}}|\mathcal{Z}) $.
All results are averaged over $100$ experiments.

\subsubsection{Error vs number of communication blocks}

\begin{figure}[t]
\centering
\includegraphics[width=7.5cm]{./sim_VerifyAsymp.pdf}
\caption{Second-order test error \eqref{eq: err2} versus the number of communication block, $T$, for wireless CMC ($\mathsf{SNR}=5$~dB, homogeneous subposteriors).}
\label{Fig: comm block}
\end{figure}

In Fig. \ref{Fig: comm block}, we consider homogeneous subposteriors, and plot the test error \eqref{eq: err2} versus the number of communication blocks. We recall that the number of samples is $S=T/K$ for OMA and  $S=T$ for NOMA.
WVCMC with NOMA is seen to be the most communication-efficient approach, as its test error vanishes for a number of blocks $T\simeq 7$,  while the other schemes require a large number of communication blocks. Validating the analysis in Sec. \ref{sec: WGCMC},  WGCMC, under both OMA and NOMA, is observed to attain  optimality as the number of communication blocks, and hence of samples, increases. We also note that, as mentioned in Sec. \ref{sec: WGCMC}, the conventional GCMC scheme is not asymptotically optimal with a large number of communication blocks.

\subsubsection{Error vs SNR}
\begin{figure}[t]
\centering
\includegraphics[width=7.5cm]{./sim_toy_SNR.pdf}
\caption{Second-order test error \eqref{eq: err2} versus SNR for wireless CMC ($T=2000$, heterogeneous subposteriors).}
\label{Fig: snr}
\end{figure}

We now study the impact of SNR in Fig. \ref{Fig: snr} under the assumption of heterogeneous subposteriors, i.e., subposteriors with different  covariance matrices $\{{\mathbf{C}}_k\}_{k=1}^K$ as described above.  WVCMC is seen to successfully adapt to different SNR levels, outperforming all the other schemes. In fact, the error of WVCMC is essentially unaffected by the SNR level in the considered range. This observation points to a key property of the proposed channel-driven MC sampling: If properly accounted for, the randomness injected by the channel noise may not affect the performance of CMC algorithms. Channel noise can be considered as a component of the sample generation process, rather than a nuisance to be mitigated.
Importantly, the results in Fig. \ref{Fig: snr} suggest that it is critical to properly optimize the aggregation function in order to leverage this unique property of MC sampling over noisy channels. We also note that the noise-invariant error of WVCMC is attained in the ideal case studied in this example, for which the subposterior and channel noise follow Gaussian distributions.
As two further comments on Fig. \ref{Fig: snr}, WGCMC under OMA is confirmed to be  advantageous over conventional GCMC, with the performance of two schemes converging at high SNR. Moreover, WGCMC under NOMA is seen to be worse than that of OMA, due to the performance loss counted by the heterogeneity of the subposteriors (see Sec. \ref{sec: WGCMC noma} for further discussion). Note that this is unlike the proposed WVCMC NOMA scheme, which remains advantageous over OMA.

\subsection{Bayesian Probit Regression}\label{sec: sim synthetic data}
We turn to focus on the more complex problem of Bayesian probit regression. Accordingly, given observations $\{{\mathbf{z}}_n=({\mathbf{u}}_n,v_n)\}_{n=1}^N$, the likelihood is given by the following binary regression model
\begin{align}\label{eq: probit model}
p(v_n=1|{{\mathbf{m}} \vartheta},{\mathbf{u}}_n)=\Phi({{\mathbf{m}} \vartheta}^{\sf T} {\mathbf{u}}_n),
\end{align}
where  $\Phi(\cdot)$ is the standard normal cumulative distribution function.
Furthermore, the prior is the Gaussian distribution $\mathcal{N}(0,\sigma^2 {\mathbf{I}})$ for a fixed variance $\sigma^2$.

We first consider a synthetic data set following the model described in \eqref{eq: probit model}. This is comprised of $N=8500$ labeled data points $\{{\mathbf{z}}_n=({\mathbf{u}}_n,v_n)\}_{n=1}^N$. The covariates ${\mathbf{u}}_n \in \mathbb{R}^{d}$ are drawn i.i.d. from Gaussian distribution $\mathcal{N}(0, {\mathbf{I}}_d)$ where $d=5$, and the binary label $v_n$ is drawn as per \eqref{eq: probit model} under the ground-truth model parameter ${{\mathbf{m}} \vartheta}^*=[0.1103, -0.5832, 0.6417, 1.8279, 0.4968]^{\sf T}$. {\color{black}Following \cite{rabinovich2015variational},  we also consider a higher-dimensional data set comprised by $N=15000$ data points with dimension $d=300$ generated in the same manner.}

The channels ${\mathbf{H}}_k$ are  $2d \times (2d+2)$ random matrices with i.i.d  standard Gaussian entries, and thus we have $\mathsf{E}\big[ \big({\mathbf{H}}_k {\mathbf{H}}_k^{\sf T})\big]^{-1}={\mathbf{I}}_{2d}$.  We implement analog repetition coding, i.e., we set the encoding matrices as ${\mathbf{E}}_k=\sqrt{P_k}[{\mathbf{I}}_d, {\mathbf{I}}_d]^{\sf T}$ for OMA and ${\mathbf{E}}=\min_k \sqrt{P_k}[{\mathbf{I}}_d, {\mathbf{I}}_d]^{\sf T}$ for NOMA, where  $P_k= PS/\sum_{s=1}^S2 \|{{\mathbf{m}} \vartheta}_k^{(s)}\|^2$.  Unless stated otherwise, the data set is divided equally among $K=20$ workers;
the number of communication blocks is set to $T=1000$;
the number of per-output sample iterations  to $t_m=50$ for NOMA and OMA;
 the learning rate to $\eta=10^{-7}$ for NOMA and $\eta=10^{-6}$ for OMA upon separate optimizations;
 and the initial point used for WVCMC is ${\mathbf{W}}^0=1/K\mathsf{E}_k ^\dagger$ for NOMA while in OMA it is set as GCMC solution as in the previous example.
 We use full mini-batches, i.e., $N_b=N$ in  \eqref{eq: gradient Dk}  and \eqref{eq: gradient Dk noma}, and details about gradient computation and about the implementation of Gibbs sampling from (sub)posteriors are given in Appendix \ref{appendix: implementation}. To show the effect of the partitions of local data set and the necessity of consensus, we also include the performance of the best single worker in OMA for reference. In this case, the server only uses the samples from the worker that yields the smallest test error.
Unless stated otherwise, the results  are averaged over $40$ experiments.

\subsubsection{Computational Efficiency at the Server} In this experiment, we compare different solutions as a function of the total number of  computed gradients at the server. As discussed, we are specifically interested in the comparison with SGLD. In Fig. \ref{Fig: CompEfficiency}, we plot the test error  \eqref{eq: err2} versus the number of computed gradients $\nabla_{{\mathbf{m}} \vartheta}\log p({{\mathbf{m}} \vartheta},{\mathbf{z}}_n)$, with the shaded area representing $90\%$ high-probability interval to illustrate both mean and variance during the whole training process. We recall that the number of computed gradients at the server is $NSt$ for WVCMC and $N_b t$ for SGLD where $t$ is the iteration round. Note that we highlight high-probability intervals here and not in other figures since the other plots only show the performance at convergence.

For WVCMC, we use $S=50$ samples, corresponding to the number of communication blocks  $T=1000$ for OMA and $T=50$ for NOMA, and the $\mathsf{SNR}=P/(2d N_0)$ is set as $15$ dB. For SGLD,  the initialization is randomized using a sample from the prior;  mini-batch size in \eqref{eq: update SGLD} is $N_b=500$; the learning rate is  $\eta^{(t)}= \alpha(\beta+t)^{-\gamma}$, with $\alpha=0.01$, $\beta=1$, $\gamma=0.7$; and number of burn-in samples $t_b=10^4$. All results are averaged over $100$ experiments.

The most important observation is that, by leveraging the computing power of $K$ parallel workers, WVCMC can vastly outperform SGLD for the same number of computed gradients at the server.
WVCMC with NOMA is seen to achieve the best accuracy when the number of computed gradients are larger than $5\times10^6$, and such a performance is attainable with least variance among all the schemes. The performance of WVCMC with OMA has a larger variance which is caused by the additive channel noise accumulated across multiple communication blocks due to the more inefficient orthogonal access.

\begin{figure}[t]
\centering
\includegraphics[width=7.5cm]{./sim_NumGrad_Probit.pdf}
\caption{Second-order test error \eqref{eq: err2} versus the number of computed gradients $\nabla_{{\mathbf{m}} \vartheta}\log p({{\mathbf{m}} \vartheta},{\mathbf{z}}_n)$ for wireless VCMC ($\mathsf{SNR}=15$~dB, $S=50$) and centralized SGLD.}
\label{Fig: CompEfficiency}
\end{figure}

\subsubsection{Error vs SNR} Fig. \ref{Fig: snr probit} plots the test error in \eqref{eq: err2} versus the SNR.
WVCMC is again confirmed to outperform all the other schemes. However, unlike the Gaussian example, for which the performance of WVCMC is unaffected by the SNR level, for the non-Gaussian subposteriors considered here, the performance degrades slightly in the low-SNR regime. Unlike WVCMC, the performance of WGCMC schemes is significantly decreased by the non-Gaussianity of the subposteriors.

\begin{figure}[t]
\centering
\includegraphics[width=7.5cm]{./sim_SNR_Probit.pdf}
\caption{Second-order test error \eqref{eq: err2} versus SNR for wireless CMC schemes ($T=1000$, $K=20$).}
\label{Fig: snr probit}
\end{figure}

\subsubsection{Error vs Number of workers} We now study the impact of the number of workers $K$ for $\mathsf{SNR}=35$ dB. The learning rate is set to adapt to the number of workers as $\eta=2\times10^{-5}/K$ for OMA and  $\eta=2\times10^{-6}/K$ for NOMA.
Fig. \ref{Fig: NumWorker} shows that increasing the number of workers has a negative effect on all the schemes. This degradation is, in part, caused by the reduction in the size of local data sets, and by the resulting need to combine more local samples \cite{rabinovich2015variational,scott2016bayes}. For OMA, increasing $K$ has the addition drawback of reducing the number of received samples for, which is seen to further degrade its performance for $K\geq20$.

\begin{figure}[t]
\centering
\includegraphics[width=7.5cm]{./sim_NumWorker_Probit.pdf}
\caption{Second-order test error \eqref{eq: err2} versus the number of workers for wireless CMC ($T=1000$, $\mathsf{SNR}=35$ dB).}\label{Fig: NumWorker}
\end{figure}

\subsubsection{Data Heterogeneity}
Fig. \ref{Fig: Het} plots the error versus a measure of the heterogeneity of the local data sets. To this end, the global data set is partitioned as follows. The $k$-th worker contains a fraction $\frac{1/k^{\zeta }}{\sum_{k=1}^K 1/k^{ \zeta }  } $ of data points from class `1' and a fraction $\frac{1/(K-k+1)^{\zeta}}{\sum_{k=1}^K 1/k^{ \zeta}  } $ of data points from class `0'.  The heterogeneity is described by the exponent  $\zeta$, with $\zeta=0$ representing homogeneous local data sets. As shown in Fig. \ref{Fig: Het},
increasing the heterogeneity of local data sets generally affects negatively to all schemes. Nevertheless, the heterogeneity of local data sets has a stronger impact for NOMA than for OMA, since, as discussed in Sec. \ref{sec: WGCMC noma}, NOMA ideally requires  homogeneous subposteriors to achieve optimality.

\begin{figure}[t]
\centering
\includegraphics[width=7.5cm]{./sim_Het_Probit.pdf}
\caption{Second-order test error \eqref{eq: err2} versus the skewness of data heterogeneity for wireless CMC ($T=1000$, $\mathsf{SNR}=20$ dB).}
\label{Fig: Het}
\end{figure}

{\color{black}
\subsubsection{Higher-Dimensional Data}
\begin{figure}[t]
\centering
\includegraphics[width=7.5cm]{./sim_SNR_Dim300.pdf}
\caption{Second-order test error \eqref{eq: err2} versus SNR for wireless CMC schemes on high-dimensional data ($T=10000$, $K=20$).}
\label{Fig: snr high dim}
\end{figure}

 We now consider the problem of Bayesian probit regression on higher-dimensional data, matching the scale of the experiments in \cite{rabinovich2015variational}.  The number of  iterations is selected as a function of the SNR  level  to guarantee the convergence, and the learning rate is set to $\eta=5\times10^{-9}$ for NOMA and $\eta=10^{-7}$ for OMA.  All the other settings are same as in the rest of this subsection, and the results are averaged over three realizations. Fig. \ref{Fig: snr high dim} plots the test error in \eqref{eq: err2} versus the SNR.  Confirming the results in Fig. \ref{Fig: snr probit}, WVCMC is seen to significantly outperform all the other schemes. Furthermore, as compared to the case of lower-dimensional data, OMA is seen to be generally advantageous with respect to NOMA, suggesting that the restriction to equal weights matrices in NOMA becomes more problematic as the problem dimension increases.

}

\subsection{MNIST Data Set}
We now consider the problem of Bayesian probit regression on the MNIST data set by using  $N=12665$  data points from the classes of $0$ and $1$. The original  data with dimension $d=784$ is projected to the lower dimension $d=30$ via principal component analysis (PCA). {\color{black}We note that the application of Bayesian learning to lower-dimensional data embeddings is a common approach to enhance scalability \cite{binois2020choice, nayebi2019framework, moriconi2020high, kristiadi2020being}.}
The prior distribution, channel statistics, and initial points used for learning are same as in Sec. \ref{sec: sim synthetic data}.  We use the mini-batch size $N_b=N/20$ in  \eqref{eq: gradient Dk}  and \eqref{eq: gradient Dk noma}.
Unless stated otherwise, the data set is divided equally among $K=10$ workers;
the number of communication blocks is set to $T=500$;
the number of per-output sample iterations  to $t_m=250$ for NOMA and OMA;
and  the learning rate to $\eta=2\times 10^{-7}$ for NOMA and $\eta=10^{-6}$ for OMA upon separate optimizations.

\subsubsection{Computational Efficiency at the Server} We start by validating the experiment in Fig.  \ref{Fig: CompEfficiency} on this data set. To this end, in this experiment, we use $S=50$ samples for WVCMC, corresponding to the number of communication blocks  $T=500$ for OMA and $T=50$ for NOMA, and the $\mathsf{SNR}=30$ dB.  For SGLD, mini-batch size in \eqref{eq: update SGLD} is $N_b=N/20$; the learning rate is  $\eta^{(t)}= \alpha(\beta+t)^{-\gamma}$, with $\alpha=0.01$, $\beta=1$, $\gamma=0.52$; and number of burn-in samples $t_b=8\times10^4$. All results are averaged over $100$ experiments. In Fig. \ref{Fig: CompEfficiency probit MNIST}, we plot the test error  \eqref{eq: err2} versus the number of computed gradients $\nabla_{{\mathbf{m}} \vartheta}\log p({{\mathbf{m}} \vartheta},{\mathbf{z}}_n)$, with the shaded area representing $90\%$ high-probability interval. Confirming the results in Fig. \ref{Fig: CompEfficiency}, WVCMC can vastly outperform SGLD for the same number of computed gradients at the server. However, unlike in the experiment on synthetic data set, for MNIST,  WVCMC with OMA is seen to achieve the best accuracy in terms of mean and variance performance. This is due to the heterogeneity of the subposteriors entailed by the use of real data, which degrades the performance of NOMA.

\begin{figure}[t]
\centering
\includegraphics[width=7.5cm]{./sim_NumGrad_MNIST.pdf}
\caption{Second-order test error \eqref{eq: err2} versus the number of computed gradients $\nabla_{{\mathbf{m}} \vartheta}\log p({{\mathbf{m}} \vartheta},{\mathbf{z}}_n)$ for wireless VCMC ($\mathsf{SNR}=30$~dB, $S=50$) and centralized SGLD on MNIST dataset.}
\label{Fig: CompEfficiency probit MNIST}
\end{figure}

\subsubsection{Error vs SNR}

\begin{figure}[t]
\centering
\includegraphics[width=7.5cm]{./sim_SNR_MNIST.pdf}
\caption{Second-order test error \eqref{eq: err2} versus SNR for wireless CMC schemes on MNIST dataset ($T=500$, $K=10$). }
\label{Fig: snr MNIST}
\end{figure}

Finally, Fig. \ref{Fig: snr MNIST} plots the test error in \eqref{eq: err2} versus the SNR.
WVCMC is confirmed to outperform all the other schemes, and OMA outperforms NOMA due to the heterogeneous subposteriors. Note that, for GCMC and WGCMC, the error is seen to increase in high SNR regime, suggesting that the presence of Gaussian noise in low SNR is actually advantageous to schemes that are designed under the assumption of Gaussian subposteriors.  Furthermore, the limited received samples in OMA cause the inaccurate estimation of covariance matrix, and WGCMC with OMA and GCMC is seen to perform worse than the single-worker benchmark.

{\color{black}
\subsubsection{Evaluation of Inference Performance}

\begin{figure}[t]
\centering
\includegraphics[width=7.5cm]{./sim_SNR_MNIST_KL.pdf}
\caption{KL divergence \eqref{eq: KL divergence} measuring the accuracy of the ensemble prediction \eqref{eq: pred dist1} versus SNR for wireless CMC schemes on MNIST dataset ($T=500$, $K=10$). }
\label{Fig: snr KL}
\end{figure}

We now study the inference performance of the ensemble prediction \eqref{eq: pred dist1} via the KL divergence \eqref{eq: KL divergence} evaluated on a test set of $N_t=2115$ PAC-pre-processed MNIST images. The results in Fig. \ref{Fig: snr KL} confirm the conclusions shown from Fig. \ref{Fig: snr probit} and Fig. \ref{Fig: snr high dim} in terms of the advantages with respect to all benchmarks.

}

\section{Conclusions}\label{sec: conclusions}

In this paper,  we have considered, for the first time, consensus Monte Carlo (CMC) solutions for distributed one-shot  Bayesian learning in wireless data centers.  Two main ideas were leveraged by the proposed designs. The first is the novel concept of channel-driven sampling: under uncoded wireless transmission of  samples from workers to server, channel noise can contribute to MC sampling, and, if properly accounted for,  is not necessarily a nuisance to be mitigated. In particular, we have argued that channel-driven sampling is asymptotically optimal with respect to the number of samples under the assumption that the  subposteriors are Gaussian. The second idea we build on is over-the-air computing in order to enhance spectral efficiency. This is in line with previous works on federated (frequentist) learning in wireless systems \cite{zhu2019broadband,amiri2020machine,yang2020federated,liu2020privacy}.  although the goal here is synthesizing samples from the posterior and not computing an average.
 For arbitrary subposteriors, building on these two ideas, we have proposed a general variational strategy termed wireless variational CMC (WVCMC) that aims at minimizing the KL divergence between the distributions of the output of the aggregation function at the server and of the global posterior.  Simulation results show that the proposed approach can vastly outperform existing solutions.

{\color{black}We note that  performance of the proposed WCMC schemes depends on the  aggregation function. Although only linear mappings are considered in this paper in line with prior work on the topic, the introduced variational framework can be directly applied to other functions. }
{\color{black}As  for future work, one may consider studying channel-driven sampling for digital transmission via lattice coding \cite{zamir2014lattice}.} It would also be interesting to alleviate the assumption that workers can generate samples from local subposteriors. A possible approach is to carry out approximate MC sampling at workers, e.g., distributed SGLD \cite{ahn2014distributed} or distributed SVGD \cite{kassab2020federated}.   As another possible extension of the current work, optimal analog encoding and analog compression in \cite{abdi2020analog, abdi2020quantized}, could be  tailored for wireless CMC strategies so as to further enhance the communication efficiency. {\color{black} One could also enforce a short-term power constraint in each communication block. A possible way to do this involves introducing a power control threshold, in a manner similar to \cite{cao2020optimized,liu2020privacy}. Accordingly, workers whose channel gain is too small refrain from transmitting. The optimization of the thresholds is an interesting open problem.}

redpajama_set_name: RedPajamaArXiv


Still a family-owned business, this Danforth location of Leon's Furniture Ltd. Has a huge showroom. Furnishings and décor for the whole home, including appliances. Some furniture can be customized with fabric choices. Founded in 1909 by Ablan Leon, the A. Leon Company started out as a general merchandise store in the small town of Welland in Southern Ontario. Today, Leon's is one of Canada's largest retailers selling a wide range of merchandise including furniture, major appliances, and home electronics.


text:
Eugen Wirth (* 12. Mai 1925 in Würzburg; † 15. Mai 2012 in Erlangen) war ein deutscher Geograph. Sein Hauptforschungsgebiet war die Kulturgeographie des Orients.

Leben
Wirth studierte von 1945 bis 1952 an den Universitäten in Erlangen, Göttingen, Freiburg im Breisgau und an der Sorbonne in Paris. Er wurde 1952 an der Universität Freiburg in den Fachbereichen Soziologie und Philosophie promoviert und 1959 an der Universität Hamburg unter Erich Otremba in Geographie habilitiert. Von 1964 bis 1991 leitete Wirth das Institut für Geographie an der Universität Erlangen. Von 1980 bis 1993 war Wirth Vorsitzender der Fränkischen Geographischen Gesellschaft. Außerdem war er Fachgutachter für Anthropogeographie der Deutschen Forschungsgemeinschaft. Wirth war evangelisch und seit 1958 mit Ingeborg, geb. Holland, verheiratet. Sie hatten vier gemeinsame Kinder (zwei Töchter, zwei Söhne).

Ehrungen
Wirth war Mitglied der Academia Europaea (1990) in London und des Deutschen Archäologischen Instituts sowie der Österreichischen Akademie der Wissenschaften. 1996 verlieh ihm die Goethe-Universität Frankfurt die Ehrendoktorwürde.

Veröffentlichungen (Auswahl)
 Agrargeographie des Irak. (Hamburger geographische Studien 13) Institut für Geographie und Wirtschaftsgeographie der Universität Hamburg, Hamburg 1962 (= Habilitationsschrift)
 Syrien, eine geographische Landeskunde. Wissenschaftliche Buchgesellschaft, Darmstadt 1971
 Zum Problem des Basars (Suq, Carsi). Versuch einer Begriffsbestimmung und Theorie des traditionellen Wirtschaftszentrums der orientalisch-islamischen Stadt. In: Islam, Band 52, Heft 1, 1975, S. 203–260
 Theoretische Geographie: Grundzüge einer Theoretischen Kulturgeographie. Teubner, Stuttgart 1979.
 Zur Konzeption der islamischen Stadt: Privatheit im islamischen Orient versus Öffentlichkeit in Antike und Okzident. In: Die Welt des Islams, Band 31, Heft 1, 1991, S. 50–92
 Handlungstheorie als Königsweg einer modernen Regionalen Geographie – Was 30 Jahre Diskussion um die Länderkunde gebracht haben. In: Geographische Rundschau, Band 51, Heft 1, 1999, S. 57–64
 Die orientalische Stadt im islamischen Vorderasien und Nordafrika. Städtische Bausubstanz und räumliche Ordnung, Wirtschaftsleben und soziale Organisation. 2 Bände. Philipp von Zabern, Mainz 2000, ISBN 3-8053-2709-9.

Literatur

Einzelnachweise

Geograph (20. Jahrhundert)
Hochschullehrer (Friedrich-Alexander-Universität Erlangen-Nürnberg)
Ehrendoktor der Goethe-Universität Frankfurt am Main
Mitglied des Deutschen Archäologischen Instituts
Mitglied der Österreichischen Akademie der Wissenschaften
Mitglied der Academia Europaea
Deutscher
Geboren 1925
Gestorben 2012
Mann

redpajama_set_name: RedPajamaWikipedia


Performance audits are used to assess the monitoring system (including the on site monitoring station and laboratory) and the quality at which it is functioning. Auditing personnel will test and grade the monitoring instruments, standard gases, sample analysis techniques, and data processing procedures. This is a quantitative assessment task and auditing methods will differ with various instruments and equipment.
Standard gases of various densities will be injected into the analysis instrument during audit with density range being three to four different densities between zero and the maximum range on the instrument. After the standard gases are diluted they will be injected into the analysis instrument and compared to the numbers of the known densities to calculate the deviation percentage.
Also, the standard gas ozone during audit is taken from the O3 generator in the calibrator which produces ozones of various densities which are injected into a certified analyzer used as a transfer standard which is taken to the audited monitoring station and compared to the ozone analysis instrument in the station to assess the precision of the ozone analysis instrument in the station.
The standard gases used by this administration during audits include CO, NO, SO2, C3H8 in high density gas cylinders which are diluted through calibrators into standard gases of various densities. The standard gases listed above must follow the highest standards such as those listed by America's National Institute of Standards and Technology (NIST) and their standard reference materials (SRM) or the their EPA's national traceable reference materials (NTRM). Also, gases produced by the gas manufacturers intermediate standard (GMIS) can also be used as standard materials.
Our testing stations include analysis instruments for sulfur nitrate (SO2), nitrogen oxides (Nox, NO, NO2), ozone (O3), carbon monoxide (CO), and hydrocarbons (HC) and performance audits are held once annually (table 1). Table 2 shows the accepted audit result standards of the administration --- automatic monitoring facilities data quality goals have been amended according to foreign standards, our monitoring instrument specifications, and opinions from scholars and specialists.
Gas dilution calibration system (air) ±5% ≧0.9950 N/A ±3%F.S.
Gas dilution calibration system (gas) ±5% ≧0.9950 N/A ±3%F.S.
Sulfur nitrate (SO2) ±15% ≧0.9950 0.85-1.15 ±3%F.S.
Nitrogen oxide (Nox) ±15% ≧0.9950 0.85-1.15 ±3%F.S.
Carbon monoxide (CO) ±15% ≧0.9950 0.85-1.15 ±3%F.S.
Ozone (O3) ±15% ≧0.9950 0.85-1.15 ±3%F.S.
Total hydrocarbons (THC) ±15% ≧0.9950 0.85-1.15 ±3%F.S.
The auditing of suspended micro particle samplers (including TSP, highvol PM10, continuous PM10) includes two audits: one is the audit of actual flow value and its deviation with the value displayed on the instrument, and the other is an audit of how close the actual flow values are to the designed flow values in the instruments. Average commercial PM10 highvol and TSP samplers have a flow design of 40.0 CFM (cubic feet per minute), but actual flow allowed is 40.0± 10%CFM. Continuous automatic monitor PM10 samplers are designed with a flow amount of 18.9LPM (Wedding) or 16.7LPM (Anderson). Thus, the standard value when auditing is either 18.9± 10%LPM or 16.7± 10%LPM.
The main audit instrument for non-automatic monitoring PM10 highvol and TSPs are the root meter, orifice, and the difference in pressure of the water pillar in correcting curve graphs to find true flow which is then compared to the recorded flow on the sampler to calculate the deviation percentage. As automatic monitoring PM10 analysis instruments use foam flow as their transfer standard, the true flow is obtained from the station to calculate the deviation percentage with displayed flow.


Today marks the day after Thanksgiving and from this point on, it will just be Dramas I am thankful for. Ones that have impacted me or left their mark in me in some way, shape, or form. It's a series within a series. Series-ception!!!!!!!!!!!!
This drama was one of the best dramas to date. I've done a drama series before and this one will forever always pop up.
I'm sure Joong Ki and Hye Kyo are extremely thankful for the drama too.
Such a cute married couple.
This was the first drama where I really got to see Joong Ki shine.
And that's something to really be thankful for.
There's not too much I can say on this blog that I didn't say on the other one. Just wanted everyone to know I just appreciate this drama.


Nintendo has released a free update for the Nintendo Switch version of Captain Toad: Treasure Tracker which adds a new two-player co-op mode.
Guess who is still ready for adventure? After downloading a free Nintendo Switch version update for the Captain Toad: Treasure Tracker game that is available starting today, all courses in the game will support a new two-player co-op mode featuring a second Toad to explore simultaneously with a friend.
And on March 14, new paid DLC arrives that adds 18 new challenges to the game, including five new courses and new Shiny Crowns, which can be obtained in each course.
New players can get the Captain Toad: Treasure Tracker game together with the updates and all the new Special Episode content in a digital bundle available for purchase on Nintendo eShop beginning February 13.
Players who purchase the bundle or the Captain Toad: Treasure Tracker – Special Episode downloadable content will receive early access to one of the courses in Captain Toad: Treasure Tracker – Special Episode, with the remainder of the courses launching on March 14.