\section{Introduction}
\label{Sec:introduction}

An inverse problem consists in inferring an unknown variable $\overline{\xb}\in {\ensuremath{\boldsymbol{\mathsf H}}} $ from degraded observations $\yb\in {\ensuremath{\boldsymbol{\mathsf H}}} $.
In the remainder, we will consider that ${\ensuremath{\boldsymbol{\mathsf H}}} = \ensuremath{\mathbb{R}}^{\overline{N}}$ and ${\ensuremath{\boldsymbol{\mathsf G}}} = \ensuremath{\mathbb{R}}^{\overline{M}}$.
Observations and parameters are typically related by a model of the form
\begin{equation} \label{eq:full_model}
    \yb = \ensuremath{\mathcal{D}}(\Ab \overline{\xb}),
\end{equation}
where the linear operator $\Ab $ models the acquisition process, and $\ensuremath{\mathcal{D}}$ models random perturbations damaging the observations. Bayesian inference relies on the posterior distribution of the random variable $\xb$ to estimate the true value $\overline{\xb}$. The posterior distribution combines information from the likelihood -- related to the observations $\yb$ -- and the prior. For instance, in image processing, a usual choice consists in promoting sparsity in a selected basis, e.g., a gradient basis leading to the total variation (TV) regularization~\cite{Rudin1992}, or a wavelet basis~\cite{Mallat2009}. Prior information can also encompass constraints such as nonnegativity, e.g., for intensity images~\cite{Cai2018, Thouvenin2021}, or more evolved constraints based on the physics of the data acquisition process, e.g., polarization constraints~\cite{Birdi2018}.
This work focuses on inverse problems whose statistical description is encoded by a generic posterior distribution of the form
\begin{equation} \label{eq:dist-gen-lin}
    \dist{}{\xb} \propto \exp ( - h(\xb) -f(\xb) - g(\Db \xb) ),
\end{equation}
where $h \colon {\ensuremath{\boldsymbol{\mathsf H}}} \to ]-\infty,+\infty]$ is a Lipschitz-differentiable function, $f \colon {\ensuremath{\boldsymbol{\mathsf H}}} \to \ensuremath{\left]-\infty,+\infty\right]}$, $g\colon {\ensuremath{\boldsymbol{\mathsf G}}} \to ]-\infty,+\infty]$, and $\Db \colon {\ensuremath{\boldsymbol{\mathsf H}}} \to {\ensuremath{\boldsymbol{\mathsf G}}}$ is a linear operator.
A classical inference approach is to form the maximum \emph{a posteriori} (MAP) estimator, defined as a solution to
\begin{equation}	\label{eq:min-gen-lin}
	\minimize{\xb \in {\ensuremath{\boldsymbol{\mathsf H}}}} h(\xb) +f(\xb) + g(\Db \xb).
\end{equation}
The problem can be efficiently solved with proximal primal-dual methods \cite{Komodakis2015, Condat2013, Vu2013}. These approaches recast the problem~\eqref{eq:min-gen-lin} as a saddle point problem in order to simplify computations associated with the composite function $g\circ \Db$.
Splitting proximal algorithms benefit from many acceleration techniques (e.g., inertia \cite{Beck2009, Nesterov1983}, preconditioning \cite{Chouzenoux2014, Combettes_Vu_2014}), they are versatile, and scalable. In particular, they are highly parallelizable, and can be efficiently distributed to split the computational cost per iteration, under well-established theoretical guaranties~\cite{Bauschke2017CAMOTHS, Combettes2011, Komodakis2015}.
However, these algorithms only provide a point estimate, without additional information.
In absence of ground truth, these approaches do not allow the uncertainty about the estimate to be directly quantified.

Markov chain Monte Carlo (MCMC) algorithms are generic approaches providing estimates with associated credibility intervals~\cite{Robert2010}. They aim to generate a Markov chain that yields samples from the target posterior distribution~\eqref{eq:dist-gen-lin} in the stationary regime. Nevertheless, they are often considered computationally too expensive to handle high dimensional problems, especially when composite functions are involved.
Over the last decade, multiple authors have proposed more versatile and scalable optimization-inspired MCMC algorithms~\cite{Durmus2018, Pereyra2016,Salim_Richtarik_2020}.
These approaches exploit quantities repeatedly used in optimization, such as gradients and proximity operators, to efficiently explore high dimensional parameter spaces.
In particular, a splitting approach based on an \emph{asymptotically exact data augmentation} (AXDA) scheme has recently been proposed in~\cite{Vono_etal_2019, Vono_etal_2020}. Reminiscent of primal-dual optimization approaches, AXDA introduces auxiliary variables to split composite distributions. The distribution with density~\eqref{eq:dist-gen-lin} is approximated by
\begin{equation} \label{eq:dist-SPA}
    \dist{(\alpha, \beta)}{\xb,\zb,\ub}
    \propto \exp \big( - h(\xb) -f(\xb) -g(\zb) - \phi_\alpha(\Db\xb,\zb-\ub) - \psi_\beta(\ub) \big),
\end{equation}
where $\phi_{\alpha} \colon {\ensuremath{\boldsymbol{\mathsf G}}} \times {\ensuremath{\boldsymbol{\mathsf G}}} \to ]-\infty, +\infty]$, $\psi_{\beta} \colon {\ensuremath{\boldsymbol{\mathsf G}}} \to ]-\infty, +\infty]$, and $(\alpha, \beta)$ controls the discrepancy between $\Db\xb$ and $\zb$. For appropriate choices of $\phi_\alpha$ and $\psi_\beta$, the marginal distribution of $\xb$ with respect to~\eqref{eq:dist-SPA} converges to the target distribution \eqref{eq:dist-gen-lin} as $(\alpha,\beta) \to (0,0)$. A block Gibbs sampler is proposed in~\cite{Vono_etal_2019, Vono_etal_2020} to draw samples from~\eqref{eq:dist-SPA}, referred to as the Split-and-Augment (SPA) Gibbs sampler.

Designing efficient algorithms to handle distributions of the form~\eqref{eq:dist-gen-lin} becomes more challenging when multiple composite functions are considered, i.e.,
\begin{equation}    \label{eq:g-multi}
    (\forall \xb \in {\ensuremath{\boldsymbol{\mathsf H}}})\quad
    g(\Db \xb) = \sum_{i=1}^I g_i(\Db_i \xb),
\end{equation}
and when the dimension of $\xb$ and $\yb$ increases.
On the one hand, primal-dual proximal algorithms~\cite{Condat2013, Vu2013, Komodakis2015, Pesquet2015, Combettes2018} enable composite terms~\eqref{eq:g-multi} to be handled in parallel in the dual domain. In this context, the dual variables can be distributed on multiple workers. This strategy has also been adopted in \cite{Vono_etal_2019} to design an SPA Gibbs sampler tailored to distributions involving composite terms~\eqref{eq:g-multi}.
On the other hand, to address problems with high dimensional variable $\xb$, a usual strategy in optimization consists in splitting $\xb$ into blocks $(\xb_k)_{1 \le k \le K}$ and either alternate between the blocks \cite{Bolte2014, Chouzenoux2016, Tseng92, TsengBCD2001}, or distribute the blocks over multiple workers \cite{Combettes2015, Pesquet2015}.
In particular, in \cite{Pesquet2015}, the authors combine such block-coordinate approaches with primal-dual splitting techniques to parallelize and distribute both the primal and dual variables.
Then, the associated minimisation problem is of the form
\begin{equation}	\label{eq:min-gen-multi}
	\minimize{\xb =(\xb_k)_{1\le k \le K} \in {\ensuremath{\boldsymbol{\mathsf H}}}} \sum_{k=1}^K h_k(\xb_k) +f_k(\xb_k) + \sum_{i=1}^I g_i(\sum_{k=1}^K \Db_{i,k} \xb_k),
\end{equation}
where, for every $k\in \{1, \ldots, K\}$, $h_k$ and $f_k$ only act on the $k$-th block of the variable $\xb \in {\ensuremath{\boldsymbol{\mathsf H}}}$.
In \cite{Pesquet2015}, the authors also combine the resulting algorithms with consensus constraints, imposing that all $(\xb_k)_{1 \le k \le K}$ are equal, to develop asynchronous distributed algorithms over hypergraph structures, allowing high flexibility in the choice of the distribution architecture.

Only a few distributed samplers have been proposed in the literature~\cite{Rendell2021, Vono_etal_2019, Vono_etal_2020}.
However, these distributed samplers only focus on distributing the dual variables associated with the composite functions \eqref{eq:g-multi}, without splitting the global variable into blocks. In particular, in~\cite{Rendell2021}, a consensus-based approximate posterior distribution has been proposed, addressed with a distributed Metropolis-within-Gibbs sampler relying on a \emph{client-server} architecture. The sampler exploits the conditional independence between some variables to parallelize computations. This is especially relevant for data-distributed applications, in which a shared parameter value needs to be inferred from a dataset distributed over multiple workers.
However, this setting has several drawbacks.
First, the consensus constraint necessitates the variables of interest to be duplicated on the different workers, which can be a critical issue in a high dimensional setting.
Second, the client-server architecture may induce communication bottlenecks, as all (\emph{clients}) workers need to communicate with the server.
It also exhibits limitations in terms of distribution flexibility, since it only applies to cases with more than two composite functions, without splitting the global variable of interest $\xb$ into blocks.

\begin{figure}[!b]
\centering
\includegraphics[width=0.5\textwidth,keepaspectratio]{figures/illustrations/blocksparse.png}
\vspace*{-0.3cm}
\caption{\label{fig:blocksparse}\small
Example of the block-sparse structure of the matrix $\Db$ involved in~\eqref{eq:dist-gen-lin}.
The columns of $\Db$ are split into $K$ contiguous blocks (black dashed lines), with a small overlap compared to the size of the blocks.
The orange dashed rectangles highlight the subparts of $\Db$ implemented on each worker $k$. These act on a parameter block stored on worker $k$, and a few parameters stored on worker $k+1$.
}
\end{figure}

In this work, we introduce a distributed block-coordinate SPA-Gibbs sampler, splitting the global variable of interest into blocks $(\xb_k)_{1 \le k \le K}$, and able to handle multiple composite terms of the form~\eqref{eq:g-multi}. Each block $\xb_k$ is associated with a worker, and communications between the workers are characterized by the structure of the linear operators $(\Db_i)_{1 \le i \le I}$. Precisely, the communications are characterized by the hyperedges of hypergraphs whose adjacency matrices are prescribed by each of the linear operators.
For the sake of clarity and simplicity of the presentation, the main part of the paper presents the proposed distributed algorithm for the probability distribution~\eqref{eq:dist-gen-lin}. The proposed approach is then extended in the appendix to distributions involving multiple linear operators, as in \eqref{eq:g-multi}.
Note that the hypergraph structure of the proposed algorithm is different from the one used in~\cite{Pesquet2015}. Indeed, we propose to use the hypergraph to distribute the communications induced by $\Db$ with respect to both the primal and dual variable, while in~\cite{Pesquet2015}, the hypergraph structure is used when multiple $(\Db_i)_{1 \le i \le I}$ are considered without splitting the primal variable into blocks $(\xb_k)_{1 \le k \le K}$.

The resulting proposed algorithm accommodates several distributed architectures, and is particularly suitable for a Single Program Multiple Data (SPMD) architecture~\cite{Darema2001}. In contrast with a client-server configuration, all the workers involved in an SPMD architecture execute similar tasks on a subpart of all the variables. An SPMD architecture can drastically reduce the communication costs compared to a client-server architecture when a small number of workers are involved in the communications, and especially when the hypergraph structure of $\Db$ induces \emph{localized} couplings between parameters. A particular case consists in considering a block-sparse matrix $\Db$, as illustrated in Figure~\ref{fig:blocksparse}. This case is encountered in applications such as image deconvolution or inpainting, or when considering local priors such as the TV norm.
When multiple composite terms are considered as in~\eqref{eq:g-multi}, the proposed algorithm is also amenable to a parallel implementation based on a client-server architecture, exploiting conditional independence between blocks of variables. As previously emphasized, this configuration may however suffer from  communication bottlenecks, with a number of workers limited to the number of conditionally independent blocks of variables in the model. In contrast, the SPMD approach allows a larger number of workers to be used by exploiting the structure of the hypergraphs.

The remainder of the paper is organized as follows.  The AXDA approach and SPA-Gibbs sampler~\cite{Vono_etal_2019, Vono_etal_2020} are recalled in~\Cref{Sec:SOA}. The proposed hypergraph model is introduced in \Cref{Sec:Model}, and the associated distributed block-coordinate SPA-Gibbs sampler is given in \Cref{Sec:dist-single-algo}. \Cref{Sec:application} describes an SPMD implementation of the proposed method for a large scale image deconvolution problem. Conclusion and perspectives are given in~\Cref{Sec:conclusion}.
Eventually, in~\Cref{appendix:SPA-PSGLA-multi} the proposed method is extended to distributions involving multiple composite terms.

\section{AXDA and SPA Gibbs samplers}
\label{Sec:SOA}

In this section, we summarize the approximation results at the basis of the AXDA approach~\cite{Vono_etal_2019}. We also describe the SPA Gibbs sampler to approximately draw samples from the distribution~\eqref{eq:dist-gen-lin}.

\subsection{AXDA splitting approach}
\label{Ssec:SPA-simple}

Consider the distribution $\dist{(\alpha, \beta)}{\xb,\zb,\ub}$ introduced in~\eqref{eq:dist-SPA}. Under technical conditions on the functions $\phi_\alpha$ and $\psi_\beta$, \cite{Vono_etal_2019} showed that the marginal distribution of $\xb$ can be made arbitrarily close to the original distribution $\dist{}{\xb}$ given in~\eqref{eq:dist-gen-lin}, depending on the choice of $(\alpha, \beta)$. This result is summarized in the following proposition.
\begin{proposition} \label{Prop:cv-SPA-simple}\cite[Thm. 1 and 2]{Vono_etal_2019}
Let $(\alpha ,\beta) \in \ensuremath{\left]0,+\infty\right[}^2$, $\phi_{\alpha} \colon {\ensuremath{\boldsymbol{\mathsf G}}} \times {\ensuremath{\boldsymbol{\mathsf G}}} \to ]-\infty, +\infty]$ and $\psi_{\beta} \colon {\ensuremath{\boldsymbol{\mathsf G}}} \to ]-\infty, +\infty]$.
Assume that
\begin{equation} \label{prop:cv-spa-simple-cond1}
    (\forall \xb \in {\ensuremath{\boldsymbol{\mathsf H}}}) \quad \lim_{\alpha \rightarrow 0} \frac{\exp\big(- \phi_{\alpha}(\Db \xb, \zb)\big)}{\int_{{\ensuremath{\boldsymbol{\mathsf G}}}} \exp\big(- \phi_{\alpha}( \Db \xb, \zb)\big) d\zb} = \delta_{\xb} (\zb),
\end{equation}
and that, there exists $\eta_{\alpha,\beta} \in \ensuremath{\left]0,+\infty\right[}$,
\begin{equation}\label{prop:cv-spa-simple-cond2}
    (\forall (\xb, \zb) \in {\ensuremath{\boldsymbol{\mathsf H}}} \times {\ensuremath{\boldsymbol{\mathsf G}}}) \quad
    \int_{{\ensuremath{\boldsymbol{\mathsf G}}}} \exp\big(- \phi_{\alpha}(\Db \xb, \zb-\ub) - \psi_{\beta}(\ub)  \big) d\ub
    \propto \exp\big(-\phi_{\eta_{\alpha,\beta}}(\Db \xb, \zb)\big).
\end{equation}
Let, for $\xb \in {\ensuremath{\boldsymbol{\mathsf H}}}$, $\pi_{\eta_{\alpha,\beta}}(\xb) = \int_{{\ensuremath{\boldsymbol{\mathsf G}}} \times {\ensuremath{\boldsymbol{\mathsf G}}}}  \pi_{\alpha, \beta} (\xb, \zb, \ub) d\ub d\zb$.
Then, $\|\pi - \pi_{\eta_{\alpha, \beta}}\|_{\text{TV}} \to 0$ when $\eta_{\alpha,\beta} \to 0$.
\end{proposition}

\Cref{Prop:cv-SPA-simple} ensures that the marginal distribution of $\xb$ with respect to $\pi_{(\alpha, \beta)}$ converges to the target distribution $\pi$ given in~\eqref{eq:dist-gen-lin}, provided conditions~\eqref{prop:cv-spa-simple-cond1} and \eqref{prop:cv-spa-simple-cond2} hold.
In particular, in \cite{Vono_etal_2019}, these conditions are shown to be satisfied for $\eta_{\alpha, \beta}^2 = \alpha^2 + \beta^2$, with
\begin{equation}    \label{eq:phi-psi-def}
    \begin{cases}
    \phi_\alpha(\Db \xb, \zb-\ub ) = \frac{1}{2 \alpha^2} \| \Db\xb - (\zb-\ub) \|^2, \\
    \psi_\beta(\ub) = \frac{1}{2 \beta^2} \| \ub \|^2.
    \end{cases}
\end{equation}
In the particular case when $f\equiv 0$, stronger results can be found in \cite[Theorem 2]{Vono_etal_2020}, including theoretical guarantees on the convergence rate.

A Gibbs sampler can be used to draw samples from \eqref{eq:dist-SPA}, approximating the target distribution~\eqref{eq:dist-gen-lin}.
The sampler successively draws one sample from each conditional distribution

\begingroup
\allowdisplaybreaks
\begin{align}
    \distc{\alpha}{\xb}{\zb,\ub}
    & \propto \exp \big( -h(\xb) -f(\xb)  - \phi_\alpha(\Db\xb,\zb-\ub) \big) ,    \label{eq:dist-SPA-cond:x}\\
    \distc{\alpha}{\zb }{ \vb,\ub }
    & \propto \exp \big( -g(\zb) - \phi_\alpha(\vb,\zb-\ub)  \big) ,    \label{eq:dist-SPA-cond:z}\\
    \distc{(\alpha,\beta)}{\ub }{ \vb,\zb }
    & \propto \exp \big( - \phi_\alpha(\vb,\zb-\ub) - \psi_\beta(\ub) \big) , \label{eq:dist-SPA-cond:u}
\end{align}
\endgroup
where $\vb = \Db \xb$.
The SPA Gibbs sampler associated with~\eqref{eq:dist-SPA-cond:x}--\eqref{eq:dist-SPA-cond:u} is thus given by
\begin{equation}    \label{algo:SPA-simple}
\begin{array}{l}
\text{let } \zb^{(0)} \in {\ensuremath{\boldsymbol{\mathsf G}}} \text{ and } \ub^{(0)} \in {\ensuremath{\boldsymbol{\mathsf G}}} \\
\text{for } t = 0,1, \ldots, T	\\
\left\lfloor
\begin{array}{l}
    \xb^{(t)} \sim \distc{\alpha}{\xb }{ \zb^{(t)},\ub^{(t)} }, \\
    \vb^{(t)} = \Db \xb^{(t)}, \\
    \zb^{(t+1)} \sim \distc{\alpha}{\zb }{ \vb^{(t)}, \ub^{(t)} }, \\
    \ub^{(t+1)} \sim \distc{(\alpha,\beta)}{\ub}{ \vb^{(t)}, \zb^{(t+1)} }.
\end{array}
\right.
\end{array}
\end{equation}

\subsection{PSGLA within SPA Gibbs sampler}
\label{Ssec:SPA-PSGLA}

Drawing samples directly from the conditional distributions \eqref{eq:dist-SPA-cond:x}--\eqref{eq:dist-SPA-cond:u} can still be difficult. This is the case for the application described in \Cref{Sec:application}.
To overcome this issue, samples can be drawn using appropriate transition kernels.

Let $\mathcal{K}_{ \xb,\alpha} (\xb, \zb, \ub)$, $\mathcal{K}_{\zb,\alpha} (\Db \xb, \zb, \ub)$ and $\mathcal{K}_{\ub, (\alpha,\beta)} ( \Db \xb, \zb, \ub)$
be transition kernels with invariant distributions $\distc{\alpha}{\xb}{\zb,\ub}$, $\distc{\alpha}{\zb }{ \xb, \ub }$, and $\distc{(\alpha,\beta)}{\ub}{ \xb, \zb }$, respectively (resp.).
A Metropolis-Hastings transition kernel is a classical choice for $\ensuremath{\mathcal{K}}_{\xb, \alpha}$, $\ensuremath{\mathcal{K}}_{\zb, \alpha}$ and $\ensuremath{\mathcal{K}}_{\ub, \alpha, \beta}$ to exactly draw samples from~\eqref{eq:dist-SPA}, leading to a conditional Metropolis-Hastings sampler~\cite{Jones2014}. Appropriate proposals include Langevin-based kernels, such as the Moreau-Yosida unadjusted Langevin algorithm (MYULA)~\cite{Durmus2018} and the proximal stochastic gradient Langevin algorithm (PSGLA) \cite{Salim_Richtarik_2020}. These kernels can handle differentiable and non-differentiable potential functions simultaneously, and have been shown suitable to address high dimensional problems. To avoid the extra cost of the accept-reject step, approximate sampling can be considered by using unadjusted kernels, as in~\cite{Durmus2018} (\emph{i.e.}, omitting the accept-reject step). This approach is adopted in the following. Deriving non-asymptotic convergence bounds and analyzing the bias between the target distribution and~\eqref{eq:dist-SPA} is, however, beyond the scope of this paper.

The choice of suitable transition kernels will be instrumental to design a distributed Gibbs sampler when direct sampling from the conditional distributions is difficult.
Technical assumptions on the functions in~\eqref{eq:dist-gen-lin}, described in \Cref{Sec:Model}, will ensure the transition kernels to be amenable to a distributed implementation, adapted to the specific structure of the distributions $\distc{\alpha}{\xb}{ \zb,\ub}$, $\distc{\alpha}{\zb}{ \xb,\ub}$ and $\distc{\alpha,\beta}{\ub}{ \xb,\zb}$.
In particular, for the proposed distributed SPA-Gibbs sampler, we consider the case when any transition kernel can be used to sample $(\zb, \ub)$, while $\xb$ is sampled using one PSGLA iteration. To this aim, we assume that $h$ is $\lambda_h$-Lipschitz-differentiable. In addition, for the sake of simplicity, $\phi_\alpha$ and $\psi_\beta$ are assumed to be $\ell_2$ norms as in~\eqref{eq:phi-psi-def}, so that $\ub$ can be directly sampled from its conditional distribution (i.e., Gaussian distribution).
In this case, \eqref{algo:SPA-simple} reads
\begin{equation}    \label{algo:SPA-simple-sampling}
\begin{array}{l}
\text{let } \xb^{(0)} \in {\ensuremath{\boldsymbol{\mathsf H}}}, \vb^{(0)} = \Db \xb^{(0)}, (\zb^{(0)}, \ub^{(0)}) \in {\ensuremath{\boldsymbol{\mathsf G}}}^2, \text{ and } \gamma \in \Big]0, \big(\lambda_h + \| \Db \|^2/\alpha^2\big)^{-1} \Big[,   \\
    \text{for } t = 0,1, \ldots, T	\\
    \left\lfloor
    \begin{array}{l}
        \xb^{(t+1)} = \ensuremath{\text{\rm prox}}_{\gamma f} \Bigg( \xb^{(t)} - \gamma \nabla h(\xb^{(t)}) - \gamma \Db^* \Big( \nabla \phi_\alpha( \cdot , \zb^{(t)}-\ub^{(t)}) (\vb^{(t)}) \Big) + \sqrt{2\gamma} \, \wb^{(t)} \Bigg), \\
        \vb^{(t+1)} = \Db \xb^{(t+1)}, \\
        \zb^{(t+1)} \sim \trans{\zb,\alpha}{\vb^{(t+1)}, \zb^{(t)}, \ub^{(t)} }, \\
        \ub^{(t+1)} \sim \ensuremath{\mathcal{N}} \Big( \frac{\beta^2}{\alpha^2+\beta^2} (\zb^{(t+1)} - \vb^{(t+1)}), \frac{\alpha^2+\beta^2}{\alpha^2 \beta^2} \ensuremath{\mathbf{I}}\, \Big),

    \end{array}
    \right.
\end{array}
\end{equation}
where $\Db^*$ denotes the adjoint of $\Db$, and $(\wb^{(t)})_{1 \le t \le T}$ is a sequence of independent and identically distributed (i.i.d) standard Gaussian random variables in~${\ensuremath{\boldsymbol{\mathsf H}}} $.

\section{Proposed distributed Gibbs sampler in a nutshell}
\label{Ssec:Back:dist}

This work aims to design an efficient distributed version of the SPA Gibbs sampler~\eqref{algo:SPA-simple-sampling} over multiple workers by exploiting the hypergraph structure of $\Db$.
In this section, we give a first overview of the main ideas underlying the proposed distributed sampler. The full splitting model is described in \Cref{Sec:Model}, and the associated distributed Gibbs sampler is given in \Cref{Sec:dist-single-algo}.

\subsection{A hypergraph structure to better communicate}

We consider a localized linear operator $\Db$ in~\eqref{eq:dist-gen-lin}, meaning that couplings between latent parameters are localized. This is for instance the case when $\Db$ is block-sparse, see \Cref{fig:blocksparse}.
Such structure is not necessary to design the proposed distributed algorithm, but will be instrumental in practice to reduce the communications between workers for an efficient implementation of the algorithm.
Linear operators $\Db$ with block-sparse structures are encountered in multiple inverse imaging problems, corresponding to, e.g., a deconvolution or inpainting measurement operator. Block-sparse operators also appear when considering local priors such as the TV norm.

The structure of the matrix $\Db$ can be described with a binary matrix that can be interpreted as the adjacency matrix of a hypergraph.
A hypergraph is a generalization of a graph, in which an edge can join any finite number of vertices. Formally, an undirected hypergraph $\ensuremath{\mathcal{H}}$ is a pair $(\xb, \eb)$ made of vertices $\xb = (x_n)_{1 \le n \le N}$ and hyperedges $\eb = (e_m)_{1\le m \le M}$, where each hyperedge is a set of connected vertices.

The main idea is to distribute over $K$ workers the computation of various quantities in algorithm~\eqref{algo:SPA-simple-sampling}, related to distribution~\eqref{eq:dist-SPA}, and in particular the product $\Db\xb$.
To this aim, the variable $\xb$ is divided into $K$ blocks $(\xb_k)_{1 \le k \le K}$, each assigned to a single worker.
The hyperedges will characterize the necessary communications between the workers. A worker $k \in \{1, \ldots, K\}$, that stores $\xb_k$, may need to access some coefficients of $\xb_{k'}$ stored on another worker $k'\in \{1, \ldots, K\} \setminus \{k\}$ to carry out its computations.
Note that the operator $\Db=(D_{m,n})_{1 \le m \le M, \, 1 \le n \le N}$ will be split once and for all in an adequate manner over the $K$ workers as well, but will not need to be communicated between the workers.

Similarly to the vertices, the quantity $\vb = (v_m)_{1 \le m \le M} = \Db\xb $ computed in \eqref{algo:SPA-simple-sampling} will also be distributed among the $K$ workers.
To this aim, each hyperedge $m\in \{1, \ldots, M\}$ is associated to a worker denoted by $k_m \in \{1, \ldots, K\}$.
Then, each subpart $v_m$ will be computed and stored on worker $k_m$.
In practice, only the vertices $x_n$ that correspond to non-zero blocks $D_{m, n}$ in $\Db$ will be necessary to compute $v_m$.
In addition, only a subset of $x_n$ will be used to compute $v_m$, but will not be used to compute any other $v_{m'}$. This is illustrated in \Cref{fig:blocksparse}.
Finally, the worker $k=k_m$ that is in charge of the computation of $v_m$ will also store $\xb_k$ to reduce communications.

Communications between workers will be necessary when the worker $k_m$ requires a subpart of $\xb_{k'}$ stored on another worker $k'\neq k_m$ to carry out the computation of $v_m$.
In the particular case when $\Db$ is block-sparse, the communication cost will remain small as long as the subpart of $\xb_{k'}$ required by worker $k_m$ remains small.

\subsection{Conditions for an efficient distributed implementation}
\label{Sssec:Model:cond-eff-dist}

As previously mentioned, although the proposed approach detailed in \Cref{Sec:Model} and \Cref{Sec:dist-single-algo} is very general, some particular structures of matrix $\Db$ can ensure limited communications, which are often the bottleneck of distributed methods.

Computing costs of algorithms can roughly be divided into three terms. First, a computation term $\tau_{\texttt{flop}}$ reflects the time to perform a single floating point operation. Second, a communication term $\tau_{\texttt{bandwidth}}$, defined as the inverse of the communication bandwidth, quantifies the time necessary to send a single value. Third, a latency term $\tau_{\texttt{latency}}$ represents the cost incurred by establishing a communication.
In practice, $\tau_{\texttt{flop}} \ll \tau_{\texttt{bandwidth}} \ll \tau_{\texttt{latency}}$. This implies that the number of communications and the size of the messages need to be sufficiently small to ensure the overall communication time to be negligible compared to the computation time.

Consequently, to ensure a higher efficiency of the proposed distributed SPA Gibbs sampler, both in terms of computations and communications, we can identify two conditions on the structure of $\Db$.
These conditions are typically satisfied when $\Db$ is block-sparse, as the operators considered in \Cref{Sec:application}.

The first condition consists in ensuring that, each worker only needs to communicate with a small amount of other workers compared to the total number of workers $K$.

The second condition is that the quantity of variables that two communicating workers need to exchange remains small compared to the variables already stored on each worker.
In particular, for two workers $(k, k')\in \{1, \ldots, K\}^2$, $k\neq k'$, the subpart of $\xb_k$ that need to be sent from $k$ to $k'$ needs to be small compared to $\xb_k$ and $\xb_{k'}$.

\section{Model and hypergraph structure}
\label{Sec:Model}

\input{section3_hypergraph_model_v3}

\section{Distributed SPA Gibbs sampler}
\label{Sec:dist-single-algo}

\input{section5_v1}

\section{Application to supervised image deconvolution}
\label{Sec:application}

\input{section6_v1}

\section{Conclusion}
\label{Sec:conclusion}

In this paper, a distributed block-coordinate SPA Gibbs sampler has been introduced to efficiently solve large scale imaging inverse problems.
The approach leverages the approximate data augmentation scheme AXDA~\cite{Vono_etal_2020,Rendell2021} to efficiently handle composite functions involving linear operators. A block-coordinate approach is adopted to split and distribute all the variables over multiple workers.
The proposed method exploits the hypergraph structure of the linear operators to design a versatile distributed block-coordinate sampler.

Experiments on a supervised image deblurring problem show that the proposed approach forms reliable estimates with quantified uncertainty in a significantly reduced amount of time, compared to a state-of-the-art non-distributed version of the sampler~\cite{Vono2019icassp}. In particular, the proposed sampler is shown to provide estimators at a fraction of the runtime of the serial implementation by using an increasing number of workers $K$. Processing a 1 million pixel image using our current Python implementation takes less than 2 minutes to obtain a good restoration with associated credibility intervals.

Note that the proposed distributed block-coordinate SPA Gibbs sampler is directly applicable to a much wider class of applications than the restoration problems addressed in this work, such as image inpainting, super-resolution or reconstruction. Future works include the development of an asynchronous version of the proposed approach, to further speed up the inference process, while maintaining the convergence of the Markov chain.

\subsection{Hypergraph structure of the model}
\label{Ssec:hypergraph-struct-model}

Let ${\ensuremath{\boldsymbol{\mathsf H}}} = \ensuremath{\mathbb{R}}^{\overline{N}}$ be such that ${\ensuremath{\boldsymbol{\mathsf H}}} = \ensuremath{{\mathsf H}}_1 \times \ldots \times \ensuremath{{\mathsf H}}_N$, where for every $n\in \{1, \ldots, N\}$, $\ensuremath{{\mathsf H}}_n = \ensuremath{\mathbb{R}}^{N_n}$ and $\overline{N} = \sum_{n=1}^N N_n$.
An element of ${\ensuremath{\boldsymbol{\mathsf H}}}$ is denoted by $\xb = (x_n)_{1 \le n \le N}$, where, for every $n\in \{1, \ldots, N\}$, $x_n \in \ensuremath{{\mathsf H}}_n$.
Similarly, let ${\ensuremath{\boldsymbol{\mathsf G}}} = \ensuremath{\mathbb{R}}^{\overline{M}}$ be such that ${\ensuremath{\boldsymbol{\mathsf G}}} = \ensuremath{{\mathsf G}}_1 \times \ldots \times \ensuremath{{\mathsf G}}_M$, where for every $m\in \{1, \ldots, M\}$, $\ensuremath{{\mathsf G}}_m = \ensuremath{\mathbb{R}}^{M_m}$ and $\overline{M} = \sum_{m=1}^M M_m$. Let $\vb = (v_m)_{1 \le m \le M}$ be an element of ${\ensuremath{\boldsymbol{\mathsf G}}}$ such that, for every $m \in \{1, \ldots, M\}$, $v_m \in \ensuremath{{\mathsf G}}_m$.

The distribution of the problem over the different workers will follow the topology of a hypergrah encoded by the structure of $\Db$. To this aim, we assume that $\Db$ holds some block separability structure, and consider the following hypergraph model.
\begin{model}\label{mod:hypergraph}
Let $\Db \colon {\ensuremath{\boldsymbol{\mathsf H}}} \to {\ensuremath{\boldsymbol{\mathsf G}}}$ be such that $\Db = (D_{m,n})_{1 \le m \le M, 1\le n \le N}$, with, for every $(m,n) \in \{1, \ldots,M\} \times \{1, \ldots,N\}$, $D_{m,n} \colon \ensuremath{{\mathsf H}}_n \to \ensuremath{{\mathsf G}}_m$.
Let $\ensuremath{\mathcal{H}} = (\xb, \eb)$ be the hypergraph associated with $\Db$, with $N$ vertices $\xb$, and $M $ hyperedges denoted by $\eb =(e_m)_{1\le m \le M}$, such that
\begin{equation}
    (\forall m \in \{1, \ldots, M\}) \quad
    e_m = \menge{ n \in \{1, \ldots, N\}}{ D_{m,n} \neq 0},
\end{equation}
where $0_{\ensuremath{{\mathsf G}}_m \times \ensuremath{{\mathsf H}}_n}$ is the null element from $\ensuremath{{\mathsf G}}_m$ to $\ensuremath{{\mathsf H}}_n$
\end{model}
The hyperedges of $\ensuremath{\mathcal{H}}$, i.e., the connections between vertices, are described by the topology of $\Db$. Precisely, the $M$ rows of $\Db$ represent the $M$ hyperedges of $\ensuremath{\mathcal{H}}$, and the $N$ columns represent the $N$ vertices of $\ensuremath{\mathcal{H}}$. Hence, as described in \Cref{mod:hypergraph}, for each $m\in \{1, \ldots, M\}$, the hyperedge $e_m$ links nodes $(x_n)$ if $D_{m,n} \neq 0_{\ensuremath{{\mathsf G}}_m \times \ensuremath{{\mathsf H}}_n}$.
Hence, $\Db$ can be seen as a weighted incidence matrix associated with $\ensuremath{\mathcal{H}}$.
In this context, any variable in $ {\ensuremath{\boldsymbol{\mathsf G}}}$ can be seen as a hyperedge weight of the hypergraph $\ensuremath{\mathcal{H}}$.

According to \Cref{mod:hypergraph}, the computation of $\vb = \Db \xb$ in Algorithm~\eqref{algo:SPA-simple-sampling} can be decomposed and computed block-wise as
\begin{equation}
\label{eq:vm-Dmn-xn}
    \vb = (v_m)_{1 \le m \le M} \text{ with }
    (\forall m \in \{1, \ldots, M\}) \quad
    v_m = \sum_{n\in e_m} D_{m,n} x_n.
\end{equation}
In~\eqref{eq:vm-Dmn-xn}, for every $m\in \{1, \ldots, M\}$, only the non-zero blocks of $\Db$ are taken into account, i.e., only using the vertices contained in the hyperedge $e_m$.
In the particular case when all hyperedges are disjoint (i.e., disconnected hypergraph), a very simple distributed algorithm would distribute the computation of the quantities $v_m$ on independent workers. This is not the case in general, and hyperedges $e_m$ can share some vertices $x_n$. A distributed algorithm will thus need to carry out communications between some workers. The next section exploits the hypergraph model described here to distribute the computation of $\Db\xb$, while limiting the communication cost.

\subsection{Distribution of the hypergraph nodes over the workers}
\label{Ssec:Model:dist-hyp-notation}

Let $K \in \ensuremath{\mathbb N}^*$, $K \le N$, be the number of workers available to the user to parallelize the algorithm.
On each worker $k \in \{1, \ldots, K\}$, we store a subpart $\xb_k$ of the nodes $\xb$ of the hypergraph $\ensuremath{\mathcal{H}}$.
The pattern for the distribution of the nodes $\xb$ over the $K$ workers will be driven by the topology of the matrix~$\Db$, as explained in \Cref{Ssec:hypergraph-struct-model}. This pattern needs to be fixed before designing the algorithm.
Note that the blocks $(D_{m,n})_{1 \le m \le M, \, n \in e_m}$ of $\Db$ will be distributed over the workers as well, but will not require to be communicated between the workers.
In particular, we provide below the notation to split the hypergraph (i.e., the operator $\Db$) over the $K$ workers. The notation for the splitting of the nodes and the hyperedge weights are illustrated in \Cref{fig:opsplitting} for a simple example of a block-sparse matrix $\Db$.

\subsubsection{Distributing the hypergraph over the workers}
\label{Sssec:distributed-hypergraph}

Let $(\ensuremath{\mathbb{V}}_k)_{1 \le k \le K}$ be the partition of the set of vertex indices $\{1, \ldots, N\}$ such that, for every $k\in \{1, \ldots, K\}$, $\ensuremath{\mathbb{V}}_k $ is the set of vertex indices handled on worker $k$, and for every $\xb \in {\ensuremath{\boldsymbol{\mathsf H}}}$,  $\xb_k = (x_n)_{n \in \ensuremath{\mathbb{V}}_k}$.
Recall that $k_m \in \{1, \ldots, K\}$ denotes the worker associated with the $m$-th hyperedge $e_m$. The association of hyperedges to workers is a choice left to the user. Then, we must have  $e_m \cap \ensuremath{\mathbb{V}}_{k_m} \neq \ensuremath{{\varnothing}}$, i.e., worker $k_m$ must contain at least one vertex belonging to hyperedge $e_m$.

The set of all the hyperedges associated with worker $k$ is denoted by
\begin{equation}    \label{def:Ebar}
\overline{\ensuremath{\mathbb{E}}}_k = \menge{m \in \{1, \ldots, M\}}{k_m = k}.
\end{equation}
Note that $\big( \overline{\ensuremath{\mathbb{E}}}_k \big)_{1 \leq k \leq K}$ defines a partition of $\{1,\ldots, M\}$ over the $K$ workers.
Note that optimizing the configuration for a specific operator $\Db$ is on its own a resource allocation problem~\cite{Pilla2021} that is out of the scope of the current work.
Using notation~\eqref{def:Ebar}, for any element $\vb \in {\ensuremath{\boldsymbol{\mathsf G}}}$, we can denote by $\vb_k = (v_m)_{m\in \overline{\ensuremath{\mathbb{E}}}_k}$ the hyperegde weights stored on worker $k$.

We also define $\ensuremath{\mathbb{E}}_k$, the set of hyperedge indices only containing vertices stored on worker $k$, as
\begin{equation}    \label{eq:def:Ek}
    \ensuremath{\mathbb{E}}_k = \menge{ m \in \{1, \ldots, M\} }{ n \in e_m \; \Leftrightarrow \; n \in \ensuremath{\mathbb{V}}_k}.
\end{equation}
The sets $(\ensuremath{\mathbb{E}}_k)_{1 \le k \le K}$ identify the rows of $\Db$ whose non-zero elements are multiplied with vertices that are stored on a single worker, i.e., the rows that can be used with no communication between two different workers. In contrast, the computation of $v_m$ for $m\in\overline{\ensuremath{\mathbb{E}}}_k\setminus\ensuremath{\mathbb{E}}_k$ will call for communications.

\subsubsection{Hyperedges overlapping over workers ensuring communications}
\label{Ssec:communication-overlap-hyperedge}

The hyperedges overlaping over multiple workers will require vertices stored on different workers to be communicated between each other.
They correspond to hyperedges $e_m$, with $ m \in \{1, \ldots, M\} \setminus (\bigcup_{k=1}^K \ensuremath{\mathbb{E}}_k)$.

For a fixed hyperedge $e_m$, let $\ensuremath{\mathcal{W}}_m \subset \{1 , \ldots, K\}\setminus \{k_m\}$ be the set of workers other than $k_m$ containing vertices $x_n$ belonging to the same hyperedge $e_m$, i.e.,
\begin{equation*}
\ensuremath{\mathcal{W}}_m = \menge{k \in \{1 , \ldots, K\}\setminus \{k_m\} }{ \exists n \in e_m \text{ such that } n \in \ensuremath{\mathbb{V}}_{k}}.
\end{equation*}
Using this notation, we can give an equivalent definition of \eqref{eq:def:Ek} as, for every $k\in \{1, \ldots, K\}$,
$\ensuremath{\mathbb{E}}_k = \menge{m \in \{1, \ldots, M\}}{k_m = k \text{ and } \ensuremath{\mathcal{W}}_m = \ensuremath{{\varnothing}}}$ (i.e., no overlap, no communication).
For completeness, we also introduce the notation $\overline{\ensuremath{\mathcal{W}}}_m = k_m \cup \ensuremath{\mathcal{W}}_m$.

For every $m \in \{1, \ldots, M\} \setminus (\bigcup_{k=1}^K \ensuremath{\mathbb{E}}_k)$, workers $k' \in \ensuremath{\mathcal{W}}_m$ will need to send some vertex values from $\xb_{k'}$ to the worker $k_m$ so that it can compute $v_m$, as described in \eqref{eq:vm-Dmn-xn}.
For every $k\in \{1, \ldots, K\}$, we denote by $\ensuremath{\mathcal{R}}_k  \subset \{1, \ldots, K\}\setminus \{k\}$ the set of workers $k'$  from which worker $k$ receives vertices. These workers store vertices belonging to a hyperedge of index $m\in \{1, \ldots, M\}$ such that $k_m = k$ but that are not stored on worker $k$.
Similarly, we denote by $\ensuremath{\mathcal{S}}_k \subset \{1, \ldots, K\}\setminus \{k\}$ the set of workers $k'$ to which worker $k$ is sending vertices, i.e., all the workers $k'$ such that $k\in \ensuremath{\mathcal{R}}_{k'}$.

Communications will occur through hyperedges connecting different workers.
For every $k'\in \ensuremath{\mathcal{R}}_k$, let
\begin{equation*}
\ensuremath{\mathbb{E}}_{(k,k')} = \menge{m\in \{1, \ldots, M\}}{k_m = k \text{ and } k' \in \ensuremath{\mathcal{W}}_m },
\end{equation*}
be the set of hyperedges containing the vertices from worker $k'$ that will be needed by worker $k$ to compute $v_m$, so that $\ensuremath{\mathcal{R}}_k = \menge{k' \in \{1, \ldots, K\} \setminus \{k\}}{\exists m \in \ensuremath{\mathbb{E}}_{(k,k')} }$.
Thus, the set of all the hyperedges that will carry out some communication to worker $k$ is
\begin{equation*}
\ensuremath{\mathbb{E}}_{\ensuremath{\mathcal{R}}_k} = \bigcup_{k'\in \ensuremath{\mathcal{R}}_k} \ensuremath{\mathbb{E}}_{(k,k')} = \menge{m \in \{1, \ldots, M\}}{k_m = k \text{ and } \ensuremath{\mathcal{W}}_m \neq \ensuremath{{\varnothing}}}.
\end{equation*}
As a result, $\ensuremath{\mathbb{E}}_k$ and $\ensuremath{\mathbb{E}}_{\ensuremath{\mathcal{R}}_k}$ form a partition of $\overline{\ensuremath{\mathbb{E}}}_k$,
where $\ensuremath{\mathbb{E}}_k$ corresponds to the set of weights $v_m$ that can be computed locally on worker $k$ (i.e. hyperedges inducing no communication), and $\ensuremath{\mathbb{E}}_{\ensuremath{\mathcal{R}}_k}$ corresponds to the weights $v_m$ that necessitate vertices provided by other workers, and therefore communications.

Each hyperedge $e_m$ potentially sends vertices towards the corresponding worker $k_m$. To identify vertices that are communicated between workers, let
\begin{equation*}
\ensuremath{\mathbb{V}}_{(k,k')} = \menge{n \in \ensuremath{\mathbb{V}}_{k'}}{\exists m \in \ensuremath{\mathbb{E}}_{(k,k')} \text{ such that } n \in e_m}
\end{equation*}
be the indices of vertices that are received by worker $k$ from $k'$.
The set of vertices sent to worker $k_m$ by all workers in $\ensuremath{\mathcal{W}}_m$ is denoted by $\ensuremath{\mathbb{V}}_{\ensuremath{\mathcal{W}}_m} = \bigcup_{k' \in \ensuremath{\mathcal{W}}_m} \ensuremath{\mathbb{V}}_{(k_m, k')}$. For completeness, we also introduce $\ensuremath{\mathbb{V}}_{\overline{\ensuremath{\mathcal{W}}}_m} = \ensuremath{\mathbb{V}}_{k_m} \cup \ensuremath{\mathbb{V}}_{\ensuremath{\mathcal{W}}_m}$, the set of all the vertex indices necessary to perform computations associated with the $m$-th hyperedge.

For every $\xb \in {\ensuremath{\boldsymbol{\mathsf H}}}$, we denote by $\xb_{(k,k')} = (x_n)_{n \in \ensuremath{\mathbb{V}}_{(k,k')}}$ the vertices communicated from worker $k'$ to worker $k$ through some hyperedge $e_m$ such that $k_m=k$.
For every $m\in \{1, \ldots,M\} \setminus (\bigcup_{k=1}^K \ensuremath{\mathbb{E}}_k)$, let $\xb_{\overline{\ensuremath{\mathcal{W}}}_m} = (x_n)_{n \in \ensuremath{\mathbb{V}}_{\overline{\ensuremath{\mathcal{W}}}_m}}$ be the concatenation of vertices stored on worker $k_m$ and vertices sent from all workers $k'\in \ensuremath{\mathcal{W}}_m$ to worker $k_m$ (i.e., all the vertices such that $D_{m,n} \neq 0$).

\begin{table}[t!]
    \centering\footnotesize
    \begin{tabular}{@{}p{4cm}@{}|p{8.1cm}|@{}c@{}|@{}c@{}} \toprule
        \multirow{2}{*}{Notation} & \multirow{2}{*}{Definition} & \multicolumn{2}{c}{Given by}    \\
        && operator $\Db$ & $\phantom{a}$ user $\phantom{a}$ \\
        \hline
        $n\in \{1, \ldots, N\}$ & Indices for vertices & \checkmark &  \\
        $m\in \{1, \ldots, M\}$ & Indices for hyperedges & \checkmark & \\
        $\eb = (e_m)_{1 \le m \le M}$ & Hyperedges & \checkmark & \\
        $e_m \subset \{1, \ldots, N\}$ & Vertex indices in hyperedge $m$ & \checkmark & \\
        $k \in \{1, \ldots, K\}$ & Indices for workers & & \checkmark \\
        $\ensuremath{\mathcal{R}}_k \subset \{1, \ldots, K\}\setminus\{k\}$ & Indices of workers $k' \in \{1, \ldots, K\}\setminus\{k\}$ from which worker~$k$ receives vertex values & & \checkmark \\
        $\ensuremath{\mathcal{S}}_k \subset \{1, \ldots, K\}\setminus\{k\}$ & Indices of workers $k' \in \{1, \ldots, K\}\setminus\{k\}$  to which worker~$k$ sends vertex values & & \checkmark \\
        $\ensuremath{\mathbb{V}}_k \subset \{1, \ldots, N\}$ & Indices of vertices stored on worker~$k$, such that $(\ensuremath{\mathbb{V}}_k)_{1 \le k \le K}$ is a partition of $\{1, \ldots, N\}$ & \checkmark & \checkmark \\
        $k_m \in\{1, \ldots, K\}$ & Worker associated with $m$-th hyperedge $e_m$ (chosen by the user). $k_m$ must satisfy $e_m \cap \ensuremath{\mathbb{V}}_{k_m} \neq \ensuremath{{\varnothing}}$  & & \checkmark \\
        $\ensuremath{\mathcal{W}}_m \subset \{1, \ldots, K\} \setminus \{k_m\}$ & Set of all workers but $k_m$, containing vertices from $e_m$ & \checkmark & \\
        $\overline{\ensuremath{\mathcal{W}}}_m \subset \{1, \ldots, K\} $ & $\overline{\ensuremath{\mathcal{W}}}_m = k_m \cup \ensuremath{\mathcal{W}}_m$ & \checkmark & \\
        $\ensuremath{\mathbb{V}}_{(k,k')} \subset \ensuremath{\mathbb{V}}_{k'}$ & Indices of vertices sent from worker $k'$ to worker~$k$ & \checkmark & \checkmark    \\
        $\ensuremath{\mathbb{V}}_{\ensuremath{\mathcal{W}}_m} \subset \{1, \ldots, N\}\setminus \ensuremath{\mathbb{V}}_{k_m}$ & $\ensuremath{\mathbb{V}}_{\ensuremath{\mathcal{W}}_m} = \bigcup_{k' \in \ensuremath{\mathcal{W}}_m} \ensuremath{\mathbb{V}}_{(k_m, k')}$ the set of vertex indices that will be communicated to worker $k_m$ from all workers $k'\in \ensuremath{\mathcal{W}}_m$ & \checkmark & \checkmark    \\
        $\ensuremath{\mathbb{V}}_{\overline{\ensuremath{\mathcal{W}}}_m} \subset \{1, \ldots, N\}$ & $\ensuremath{\mathbb{V}}_{\overline{\ensuremath{\mathcal{W}}}_m} = \ensuremath{\mathbb{V}}_{k_m} \cup \ensuremath{\mathbb{V}}_{\ensuremath{\mathcal{W}}_m}$ the set of vertex indices necessary to perform computations associated with $k_m$ & \checkmark & \checkmark    \\
        $\ensuremath{\mathbb{E}}_k \subset \{1, \ldots, M\}$ & Indices of hyperedges only containing vertices stored on worker~$k$ & \checkmark & \checkmark \\
        $\ensuremath{\mathbb{E}}_{(k,k')} \subset \{1, \ldots, M\}$ & Indices of hyperedges containing vertices sent from worker $k'$ to worker~$k$ & \checkmark & \checkmark \\
        $\ensuremath{\mathbb{E}}_{\ensuremath{\mathcal{R}}_k} \subset \{1, \ldots, M\}$ & $\ensuremath{\mathbb{E}}_{\ensuremath{\mathcal{R}}_k} = \cup_{k'\in \ensuremath{\mathcal{R}}_k} \ensuremath{\mathbb{E}}_{(k,k')}$ set of all hyperedges containing vertices that will be communicated to worker $k$ & \checkmark & \checkmark    \\
        $\overline{\ensuremath{\mathbb{E}}}_k \subset \{1, \ldots, M\}$ & $ \overline{\ensuremath{\mathbb{E}}}_k = \ensuremath{\mathbb{E}}_k \cup \ensuremath{\mathbb{E}}_{\ensuremath{\mathcal{R}}_k} $, such that $(\overline{\ensuremath{\mathbb{E}}}_k)_{1 \le k \le K}$ is a partition of $\{1, \ldots, M\}$ & \checkmark & \checkmark \\
        \bottomrule
    \end{tabular}
    \caption{\label{table:notation-hypergraph}\small
    Sets used to define the hypergraph structure of the problem and the distributed algorithm.
    }
\end{table}

\begin{remark}
Let $k \in \{1, \ldots, K\}$.
For every $k'\in {\ensuremath{\mathcal{S}}}_k$, since $\ensuremath{\mathbb{V}}_{(k,k')}$ denotes the set of vertices that need to be communicated from worker $k'$ to worker $k$, the (reciprocal) set $\ensuremath{\mathbb{V}}_{(k',k)} \subset \ensuremath{\mathbb{V}}_k$ corresponds to the set of vertices that will be sent back from worker $k$ to $k'$ in the distributed implementation of the adjoint operator $\Db^*$.
Similarly, with $\ensuremath{\mathbb{E}}_{(k,k')}$ the set of hyperedges containing vertices that will be communicated from worker $k'$ to worker $k$, the set $\ensuremath{\mathbb{E}}_{(k',k)} \subset \{1, \ldots, M\}$ corresponds to the hyperedge indices containing vertices that will be sent back from worker $k$ to worker $k'$.
\end{remark}

\subsubsection{Splitting of $\Db$ over the workers} \label{Sssec:splitting_D}

Eventually, we will introduce some notation to split the matrix $\Db$ over the workers.
In practice, the input space ${\ensuremath{\boldsymbol{\mathsf G}}}$ and the output space ${\ensuremath{\boldsymbol{\mathsf H}}}$ are partitioned using the sets $(\overline{\ensuremath{\mathbb{E}}}_k)_{k\in \{1, \ldots, K\}}$ and $(\ensuremath{\mathbb{V}}_k)_{k\in \{1, \ldots, K\}}$, respectively. For every $k \in \{1, \ldots, K\}$, let ${\ensuremath{\boldsymbol{\mathsf G}}}_k = \ensuremath{\raisebox{-0.5mm}{\mbox{\Large{$\times$}}}}_{m \in \overline{\ensuremath{\mathbb{E}}}_k} \ensuremath{{\mathsf G}}_m$ and ${\ensuremath{\boldsymbol{\mathsf H}}}_k=\ensuremath{\raisebox{-0.5mm}{\mbox{\Large{$\times$}}}}_{n \in \ensuremath{\mathbb{V}}_k} \ensuremath{{\mathsf H}}_n$ such that ${\ensuremath{\boldsymbol{\mathsf G}}} = \ensuremath{\raisebox{-0.5mm}{\mbox{\Large{$\times$}}}}_{1\leq k \leq K} {\ensuremath{\boldsymbol{\mathsf G}}}_k$ and ${\ensuremath{\boldsymbol{\mathsf H}}} = \ensuremath{\raisebox{-0.5mm}{\mbox{\Large{$\times$}}}}_{1\leq k \leq K} {\ensuremath{\boldsymbol{\mathsf H}}}_k$.

For every $k\in \{1, \ldots, K\}$, the subparts $\vb_k$ of $\Db\xb$ that belong to $ {\ensuremath{\boldsymbol{\mathsf G}}}_k$ are stored on worker~$k$. The subparts of $\Db$ acting on ${\ensuremath{\boldsymbol{\mathsf H}}}_k$ only are associated with vertices in $\ensuremath{\mathbb{V}}_k$. For every $m \in \ensuremath{\mathbb{E}}_k$, they will be denoted by $\Db_{m,k} = (D_{m,n})_{n \in \ensuremath{\mathbb{V}}_k}$. These subparts of $\Db$ are involved in purely local computations only.
Similarly, for every $k' \in \ensuremath{\mathcal{R}}_k$, for every $m \in \ensuremath{\mathbb{E}}_{(k,k')}$, the subparts of $\Db$ acting on vertices $\ensuremath{\mathbb{V}}_{(k,k')}$ will be denoted by $\Db_{m,(k,k')} = (D_{m,n})_{n \in \ensuremath{\mathbb{V}}_k \cup \ensuremath{\mathbb{V}}_{(k,k')}} $.
In addition, for every $m \in \ensuremath{\mathbb{E}}_{\ensuremath{\mathcal{R}}_k}$, we will denote by $\Db_{m,\overline{\ensuremath{\mathcal{W}}}_m} = (D_{m,n})_{n \in \ensuremath{\mathbb{V}}_{\overline{\ensuremath{\mathcal{W}}}_m}}$ the subparts of $\Db$ acting on vertices stored either on worker $k_m$ or sent to worker $k_m$ by all other workers $k'\in \ensuremath{\mathcal{W}}_m$. Hence
\begin{equation}
    \label{eq:vm-formula}
    v_m =
    \begin{cases}
    \Db_{m,k} \xb_k, &  \text{if } m \in \ensuremath{\mathbb{E}}_k, \\
    \Db_{m,\overline{\ensuremath{\mathcal{W}}}_m} \xb_{\overline{\ensuremath{\mathcal{W}}}_m}, &   \text{otherwise}.
    \end{cases}
\end{equation}
As a result, the subparts of $\Db$ which must be stored on worker $k$ are the $\Db_{m,k}$, for every $m\in\ensuremath{\mathbb{E}}_k$, and the $\Db_{m,\overline{\ensuremath{\mathcal{W}}}_m}$, for every $m \in \ensuremath{\mathbb{E}}_{\ensuremath{\mathcal{R}}_k}$.

The notation given above are summarized in~\Cref{table:notation-hypergraph} and \Cref{table:notation2-hypergraph}. A simple example for block-sparse matrices is also provided below, and illustrated in~\Cref{fig:opsplitting}. This example will be reused later on to illustrate the communications between the workers (see~\Cref{fig:comp-comm-D}).

\begin{figure}
\begin{center}
\includegraphics[width=0.5\textwidth,keepaspectratio]{figures/illustrations/operator_notation_matrix_comp.png}
\end{center}

\caption{\label{fig:opsplitting}\small
Distribution of vertices (horizontal lines) and hyperedges (vertical lines) of the hypergraph over $K=3$ workers, for a block-sparse matrix $\Db$ corresponding to~\Cref{ex:block-sparse}.
}
\end{figure}

\begin{example}\label{ex:block-sparse}
To illustrate some of the notation introduced above, consider the block-sparse matrix $\Db$ shown in~\Cref{fig:opsplitting}, with $K=3$ workers.
For the vertices, $\ensuremath{\mathbb{V}}_1$ (resp. $\ensuremath{\mathbb{V}}_2$ and $\ensuremath{\mathbb{V}}_3$) contains the vertices handled on worker $k=1$ (resp. $k=2$ and $k=3$).
In addition, $\ensuremath{\mathbb{V}}_{(1, 2)}$ (resp. $\ensuremath{\mathbb{V}}_{(2, 3)}$) identifies the vertices of $\xb_2$ (resp. $\xb_3$) that will need to be communicated to worker $k=1$ (resp. $k=2$).
For the hyperedges, $\overline{\ensuremath{\mathbb{E}}}_1$ (resp. $\overline{\ensuremath{\mathbb{E}}}_2$ and $\overline{\ensuremath{\mathbb{E}}}_3$) identifies the hyperedge weights stored on worker $k=1$ (resp. $k=2$ and $k=3$).
The set $\ensuremath{\mathbb{E}}_1$ (resp. $\ensuremath{\mathbb{E}}_2$ and $\ensuremath{\mathbb{E}}_3$) identifies the hyperedges fully contained in worker $k=1$ (resp. $k=2$ and $k=3$).
Hyperedges in $\ensuremath{\mathbb{E}}_{(1, 2)}$ (resp. $\ensuremath{\mathbb{E}}_{(2, 3)}$) contain the vertices that will be communicated from worker $k=2$ to worker $k=1$ (resp. $k=2$ to $k=3$).
In this example, only the blocks $\Big( (\Db_{m,k})_{m \in \ensuremath{\mathbb{E}}_k} \Big)_{1 \le k \le 3}$, $(\Db_{m, (1,2)})_{m \in \ensuremath{\mathbb{E}}_{(1,2)}}$ and $(\Db_{m, (2,3)})_{m \in \ensuremath{\mathbb{E}}_{(2,3)}}$ are non-zero (e.g., convolution operator).
Blocks $\Big( (\Db_{m,1})_{m\in \ensuremath{\mathbb{E}}_1}, (\Db_{m, (1,2)})_{m\in \ensuremath{\mathbb{E}}_{(1,2)}} \Big)$
(resp. $\Big( (\Db_{m,2})_{m \in \ensuremath{\mathbb{E}}_2}, (\Db_{m, (2,3)})_{m \in \ensuremath{\mathbb{E}}_{}(2, 3)} \Big)$
and $(\Db_{m,3})_{m \in \ensuremath{\mathbb{E}}_3}$) are stored on worker $k=1$ (resp. $k=2$ and $k=3$).
For worker $k=1$ (resp. $k=2$ and $k=3$), workers sending vertex values to this worker is given by $\ensuremath{\mathcal{R}}_1 = \{2\}$ (resp. $\ensuremath{\mathcal{R}}_2 = \{3\}$ and $\ensuremath{\mathcal{R}}_3 = \ensuremath{{\varnothing}}$). Similarly, workers receiving vertex values from this worker is given by $\ensuremath{\mathcal{S}}_1 = \ensuremath{{\varnothing}}$ (resp. $\ensuremath{\mathcal{S}}_2 = \{1\}$ and $\ensuremath{\mathcal{S}}_3 = \{2\})$.
For every $m\in \overline{\ensuremath{\mathbb{E}}}_1$ (resp. $\overline{\ensuremath{\mathbb{E}}}_2$ and $\overline{\ensuremath{\mathbb{E}}}_3$), $k_m = 1$ (resp. $k_m = 2$ and $k_m = 3$).
For every $m \in \ensuremath{\mathbb{E}}_{(1, 2)}$, $\ensuremath{\mathcal{W}}_m = \ensuremath{\mathcal{R}}_1 = \ensuremath{\mathcal{S}}_2 = \{2\}$, and for every $m\in \ensuremath{\mathbb{E}}_{(2,3)}$, $\ensuremath{\mathcal{W}}_m = \ensuremath{\mathcal{R}}_2 = \{3\}$.
Finally, for every $m\in \ensuremath{\mathbb{E}}_{\ensuremath{\mathcal{R}}_1}= \ensuremath{\mathbb{E}}_{(1,2)}$ (resp. $ \ensuremath{\mathbb{E}}_{\ensuremath{\mathcal{R}}_2}= \ensuremath{\mathbb{E}}_{(2,3)}$), $\Db_{m, \overline{\ensuremath{\mathcal{W}}}_m} = (D_{m,n})_{n \in \ensuremath{\mathbb{V}}_1 \cup \ensuremath{\mathbb{V}}_{(1,2)}}$ (resp.  $\Db_{m, \overline{\ensuremath{\mathcal{W}}}_m} = (D_{m,n})_{n \in \ensuremath{\mathbb{V}}_2 \cup \ensuremath{\mathbb{V}}_{(2,3)}}$). These last notation correspond to the orange rectangles displayed in~\Cref{fig:blocksparse}. In this example, $\ensuremath{\mathbb{E}}_{\ensuremath{\mathcal{R}}_3} = \ensuremath{{\varnothing}}$.
\end{example}

\subsection{Separability assumptions}
\label{Ssec:Model:separability}

We can now give assumptions that will be used to design an efficient distributed SPA Gibbs sampler.
Precisely, we will consider distribution~\eqref{eq:dist-SPA}, where the functions are assumed to satisfy some separability conditions.
\begin{assumption}  \label{ass:gen}\
\begin{enumerate}
    \item\label{ass:gen:i}
    $h \colon {\ensuremath{\boldsymbol{\mathsf H}}} \to \ensuremath{\mathbb{R}}$ is $\lambda_h$-Lipschitz differentiable, with $\lambda_h>0$.

    \item \label{ass:gen:ii}
    $h \colon {\ensuremath{\boldsymbol{\mathsf H}}} \to \ensuremath{\mathbb{R}}$ and $f  \colon {\ensuremath{\boldsymbol{\mathsf H}}} \to \ensuremath{\left]-\infty,+\infty\right]}$ are additively separable on the workers, i.e.
    \begin{equation}
        (\forall \xb \in {\ensuremath{\boldsymbol{\mathsf H}}}) \quad
            h(\xb) = \sum_{k=1}^K h_k(\xb_k),
            \text{ and }
            f(\xb) = \sum_{k=1}^K f_k(\xb_k),
    \end{equation}
    where, for every $k \in \{1, \ldots, K\}$, $h_k \colon {\ensuremath{\boldsymbol{\mathsf H}}}_k \to \ensuremath{\mathbb{R}}$ and $f_k \colon {\ensuremath{\boldsymbol{\mathsf H}}}_k \to \ensuremath{\left]-\infty,+\infty\right]}$.

    \item\label{ass:gen:iii}
    $g \colon {\ensuremath{\boldsymbol{\mathsf G}}} \to \ensuremath{\left]-\infty,+\infty\right]}$ is additively separable, i.e.
    \begin{equation}
        (\forall \vb  \in {\ensuremath{\boldsymbol{\mathsf G}}}) \quad
            g(\vb) = \sum_{m=1}^M g_m(v_m), \label{eq:g-add}
    \end{equation}
    with, for every $m\in \{1, \ldots,M\}$, $g_m \colon \ensuremath{{\mathsf G}}_m \to \ensuremath{\left]-\infty,+\infty\right]}$.

    \item\label{ass:gen:iv}
    $\phi_\alpha$ and $\psi_\beta$ are given by
    \begin{align}
    (\forall (\vb, \ub) \in {\ensuremath{\boldsymbol{\mathsf G}}}^2)\quad
    &  \phi_\alpha(\vb, \ub) = \sum_{m=1}^M \phi_{m, \alpha} (v_m, u_m) = \sum_{m=1}^M  \frac{1}{2\alpha^2} \| v_m - u_m \|^2, \label{eq:phi-add} \\
    &   \psi_\beta(\ub) = \sum_{m=1}^M \psi_{m, \beta} (u_m) = \sum_{m=1}^M  \frac{1}{2\beta^2} \| u_m \|^2 .  \label{eq:psi-add}
    \end{align}
\end{enumerate}
\end{assumption}

The above assumptions can be rewritten to highlight the separability over the workers. Indeed, there exists a permutation $\varrho \colon {\ensuremath{\boldsymbol{\mathsf G}}} \to {\ensuremath{\boldsymbol{\mathsf G}}}$ such that
\begin{align}
    (\forall \xb \in {\ensuremath{\boldsymbol{\mathsf H}}}) \quad
    \Db \xb
    = \left( \left( D_{m,n} \right)_{1 \le n \le N} \xb \right)_{1 \le m \le M}
    &= \varrho \left(\left( \begin{matrix}
    \left(  \Db_{m,k}\right)_{m \in \ensuremath{\mathbb E}_k}  \xb_k     \\[0.1cm]
    \left(  \Db_{m,\overline{\ensuremath{\mathcal{W}}}_m}  \xb_{\overline{\ensuremath{\mathcal{W}}}_m} \right)_{m \in \ensuremath{\mathbb{E}}_{\ensuremath{\mathcal{R}}_k}}
    \end{matrix} \right)_{1 \le k \le K} \right). \label{eq:def:perm}
\end{align}
Reorganizing the vertex indices emphasizes the distinction between local computations and those for which communications are required.
For every $(\ub, \vb) \in {\ensuremath{\boldsymbol{\mathsf G}}}^2$ and $\xb \in {\ensuremath{\boldsymbol{\mathsf H}}}$, \eqref{eq:g-add}, \eqref{eq:phi-add} and \eqref{eq:psi-add} can thus be rewritten as
\begin{align}
    g(\vb)
    &=  \sum_{k=1}^K \Big(
        \sum_{m \in \overline{\ensuremath{\mathbb{E}}}_k} g_m \left( v_m  \right)
        \Big),    \label{eq:g-dist-simple}   \\
    \phi_\alpha(\Db \xb, \ub)
    &=  \sum_{k=1}^K \Big(
        \sum_{m \in \overline{\ensuremath{\mathbb{E}}}_k} \phi_{m, \alpha} \left( [\Db \xb]_m, u_m  \right)
        \Big),  \\
    &=  \sum_{k=1}^K \Big(
        \sum_{m \in \ensuremath{\mathbb{E}}_k} \phi_{m,\alpha} \left( \Db_{m,k} \xb_k, u_m \right)
        + \sum_{m \in \ensuremath{\mathbb{E}}_{\ensuremath{\mathcal{R}}_k}}  \phi_{m,\alpha} \left( \Db_{m,\overline{\ensuremath{\mathcal{W}}}_m} \xb_{\overline{\ensuremath{\mathcal{W}}}_m}, u_m \right)
        \Big),   \label{eq:phi-dist-simple}    \\
    \psi_\beta(\ub)
    &=  \sum_{k=1}^K \Big(
        \sum_{m \in \overline{\ensuremath{\mathbb{E}}}_k} \psi_{m,\beta} \left( u_m \right)
        \Big),    \label{eq:psi-dist-simple}
\end{align}
where $[\cdot]_m$ denotes the $m$-th element of its argument.

\begin{table}[t!]
    \centering\footnotesize
    \begin{tabular}{@{}p{4.5cm}@{}|@{}p{11.1cm}@{}} \toprule
        Notation & Definition  \\
        \hline
        $\xb = (x_n)_{1 \le n \le N} \in {\ensuremath{\boldsymbol{\mathsf H}}}$ & Vertex values \\
        $\Db \colon {\ensuremath{\boldsymbol{\mathsf H}}} \to {\ensuremath{\boldsymbol{\mathsf G}}}$ & $\Db = (D_{m,n})_{1 \le m \le M, 1\le n \le N}$ linear operator defining the hypergraph structure \\
        %
        $\ub$, $\vb\in {\ensuremath{\boldsymbol{\mathsf G}}}$ & \multirow{2}{*}{Hyperedge weights} \\
        $\vb = (v_m)_{1 \le m \le M} \in {\ensuremath{\boldsymbol{\mathsf G}}}$ & \\
        %
        $\vb_k \in {\ensuremath{\boldsymbol{\mathsf G}}}_k$ & $\vb_k = (v_m)_{m \in \overline{\ensuremath{\mathbb{E}}}_k}$ hyperedge weights stored on worker $k$ \\
        %
        $\xb_k  \in {\ensuremath{\boldsymbol{\mathsf H}}}_k$ & $\xb_k = (x_n)_{n \in \ensuremath{\mathbb{V}}_k}$ vertex values stored on worker $k$ \\
        %
        $\Db_{m,k} \colon {\ensuremath{\boldsymbol{\mathsf H}}}_k \to
        \ensuremath{{\mathsf G}}_m $ & $\Db_{m,k} = (D_{m,n})_{n \in \ensuremath{\mathbb{V}}_k}$, for $m\in \ensuremath{\mathbb{E}}_k$, subpart of $\Db$ stored on worker $k$, associated with hyperedges containing vertices only on worker $k$ \\
        %
        $\xb_{(k,k')} \in \underset{n\in \ensuremath{\mathbb{V}}_{(k,k')}}{\ensuremath{\raisebox{-0.5mm}{\mbox{\LARGE{$\times$}}}}} \ensuremath{{\mathsf H}}_n$  & $\xb_{(k,k')} = (x_n)_{n \in \ensuremath{\mathbb{V}}_{(k,k')}}$ vertex values sent from worker $k'$ to worker $k$ \\
        %
        $\Db_{m,(k,k')} \colon \underset{n\in \ensuremath{\mathbb{V}}_{(k,k')}}{\ensuremath{\raisebox{-0.5mm}{\mbox{\LARGE{$\times$}}}}} \ensuremath{{\mathsf H}}_n \to
        \ensuremath{{\mathsf G}}_m $ & $\Db_{m,(k,k')} = (D_{m,n})_{n \in \ensuremath{\mathbb{V}}_{(k,k')}}$, for $m\in \ensuremath{\mathbb{E}}_{(k,k')}$, subpart of $\Db$ stored on worker $k$, associated with hyperedges containing vertices overlapping workers $k$ and $k'$ \\
        %
        $\xb_{\overline{\ensuremath{\mathcal{W}}}_m} \in \underset{n\in \ensuremath{\mathbb{V}}_{\overline{\ensuremath{\mathcal{W}}}_m}}{\ensuremath{\raisebox{-0.5mm}{\mbox{\LARGE{$\times$}}}}} \ensuremath{{\mathsf H}}_n$  & $\xb_{\overline{\ensuremath{\mathcal{W}}}_m} = (x_n)_{n \in \ensuremath{\mathbb{V}}_{\overline{\ensuremath{\mathcal{W}}}_m}}$, concatenation of vertex values stored on worker $k_m$, and those sent from all worker $k' \in \ensuremath{\mathcal{W}}_m$ to worker $k_m$ \\
        %
        $\Db_{m,\overline{\ensuremath{\mathcal{W}}}_m} \colon \underset{n\in \ensuremath{\mathbb{V}}_{\overline{\ensuremath{\mathcal{W}}}_m} }{\ensuremath{\raisebox{-0.5mm}{\mbox{\LARGE{$\times$}}}}} \ensuremath{{\mathsf H}}_n \to
        \ensuremath{{\mathsf G}}_m $ &
        $\Db_{m,\overline{\ensuremath{\mathcal{W}}}_m} = (D_{m,n})_{n \in \ensuremath{\mathbb{V}}_{\overline{\ensuremath{\mathcal{W}}}_m}}$,
        for $m\in \ensuremath{\mathbb{E}}_{\ensuremath{\mathcal{R}}_k}$, subpart of $(D_{m,n})_{1 \le n \le N}$ corresponding to vertices stored either on worker $k_m$ or on other workers of $\ensuremath{\mathcal{W}}_m$ \\
        %
        \bottomrule
    \end{tabular}
    \caption{\label{table:notation2-hypergraph}\small
    Notation used for the variables involved in the distributed algorithm.
    }
\end{table}

\subsection{Proposed distributed algorithm}
\label{Ssec:algo}

Using the model from the previous section, we are now able to introduce a distributed block-coordinate version of Algorithm~\eqref{algo:SPA-simple-sampling}.
The term \emph{block-coordinate} is to be interpreted as in the optimization literature, where the variables from \eqref{eq:dist-SPA} are divided into $K$ blocks distributed over the $K$ workers, and updated in parallel.
The proposed distributed algorithm is the following.

\begin{proposition} \label{prop:synch-SPA-simple}
Consider a distribution~\eqref{eq:dist-SPA} satisfying \Cref{ass:gen}. Assume that the operator $\Db$ follows \Cref{mod:hypergraph}, and is split over workers $\{1, \ldots, K\}$ such that, for every $k\in \{1, \ldots, K\}$, $(\Db_{m,k})_{m \in \overline{\ensuremath{\mathbb{E}}}_k}$ is stored on worker $k$.
For every $k\in \{1, \ldots, K\}$, let $\xb_k^{(0)} \in {\ensuremath{\boldsymbol{\mathsf H}}}_k$, $ \zb_k^{(0)} \in {\ensuremath{\boldsymbol{\mathsf G}}}_k$, and $ \ub_k^{(0)} \in {\ensuremath{\boldsymbol{\mathsf G}}}_k$.
Let $(\xb^{(t)}, \zb^{(t)}, \ub^{(t)})_{1 \le t \le T}$ be samples generated by the following algorithm:
\begin{equation}    \label{algo:single-synch}
\begin{array}{l}
\text{for } k = 1, \ldots, K \\
\left\lfloor
\begin{array}{l}
    %
    \text{for } k' \in \ensuremath{\mathcal{S}}_k \\[0.1cm]
    \left\lfloor
    \begin{array}{l}
        \text{Send } (x_n^{(0)})_{n \in \ensuremath{\mathbb{V}}_{(k',k)}} \text{ to worker $k'$},  \\
    \end{array}
    \right. \\[0.1cm]
    \text{for } k' \in \ensuremath{\mathcal{R}}_k \\
    \left\lfloor
    \begin{array}{l}
        \text{Receive } (x_n^{(0)})_{n \in \ensuremath{\mathbb{V}}_{(k,k')}} \text{ from worker $k'$},  \\
    \end{array}
    \right. \\[0.3cm]
    \displaystyle \vb_k^{(0)} =
        \left(\begin{matrix}
        ( \Db_{m,k} \xb_k^{(0)} )_{m \in \ensuremath{\mathbb{E}}_k}  \\
        ( \Db_{m,\overline{\ensuremath{\mathcal{W}}}_m} \xb_{\overline{\ensuremath{\mathcal{W}}}_m}^{(0)} )_{m \in \ensuremath{\mathbb{E}}_{\ensuremath{\mathcal{R}}_k}}
        \end{matrix}\right), \\[0.2cm]
\end{array}
\right. \\[0.2cm]
\text{for } t = 0,1, \ldots, T	\\
\left\lfloor
\begin{array}{l}
    \text{for } k = 1, \ldots, K	\\
    \left\lfloor
    \begin{array}{l}
        \displaystyle

        (d_m^{(t)})_{m \in \overline{\ensuremath{\mathbb{E}}}_k}
            = \left( \phi_{m,\alpha}' (\cdot, z_m^{(t)}-u_m^{(t)}) (v_m^{(t)}) \right)_{m \in \overline{\ensuremath{\mathbb{E}}}_k}, \\[0.2cm]
        \text{for } k' \in \ensuremath{\mathcal{R}}_k \\
        \left\lfloor
        \begin{array}{l}
            \displaystyle \widetilde{\db}_{(k',k)}^{(t)} = \sum_{m \in \ensuremath{\mathbb{E}}_{(k',k)}} \Db_{m,k'}^* d_m^{(t)},
            \text{Send } \widetilde{\db}_{(k',k)}^{(t)} \text{ to worker $k'$}, \\[0.2cm]
        \end{array}
        \right.\\[0.2cm]
        \text{for } k' \in \ensuremath{\mathcal{S}}_k \\
        \left\lfloor
        \begin{array}{l}
            \text{Receive } \widetilde{\db}_{(k,k')}^{(t)} \text{ from worker $k'$}, \\[0.2cm]
        \end{array}
        \right.\\
        \displaystyle \ensuremath{\boldsymbol{\delta}}_k^{(t)} = \sum_{m \in \overline{\ensuremath{\mathbb{E}}}_k} \Db_{m,k}^* d_m^{(t)} + \sum_{k'\in \ensuremath{\mathcal{S}}_k} \widetilde{\db}_{(k,k')}^{(t)},     \\[0.1cm]
        \displaystyle \xb_k^{(t+1)} = \ensuremath{\text{\rm prox}}_{\gamma f_k}\left( \xb_k^{(t)} - \gamma \nabla h_k(\xb_k^{(t)}) - \gamma \ensuremath{\boldsymbol{\delta}}_k^{(t)} + \sqrt{2\gamma} \wb_k^{(t)} \right), \\[0.2cm]
        %
        \text{for } k' \in \ensuremath{\mathcal{S}}_k \\
        \left\lfloor
        \begin{array}{l}
            \text{Send } (x_n^{(t+1)})_{n \in \ensuremath{\mathbb{V}}_{(k',k)}} \text{ to worker $k'$},  \\
        \end{array}
        \right.\\[0.1cm]
        \text{for } k' \in \ensuremath{\mathcal{R}}_k \\
        \left\lfloor
        \begin{array}{l}
            \text{Receive } (x_n^{(t+1)})_{n \in \ensuremath{\mathbb{V}}_{(k,k')}} \text{ from worker $k'$},   \\[0.2cm]
        \end{array}
        \right. \\[0.2cm]
        %

        \displaystyle \vb_k^{(t+1)} = \rho
        \left(\begin{matrix}
        ( \Db_{m,k} \xb_k^{(t+1)} )_{m \in \ensuremath{\mathbb{E}}_k}  \\
        ( \Db_{m,\overline{\ensuremath{\mathcal{W}}}_m} \xb_{\overline{\ensuremath{\mathcal{W}}}_m}^{(t+1)} )_{m \in \ensuremath{\mathbb{E}}_{\ensuremath{\mathcal{R}}_k}}
        \end{matrix}\right), \\[0.3cm]
        %
        \zb_k^{(t+1)} \sim \trans{\zb_k,\alpha}{ \vb_k^{(t+1)}, \zb_k^{(t)}, \ub_k^{(t)} }, \\[0.2cm]
        %
        \ub_{k}^{(t+1)} \sim \ensuremath{\mathcal{N}} \Big( \frac{\beta^2}{\alpha^2+\beta^2} (\zb_{k}^{(t+1)} - \vb_{k}^{(t+1)}), \frac{\alpha^2+\beta^2}{\alpha^2 \beta^2} \ensuremath{\mathbf{I}}\, \Big),
    \end{array}
\right.
\end{array}
\right.
\end{array}
\end{equation}
where $\gamma \in ]0, (\lambda_h + \| \Db \|^2/\alpha^2)^{-1}[$ and, for every $k\in \{1, \ldots, K\}$, $(\wb^{(t)}_k)_{1 \le t \le T}$ is a sequence of i.i.d. standard Gaussian random variables in ${\ensuremath{\boldsymbol{\mathsf H}}}_k$.
In addition, the transition kernel
$\ensuremath{\mathcal{K}}_{\zb_k, \alpha}$ has the following invariant distribution
\begin{align}
       \distc{k, \alpha}{\zb_k}{\vb_k, \ub_k}
        \propto \exp \Bigg( - \sum_{m\in \overline{\ensuremath{\mathbb{E}}}_k} \left( g_m(z_m) + \phi_{m,\alpha}( v_m  , z_m - u_m ) \right) \Bigg)
\end{align}
where, for every $\xb \in {\ensuremath{\boldsymbol{\mathsf H}}}$, $\vb_k = (v_m)_{m \in \overline{\ensuremath{\mathbb{E}}}_k} = \rho\left(\begin{matrix}
( \Db_{m,k} \xb_k )_{m \in \ensuremath{\mathbb{E}}_k}  \\
( \Db_{m,\overline{\ensuremath{\mathcal{W}}}_m} \xb_{\overline{\ensuremath{\mathcal{W}}}_m} )_{m \in \ensuremath{\mathbb{E}}_{\ensuremath{\mathcal{R}}_k}}
\end{matrix}\right)$. \\
Then, Algorithm~\eqref{algo:single-synch} is equivalent to Algorithm~\eqref{algo:SPA-simple-sampling}.
\end{proposition}

All the computations described in Algorithm~\ref{algo:single-synch} are conducted simultaneously on each worker, after the necessary communication phases.

Before giving the proof of this result, we want to emphasize that communications to apply the operator $\Db$ and its adjoint $\Db^*$ are symmetric: the same workers are involved in communications, but the direction of the communications (send/receive) is reversed; the indices of the subparts that are communicated are not the same. Communications necessary to computations involving $\Db$ are received by worker $k$ from those of ${\cal R}_k$; communications necessary to computations involving $\Db^*$ are received by worker $k$ from those of ${\cal S}_k$.
\Cref{fig:comp-comm-D}(a) illustrates the communications involved in the application of the operator $\Db$ from \Cref{ex:block-sparse}. \Cref{fig:comp-comm-D}(b) illustrates the communications associated with the application of the adjoint operator $\Db^*$.

\begin{proof}
Using notation from \Cref{Sec:Model}, \eqref{eq:g-add}--\eqref{eq:psi-add} in \Cref{ass:gen} can be directly re-written as in~\eqref{eq:g-dist-simple}--\eqref{eq:psi-dist-simple}. Then, fixing $\vb = \Db\xb$, the conditional distributions~\eqref{eq:dist-SPA-cond:z}-\eqref{eq:dist-SPA-cond:u} can be rewritten as
\begin{align}
    \distc{\alpha}{\zb}{\vb, \ub}
    &   = \prod_{k=1}^K  \distc{k, \alpha}{\zb_k}{\vb_k, \ub_k}, \\
    \distc{(\alpha, \beta)}{\ub}{\vb, \zb}
    &= \prod_{k=1}^K   \distc{k, (\alpha, \beta)}{\ub_k}{\vb_k, \zb_k},
\end{align}
respectively, where, for every $k\in \{1, \ldots, K\}$,
\begin{align}
    \distc{k, \alpha}{\zb_k}{\vb_k, \ub_k}
    &   \propto \exp \left( - \sum_{m\in \overline{\ensuremath{\mathbb{E}}}_k} \left( g_m(z_m) + \phi_{m,\alpha}(v_m  , z_m - u_m ) \right)\right) ,  \\
    \distc{k, (\alpha, \beta)}{\ub_k}{\vb_k, \zb_k}
    &   \propto \exp \left( - \sum_{m\in \overline{\ensuremath{\mathbb{E}}}_k} \left( \phi_{m, \alpha}( v_m, z_m - u_m) - \psi_{m, \beta} (u_m) \right) \right),
    \\
    &= \ensuremath{\mathcal{N}} \bigg( \frac{\beta^2}{\alpha^2+\beta^2} (\zb_k - \vb_k), \frac{\alpha^2+\beta^2}{\alpha^2 \beta^2} \ensuremath{\mathbf{I}}\, \bigg), \nonumber
\end{align}
for $\vb_k = (v_m)_{m \in \overline{\ensuremath{\mathbb{E}}}_k} = ([\Db \xb]_m)_{m \in \overline{\ensuremath{\mathbb{E}}}_k}$.
Hence, for every iteration $t \in \{0, \ldots, T\}$ of the Algorithm~\eqref{algo:SPA-simple}, and for every worker $k \in \{1, \ldots, K\}$,
the sampling of $\zb_k^{(t+1)}$ and $\ub_k^{(t+1)}$ only requires partial information from $\big( \vb_k^{(t+1)}, \ub_k^{(t)} \big) $ and $\big( \vb_k^{(t+1)}, \zb_k^{(t+1)} \big) $, respectively.
Algorithm~\eqref{algo:SPA-simple-sampling} is thus equivalent to
\begin{equation}    \label{prop:synch-SPA-simple:proof:algo1}
\begin{array}{l}
\text{for } t = 0,1, \ldots, T	\\
\left\lfloor
\begin{array}{l}
    \xb^{(t+1)} = \ensuremath{\text{\rm prox}}_{\gamma f}\left( \xb^{(t)} - \gamma \nabla h(\xb^{(t)}) - \gamma \Db^* \Big( \nabla \phi_\alpha(\cdot , \zb^{(t)}-\ub^{(t)})((v_m^{(t)})_{m \in \overline{\ensuremath{\mathbb{E}}}_k}) \Big) + \sqrt{2\gamma} \wb^{(t)} \right) \\
    \text{for } k = 1, \ldots, K	\\
    \left\lfloor
    \begin{array}{l}
        \vb_k^{(t+1)}
        = ([\Db \xb^{(t+1)}]_m)_{m \in \overline{\ensuremath{\mathbb{E}}}_k}, \\
        \zb_k^{(t+1)}
            \sim \trans{\zb_k, \alpha}{\vb_k^{(t+1)}, \zb_k^{(t)}, \ub_k^{(t)}}, \\
        \ub_k^{(t+1)}
            \sim \ensuremath{\mathcal{N}} \bigg( \frac{\beta^2}{\alpha^2+\beta^2} (\zb_k^{(t+1)} - \vb_k^{(t+1)}), \frac{\alpha^2+\beta^2}{\alpha^2 \beta^2} \ensuremath{\mathbf{I}}\, \bigg),
    \end{array}
    \right. \\[0.1cm]
\end{array}
\right.
\end{array}
\end{equation}
where, for every $k\in \{1, \ldots, K\}$, $\trans{\zb_k, \alpha}{ [\Db \xb]_{m \in \overline{\ensuremath{\mathbb{E}}}_k}, \zb_k, \ub_k}$ is a transition kernel with invariant distribution $\distc{k,\alpha}{\zb_k}{[\Db \xb]_{m \in \overline{\ensuremath{\mathbb{E}}}_k}, \ub_k}$.

It remains to show that, for every $t\in \{0, \ldots, T\}$, the computation of $\xb^{(t+1)}$ can be parallelized over $k\in \{1, \ldots, K\}$.
Since $(\ensuremath{\mathbb{V}}_k)_{1 \le k \le K}$ is a partition of $\{1, \ldots, N\}$, then, for every $k\in \{1, \ldots, K\}$, $\xb_k^{(t+1)}$ will be given by the $k$-th vertex values of the random variable generated from the
PSGLA transition kernel.

The gradient of \eqref{eq:phi-dist-simple} with respect to $\xb$ is required to compute $(\xb_k^{(t+1)})_{1 \le k \le K}$.
Using the chain rule on $\xb \in {\ensuremath{\boldsymbol{\mathsf H}}} \mapsto (\phi_\alpha( \cdot, \zb-\ub) \circ \Db) (\xb)$, for $(\zb, \ub)\in {\ensuremath{\boldsymbol{\mathsf G}}}^2$ fixed, we obtain, for every $\xb \in {\ensuremath{\boldsymbol{\mathsf H}}}$,
\begin{equation}
    \nabla_{\xb} \phi_\alpha(\Db \xb , \zb-\ub)
    = \nabla_{\xb} (\phi_\alpha( \cdot, \zb-\ub) \circ \Db) (\xb)
    = \Db^* \nabla \phi_\alpha (\cdot, \zb-\ub) (\Db \xb).
\end{equation}
According to \eqref{eq:phi-dist-simple}, for every $(\vb, \zb) \in {\ensuremath{\boldsymbol{\mathsf G}}}^2$, we have
\begin{align}
    \nabla_{\vb} \phi_\alpha (\vb, \zb)

    = \varrho \left( \left( \begin{matrix}
    \Big( \phi'_{m,\alpha}(v_m, z_m) \Big)_{m \in \ensuremath{\mathbb{E}}_k} \\
    \Big( \phi'_{m, \alpha}(v_m, z_m) \Big)_{m \in \ensuremath{\mathbb{E}}_{\ensuremath{\mathcal{R}}_k}}
    \end{matrix} \right)_{1 \le k \le K} \right),
\end{align}
where $\varrho$ is the permutation operator defined in~\eqref{eq:def:perm}.
For every $k \in \{1, \ldots, K\}$, every $\xb \in {\ensuremath{\boldsymbol{\mathsf H}}}$ and $(\zb, \ub) \in {\ensuremath{\boldsymbol{\mathsf G}}}^2$, let $(d_m)_{m\in \overline{\ensuremath{\mathbb{E}}}_k}$ be defined as
\begin{align}
    d_m
    := \phi'_{m,\alpha}(\cdot , z_m-u_m) (v_m)
    &= \begin{cases}
    \phi'_{m,\alpha}(\cdot , z_m-u_m) (\Db_{m,k} \xb_k),
    &   \text{if } m \in {\ensuremath{\mathbb{E}}}_k,  \\
    \phi'_{m,\alpha}(\cdot , z_m-u_m) (\Db_{m,\overline{\ensuremath{\mathcal{W}}}_m} \xb_{\overline{\ensuremath{\mathcal{W}}}_m}),
    &   \text{if } m \in \ensuremath{\mathbb{E}}_{\ensuremath{\mathcal{R}}_k}.
    \end{cases}

\end{align}
Finally,
\begin{equation}
    \nabla_{\xb} \phi_\alpha(\Db \xb , \zb-\ub) = \Db^* \varrho \Bigg( \Big( (d_m)_{m \in \overline{\ensuremath{\mathbb{E}}}_k} \Big)_{1 \le k \le K} \Bigg). \label{proof:gradient}
\end{equation}
To extract the $k$-th block from this gradient, \eqref{proof:gradient} can be decomposed as follows:
\begin{align}
    \ensuremath{\boldsymbol{\delta}}_k
    &:= [\nabla_{\xb} \phi_\alpha(\Db \xb , \zb-\ub)]_k \nonumber \\
    &= \sum_{m \in \overline{\ensuremath{\mathbb{E}}}_k} \Db_{m,k}^* d_m + \sum_{k' \in \ensuremath{\mathcal{S}}_k} \sum_{m \in \ensuremath{\mathbb{E}}_{(k,k')}} \Db_{m, k}^* d_m.  \label{prop:synch-SPA-simple:sum-grad-glob}
\end{align}
Figure~\ref{fig:comp-comm-D} illustrates this for the simple example described in Figure~\ref{fig:opsplitting}.

The second term in \eqref{prop:synch-SPA-simple:sum-grad-glob} gathers information received from the set of neighbours $\ensuremath{\mathcal{S}}_k$ of worker $k$ that is necessary to computations involving $\Db^*$.
Then, using the notation
\begin{equation}
    (\forall k' \in \ensuremath{\mathcal{S}}_k)\quad
    \widetilde{\db}_{(k,k')} = \sum_{m \in \ensuremath{\mathbb{E}}_{(k,k')}} \Db_{m,k}^* d_m,
\end{equation}
we obtain
\begin{equation} \label{eq:aggregate_for_adjoint}
    \ensuremath{\boldsymbol{\delta}}_k
    = \sum_{m \in \overline{\ensuremath{\mathbb{E}}}_k} \Db_{m,k}^* d_m + \sum_{k' \in \ensuremath{\mathcal{S}}_k} \widetilde{\db}_{(k,k')} ,
\end{equation}
where, for every $k' \in \ensuremath{\mathcal{S}}_k$, $\widetilde{\db}_{(k,k')}$ is computed on worker $k'$ and communicated to worker $k$ to form $\ensuremath{\boldsymbol{\delta}}_k$.
Therefore the computations performed on a given worker $k \in \{1, \ldots, K\}$ include the computation of \textit{local} gradients $(\Db_{m,k}^* d_m)_{m \in \overline{\ensuremath{\mathbb{E}}}_k}$ as well as the gradients to be sent to neighbours, i.e., $(\widetilde{\db}_{(k',k)})_{k' \in \ensuremath{\mathcal{R}}_k}$.
Note that the role of the sets $\ensuremath{\mathcal{R}}_k$ and $\ensuremath{\mathcal{S}}_k$ is exchanged when the adjoint operator $\Db^*$ is considered.

Using the above notation, and using \cite[Prop. 24.11]{Bauschke2017CAMOTHS} on the additively separable functions $h$ and $f$, leads to the conclusion that \eqref{algo:single-synch} is equivalent to \eqref{prop:synch-SPA-simple:proof:algo1}, and therefore to Algorithm~\eqref{algo:SPA-simple-sampling}.
\end{proof}

\begin{figure}[ht]
\centering
\begin{tabular}{cc}\small
    \includegraphics[width=0.48\textwidth]{figures/illustrations/applyD.png}
    &   \includegraphics[width=0.4\textwidth]{figures/illustrations/adjoint_grad_glob.png}  \\% [-0.2cm]
    (a) Computation of $\Db \xb$
    &   (b) Computation of $\Db^* \db$
\end{tabular}

\caption{\label{fig:comp-comm-D}\small
Illustration of the computation of (a) $\Db \xb$ for $\xb \in {\ensuremath{\boldsymbol{\mathsf H}}}$, and (b) $\Db^* \db$ for $\db \in {\ensuremath{\boldsymbol{\mathsf G}}}$. $\Db$ corresponds to the block-sparse matrix from Figure~\ref{fig:opsplitting}, and $\Db^*$ to its adjoint. The yellow, green and blue colours identify how $\Db$, $\Db^*$, $\xb$ and $\db$ are distributed over the workers $k=1$, $k=2$, and $k=3$, respectively.
Communications are illustrated by coloured arrows over subparts of $\xb$ and $\db$, respectively.
Continuous coloured lines on $\xb$ and $\db$ identify vertices and hyperedge weights required to perform local computations. The corresponding blocks in $\Db$ and $\Db^*$ are emphasized accordingly. Vertices and hyperedge weights to be sent from a worker to another are delineated in dashed lines, with the corresponding blocks within $\Db$ and $\Db^*$ highlighted accordingly.
(a) Vertices are on columns, and hyperedges on rows.
(b) Vertices are on rows, and hyperedges on columns. The circled ``plus'' symbol emphasizes that hyperedge weights are aggregated on the reception workers.
}
\end{figure}

\begin{example}\label{ex:comp-comm-D}
Figure~\ref{fig:comp-comm-D}(a) illustrates the computation of $\Db \xb$, for $\xb \in {\ensuremath{\boldsymbol{\mathsf H}}}$ for the block-sparse matrix $\Db$ of Example~\ref{ex:block-sparse}, see also Figure~\ref{fig:blocksparse}. Figure~\ref{fig:comp-comm-D}(b) illustrates the computation of $\Db^* \db$, for $\db \in {\ensuremath{\boldsymbol{\mathsf G}}}$ with some adjoint symmetry. \\
Subparts $\ensuremath{\mathbb{V}}_{(1,2)}$ of $\xb_2$ and $\ensuremath{\mathbb{V}}_{(2,3)}$ of $\xb_3$ need to be communicated to workers $1$ and $2$, respectively, to compute $\Db \xb$.
In Figure~\ref{fig:comp-comm-D}(a), for worker $k=1$, the quantity $(\Db_{m,\overline{\ensuremath{\mathcal{W}}}_m} \xb_{\overline{\ensuremath{\mathcal{W}}}_m})_{m\in \ensuremath{\mathbb{E}}_{\ensuremath{\mathcal{R}}_1}}$ can be computed by multiplying the subpart of $\Db$ corresponding to $(\Db_{m, \overline{\ensuremath{\mathcal{W}}}_m})_{m \in \ensuremath{\mathbb{E}}_{\ensuremath{\mathcal{R}}_1}}$ (see \Cref{fig:comp-comm-D}~(a), horizontal rectangle in continuous yellow lines) with $(\xb_{\overline{\ensuremath{\mathcal{W}}}_m})_{m\in \ensuremath{\mathbb{E}}_{\ensuremath{\mathcal{R}}_1}} = (\xb_1, \xb_{(1,2)})$ (vertical rectangle in continuous yellow lines).
This quantity is computed on worker $1$, once worker $2$ has communicated $\xb_{(1,2)}$ to worker $k=1$.
Similarly, for worker $k=2$, the quantity $(\Db_{m,\overline{\ensuremath{\mathcal{W}}}_m} \xb_{\overline{\ensuremath{\mathcal{W}}}_m})_{m\in \ensuremath{\mathbb{E}}_{\ensuremath{\mathcal{R}}_2}}$ can be computed by multiplying the subpart of $\Db$ corresponding to $(\Db_{m , \overline{\ensuremath{\mathcal{W}}}_m})_{m \in \ensuremath{\mathbb{E}}_{\ensuremath{\mathcal{R}}_2}}$ (horizontal rectangle in continuous green lines) with $(\xb_{\overline{\ensuremath{\mathcal{W}}}_m})_{m\in \ensuremath{\mathbb{E}}_{\ensuremath{\mathcal{R}}_2}} = (\xb_2, \xb_{(2,3)})$ (vertical rectangle in continuous green lines).
This quantity is computed on worker $2$, once worker $3$ has communicated $\xb_{(2,3)}$ to worker $k=2$. \\
To compute $\Db^* \db$, communications are performed after computing subparts of $\Db^* \db$.
In Figure~\ref{fig:comp-comm-D}(b), for worker $k=1$, the quantity $[\Db^* \db]_1$ can be computed by multiplying the subpart of $\Db^*$ corresponding to $(\Db^*_{m, 1})_{m \in \overline{\ensuremath{\mathbb{E}}}_1}$ (continuous-line yellow rectangle, top left) with $(d_m)_{m \in \overline{\ensuremath{\mathbb{E}}}_1}$. For the computation of $[\Db^* \db]_1$, no communication with other workers is needed in this example.
For worker $k=2$, $[\Db^* \db]_2$ needs to be decomposed between parts of $\Db^*$ stored on worker $k=2$, and parts of $\Db$ that are stored on other workers, i.e., $k=1$ for this example. For the parts of $\Db^*$ stored on worker $k=2$, the subpart of $\Db^*$ corresponding to $(\Db^*_{m, 2})_{m \in \overline{\ensuremath{\mathbb{E}}}_2}$ (continuous-line green rectangle) needs to be multiplied with $(d_m)_{m \in \overline{\ensuremath{\mathbb{E}}}_2}$. For the parts of $\Db^*$ stored on $k=1$, the subpart of $\Db^*$ corresponding to $(\Db_{m, 2}^*)_{m \in \ensuremath{\mathbb{E}}_{(1,2)}})$ (dotted-line green rectangle) needs to be multiplied with $(d_m)_{m\in \ensuremath{\mathbb{E}}_{(1,2)}}$ (dotted-line green rectangle on vector $\db$). This second part is computed on worker $k=1$ (using only the yellow part of the dotted-line green rectangle), and then communicated and aggregated on worker $k=2$.
For worker $k=3$, the decomposition of the quantity $[\Db^* \db]_3$ is similar to the one taken for worker $k=2$.
\end{example}

\subsection{Distributed SPMD architecture}
\label{Ssec:algo:dist-spmd}

The distributed block-coordinate sampler described in the previous section can benefit from an implementation on an Single Program Multiple Data (SPMD) architecture.
In contrast with a client-server configuration, all the workers
execute the same task on a subset of each block of parameters and observations~\cite{Darema2001}. This enables the hypergraph structure of $\Db$ to be exploited to reduce the number and volume of the communications.

In practice, \Cref{mod:hypergraph} and \Cref{ass:gen} ensure that most of the operations of Algorithm~\eqref{algo:SPA-simple-sampling} are compatible with an SPMD architecture.
As noted in~\Cref{Sec:introduction}, the data fidelity term can be represented by either the function $f$, $g$ of $h$ in~\eqref{eq:dist-gen-lin}, depending on the nature of the noise $\ensuremath{\mathcal{D}}$ affecting the observations $\yb \in {\ensuremath{\boldsymbol{\mathsf G}}}$ in~\eqref{eq:full_model}. The separability Assumptions~\ref{ass:gen}-\ref{ass:gen:ii} and~\ref{ass:gen:iii} ensure that $K$ conditionally independent blocks of observations can be formed. This leads to a partition of the observations into $\yb = (\yb_k)_{1 \leq k \leq K}$, with $\yb_k \in {\ensuremath{\boldsymbol{\mathsf G}}}_k$ the block of observations stored on the worker $k$. In this case, no observation needs to be exchanged between the workers. The separability of $g$ further implies that the evaluation of its proximity operator is easy to parallelize. Combining the definition of $\phi_{\alpha}$ and $\psi_{\beta}$ (\Cref{ass:gen}-\ref{ass:gen:iv}) with the structure of $\Db$ finally enables Algorithm~\eqref{algo:SPA-simple-sampling} to be reformulated using blocks of parameters, each stored on (only) one of the $K$ workers. Using an SPMD architecture for Algorithm~\eqref{algo:single-synch} offers several advantages listed below.

\begin{enumerate}
    \item Load balancing: computation costs can be equally shared among the workers, as they all operate similar tasks on a subset of (overlapping) parameters and observations.
    %
    \item Parallelization flexibility: an SPMD architecture can be readily used to address~\eqref{eq:dist-gen-lin} under \Cref{ass:gen}, whereas a client-server architecture cannot (no conditionally independent blocks of variables as in~\cite{Rendell2021})
    %
    \item Memory and computing costs per worker: each worker can be assigned a block of observations and a corresponding block of parameters. This opportunity can significantly reduce the computing and memory costs per worker.
    %
    \item Communications and data locality: all the observations and most of the parameters required to perform operations on a worker can be directly stored on the same worker, and do not require to be communicated (\emph{data locality}). The conditions given in~\Cref{Sssec:Model:cond-eff-dist} guarantee that communication costs are limited: only a few elements need to be retrieved by each worker from a small number of connected workers.
\end{enumerate}

\Cref{appendix:SPA-PSGLA-multi} addresses the multi-term extension of~\eqref{eq:dist-gen-lin}. Note that a client-server architecture can also be used to address this case. However, it cannot take advantage of the structure of the hypergraphs underlying the linear operators $(\Db_i)_{1 \leq i \leq I}$. The number of workers it can accommodate is also restricted to the number of conditionally independent blocks of variables, as in~\cite{Vono_etal_2020, Rendell2021}.

\subsection{Problem statement}
\label{Ssec:application:problem_statement}

Supervised Poisson deconvolution aims at inferring an unknown image $\overline{\xb} \in \ensuremath{\mathbb{R}}^N$ from observations $\yb = (y_m)_{1\le m \le M} \in \ensuremath{\mathbb{R}}^M$ such that
\begin{equation} \label{eq:model_convolution}
    (\forall m \in \{1, \ldots, M\}) \quad
    y_m \sim \ensuremath{\mathcal{P}} ( [\Db_1 \overline{\xb}]_m ),
\end{equation}
where $\Db_1 \in \ensuremath{\mathbb{R}}^{M \times N}$ is a discrete convolution operator derived from a kernel of size $L = L_1 \times L_2 \ll N$, and $\ensuremath{\mathcal{P}}(\mu)$ is a Poisson distribution with mean $\mu$.
In this context, $\overline{N}=N$ and, for every $n\in \{1, \ldots, N\}$, $\ensuremath{{\mathsf H}}_n = \ensuremath{\mathbb{R}}$. Similarly, $\overline{M}=M$ and, for every $m\in \{1, \ldots, M\}$, $\ensuremath{{\mathsf G}}_m = \ensuremath{\mathbb{R}}$.

We propose to solve this problem with
a hybrid prior combining a non-negativity constraint and a TV regularization~\cite{Rudin1992}. Such prior has for instance been considered in \cite{Figueiredo2010,Vono2019icassp}.
The resulting posterior distribution is given by
\begin{equation}
    \dist{}{\xb} \propto \exp \left( -f(\xb) - g_1(\Db_1 \xb) - g_2(\Db_2 \xb) \right),
\end{equation}
where
$g_1\colon \mathbb{R}^M \to ]-\infty, +\infty] \colon (z_m)_{1 \le m \le M}  \mapsto \sum_{m=1}^M g_{1,m}(z_m)$ is the data-fidelity term given by
\begin{align}
    &(\forall m \in \{1, \ldots, M\}) \quad
    g_{1,m}(z_m) =  - y_m \log( z_m ) + z_m, \label{eq:data_fidelity}
\end{align}
due to the Poisson distribution, $f =\iota_{[0,+\infty[^N} \colon \ensuremath{\mathbb{R}}^N \to ]-\infty, +\infty]$ is the indicator function of the positive orthant, and $g_2 \circ \Db_2$ models the discrete isotropic TV~\cite{Rudin1992}.
Precisely, $\Db_2 \colon \ensuremath{\mathbb{R}}^N \to \ensuremath{\mathbb{R}}^{2\times N}$ is the concatenation of the vertical and horizontal discrete gradients, and $g_2: \mathbb{R}^{2 \times N} \to ]-\infty, +\infty]$ is the $\ell_{2,1}$-norm
\begin{equation}
    \begin{split}
        &(\forall \zb = (\zb_n)_{1\leq n \leq N} \in \ensuremath{\mathbb{R}}^{2 \times N}) \quad g_2(\zb) = \sum_{n=1}^{N} g_{2, n} (\zb_n), \\
        &(\forall n \in \{1, \ldots, N\}) (\forall \zb_n \in \ensuremath{\mathbb{R}}^2) \quad
        g_{2, n}(\zb_n) = \kappa \| \zb_n \|_2, \quad \text{where }\kappa > 0.
    \end{split}
\end{equation}

As a first step towards a distributed implementation, the AXDA approach (see~\Cref{appendix:SPA-PSGLA-multi} for details) is used to approximate $ \dist{}{\xb} $ by
\begin{multline} \label{eq:posterior_deconvolution}
     \dist{\ensuremath{\boldsymbol{\alpha}}, \ensuremath{\boldsymbol{\beta}}}{\xb,(\zb_i, \ub_i)_{1 \le i \le 2}}
     \!\propto\! \exp \bigg( \!\! - f(\xb) - \sum_{i=1}^2 \big( g_i(\zb_i) + \phi_{i, \alpha_i}(\Db_i \xb, \zb_i-\ub_i)  + \psi_{i, \beta_i}(\ub_i) \big) \bigg),
\end{multline}
where $(\alpha_1, \alpha_2, \beta_1, \beta_2) \in ]0, +\infty[^4$, and the functions $(\phi_{i, \alpha_i}, \psi_{i, \beta_i})_{1 \le i \le 2}$ are defined in~\eqref{eq:phi-psi-def}.

The directed acyclic graph reported in~\Cref{fig:dag} summarizes the structure of the approximate posterior distribution~\eqref{eq:posterior_deconvolution}, highlighting conditional dependencies between variables. Note that $(\zb_1, \ub_1)$ and $(\zb_2, \ub_2)$ are the only conditionally independent blocks of variables.

\begin{figure}
    \centering
    \includegraphics[keepaspectratio, width=0.45\textwidth]{figures/illustrations/dag.pdf}

    \vspace*{-0.5cm}

    \caption{\small
    Directed acyclic graph describing the factorization of the approximate posterior distribution~\eqref{eq:posterior_deconvolution}. Circled variables correspond to

    random variables. Other variables are fixed \emph{a priori}.}
    \label{fig:dag}
\end{figure}

\subsection{Proposed SPMD implementation}
\label{Ssec:application:spmd}

The proposed approach is applicable when a single operator is involved in the model (as in Section~\ref{Sec:dist-single-algo}), but also when multiple operators are involved (see Appendix~\ref{appendix:SPA-PSGLA-multi}), as in the present example.
Note that a client-server architecture could also be considered. However, it would drastically increase the communication costs, as full-size variables would need to be duplicated on all workers or exchanged (see \Cref{Ssec:algo:dist-spmd} for more details).

We propose to use the proposed distributed SPA Gibbs sampler~\eqref{algo:multi-synch} with $I=2$ linear operators. In addition, for every $t\in \{0, \ldots, T\}$ and $k\in \{1, \ldots, K\}$, we use PSGLA transitions to compute $(\zb_{1, k}^{(t+1)}, \zb_{2, k}^{(t+1)})$, while $(\ub_{1, k}^{(t+1)}, \ub_{2, k}^{(t+1)})$ are drawn from their conditional distribution. More precisely, we have, for $i \in \{1, 2\}$,
\begin{align}
    \zb_{i, k}^{(t+1)}
    &=   \ensuremath{\text{\rm prox}}_{\eta_i, g_{i,k}} \left( \zb_{i,k}^{(t)} - \eta_i\alpha_i^{-2} \big( \zb_{i,k}^{(t)} - \vb_{i,k}^{(t+1)} + \ub_{i,k}^{(t)} \big) + \sqrt{2\eta_i} \ensuremath{\boldsymbol{\xi}}_{i,k}^{(t)} \right) ,
\end{align}
where $g_{i, k} = \sum_{m\in \overline{\ensuremath{\mathbb{E}}}_{i,k}} g_{i,m}$, $\eta_i \in ]0, \alpha_i^{-2}[$, and $\xi_{i,k}^{(t)} \sim \ensuremath{\mathcal{N}}(0, \ensuremath{\mathbf{I}}\,)$.
The proximity operators involved in the resulting algorithm can be found, e.g., in~\cite{Figueiredo2010,Komodakis2015}.
The distributed implementation of $(\Db_i)_{1 \leq i \leq 2}$ and $(\Db_i^*)_{1 \leq i \leq 2}$ is detailed below, using a 2D Cartesian grid of $K$ workers.

\begin{figure
    \centering
    \begin{tabular}{@{}c@{}c@{}}
    \includegraphics[keepaspectratio, width=0.5\textwidth]{figures/illustrations/broadcast2d_sync_full.pdf}
    &	\includegraphics[keepaspectratio, width=0.5\textwidth]{figures/illustrations/aggregate2d_sync_full.pdf}	\\[-0.2cm]
    (a) distributed implementation of $\Db_1$
    &    (b) distributed implementation of $\Db_1^*$
    \end{tabular}

    \vspace*{-0.2cm}
    \caption{\small
    Communication patterns involved in the distributed implementation of (a) $\Db_1$ and (b) $\Db_1^*$.
    Each worker $k$ is required to communicate with two contiguous workers along each direction successively.
    Colored arrows indicate whether vertices (a) or hyperedge weights (b) are received from or sent to a nearby worker. Vertices (a) and hyperedge weights (b) sent and received along the horizontal and vertical communication steps are highlighted in blue and red, respectively. (b) The circled ``plus'' symbols represent contributions aggregated with the corresponding hyperedge weights of the reception worker. These operations correspond to the second term in~\eqref{eq:aggregate_for_adjoint}.
    }
    \label{fig:app:comm2D}
\end{figure}

\paragraph{Distributed implementation of $\Db_1$}

\Cref{fig:app:comm2D} illustrates the preliminary communications required by the distributed implementation of $\Db_1$. Each worker $k\in \{1, \ldots, K\}$ needs to collect $(\xb_{1,(k,k')})_{k' \in \ensuremath{\mathcal{R}}_k}$, as illustrated in~\Cref{fig:app:comm2D}(a). A first communication step occurs along the horizontal axis. Each worker $k$ sends vertices from its left-hand-side border (of width $L_2$, the horizontal width of the blur kernel) to its neighbour on the left (\Cref{fig:app:comm2D}(a), dark blue areas), and receives vertices from the neighbour on its right (\Cref{fig:app:comm2D}(a), light blue areas). Once the first step is complete, a second communication step occurs along the vertical axis with the top and bottom neighbours of the worker $k$ (\Cref{fig:app:comm2D}(a), dark and light red areas, of width $L_1$). The local operator $(\Db_{1,m,\overline{\ensuremath{\mathcal{W}}}_{1,m}})_{ m \in  \ensuremath{\mathbb{E}}_{1,\ensuremath{\mathcal{R}}_k} }$ is equivalent to considering a convolution matrix and a selection operator. The latter ensures that  the correct boundaries are considered (only the convolution outputs which have not interacted with the boundaries of $(\xb_{\overline{\ensuremath{\mathcal{W}}}_{1, m}})_{ m \in  \ensuremath{\mathbb{E}}_{\ensuremath{\mathcal{R}}_k} }$ are valid).

\paragraph{Distributed implementation of $\Db_1^*$}

The distributed implementation of $\Db_1^*$ is similar to that of $\Db_1$. A first communication step occurs along the horizontal axis. Each worker $k$ sends hyperedge weights from its right-hand-side border (of width $L_2$) to its neighbour on the right (\Cref{fig:app:comm2D}(b), dark blue areas), and aggregates weights from the neighbour on its left (\Cref{fig:app:comm2D}(b), light blue areas). A similar communication step then occurs along the vertical axis with the bottom and top neighbours of the worker $k$ (\Cref{fig:app:comm2D}(b), dark and light red areas, of width $L_1$). A convolution is then applied on each worker to all the weights locally available, followed by a selection operator to ensure correct boundaries have been used.

\paragraph{Distributed implementation of $\Db_2$}

The distributed implementation of $\Db_2$ requires the same communication pattern as $\Db_1$, successively exchanging messages with width 1 along the horizontal and vertical directions. The operator $(\Db_{2,m,\overline{\ensuremath{\mathcal{W}}}_{2,m}})_{ m \in  \ensuremath{\mathbb{E}}_{2,\ensuremath{\mathcal{R}}_k} }$ corresponds to a local discrete gradient operator, using the boundaries retrieved during the communication step. Note that, for $k \in \{1, \dotsc, K\}$ and $k' \in \ensuremath{\mathcal{R}}_k$, the set of vertices $\ensuremath{\mathbb{V}}_{2, (k',k)}$ to be communicated by the worker $k$ to $k'$ is such that $\ensuremath{\mathbb{V}}_{2, (k,k')} \subset \ensuremath{\mathbb{V}}_{1, (k,k')}$. The total number of elements to be communicated in the algorithm is thus reduced, as vertices required by the distributed implementation of $\Db_2$ already need to be communicated for $\Db_1$.

\paragraph{Distributed implementation of $\Db_2^*$}

The distributed implementation of $\Db_2^*$ requires the same communication pattern as $\Db_1^*$, successively exchanging messages with width 1 along the horizontal and vertical directions. A local discrete adjoint gradient operator is then applied on each worker $k$ to all the hyperedge weights available locally.

\subsection{Experiments}
\label{Ssec:application:experiments}

The proposed approach is evaluated in terms of estimation quality and scalability on the deconvolution problem of~\Cref{Ssec:application:problem_statement}. Results are compared with those of the reference serial SPA Gibbs sampler from~\cite{Vono2019icassp}.

\subsubsection{Simulation setting} \label{Ssec:experiments:settings}

All the experiments have been conducted on a single computer equipped with two 2.1~GHz, 18-core, Intel Xeon E5-2695 v4 series processors (36 CPU cores in total). In this setting, a worker corresponds to a process running on one CPU core. The proposed distributed sampler has been implemented in Python using the \texttt{mpi4py} library~\cite{Dalcin2021}. A short code to reproduce part of the experiments is available at \url{https://gitlab.cristal.univ-lille.fr/pthouven/dspa} and will be completed upon acceptance of the paper.

The proposed approach is compared with the SPA Gibbs sampler proposed in~\cite{Vono2019icassp}. The latter relies on a different splitting strategy compared to \Cref{Ssec:application:problem_statement} using $I = 3$ operators (see~\cite{Vono2019icassp} for further details).
In practice, MYULA transition kernels are leveraged to sample from conditional distributions involving non-smooth potential functions. This choice of transition kernel and splitting strategy requires the proximity operator of the TV norm to be evaluated at each iteration of the sampler, using a primal-dual algorithm~\cite{Chambolle2011}.

Performance assessment is conducted in terms of runtime per iteration and quality of both the \emph{minimum mean square error} (MMSE) and \emph{maximum a posteriori} (MAP) estimators. The estimators are denoted $\xb_{\text{MMSE}}$ and $\xb_{\text{MAP}}$, respectively. Reconstruction quality is quantified through the structural similarity index (SSIM)~\cite{Wang2004} and the signal-to-noise ratio (SNR)
expressed in dB. The associated 95\% confidence intervals are also reported.

The sampler~\cite{Vono2019icassp} has been systematically applied with $((\alpha_i^2, \beta_i^2)_{1 \leq i \leq 3}, \kappa) = \mathbf{1}_7$, where $\mathbf{1}_Q \in \mathbb{R}^Q$ is a vector with entries all equal to 1. The proposed approach uses $((\alpha_i^2, \beta_i^2)_{1 \leq i \leq 2}, \kappa) = \mathbf{1}_5$. For both algorithms, $N_\text{MC} = 5 \times 10^3$ samples have been generated to form $\xb_{\text{MMSE}}$ and $\xb_{\text{MAP}}$ and the 95\% credibility intervals, discarding $N_\text{bi} = 2 \times 10^3$ burn-in samples. Two series of experiments have been considered to assess the estimation quality and scalability.

\subsubsection{Experiment results} \label{Ssec:experiments:results}

\paragraph*{Estimation quality}

Using $K=1$ worker, the proposed approach is compared with~\cite{Vono2019icassp}. Ground truth images with different maximum intensity levels $\overline{x}_{\max} = \max_{1 \leq n \leq N} \overline{x}_n$ have been considered. The values $\overline{x}_{\max} \in \{20, 30\}$ have been adopted for the following datasets:
\begin{enumerate}
    \item \texttt{house} image ($N = 256^2$), using a normalized Gaussian kernel of size $L \in \{3^2, 7^2\}$;
    \item \texttt{peppers} image ($N = 512^2$) using a normalized Gaussian kernel of size $L \in \{7^2, 15^2\}$.
\end{enumerate}

The results reported in \Cref{tab:results_serial} show that the estimators formed with the proposed approach have higher quality metrics compared to the method proposed in~\cite{Vono2019icassp}. In addition, the computing time required by the proposed sampler is between 1.5 and 2 times smaller than~\cite{Vono2019icassp}. This discrepancy comes from the difference in the splitting strategies considered by the two methods. In particular, the splitting approach proposed in~\cite{Vono2019icassp} requires the evaluation of the proximal operator of the TV norm, obtained as the output of an iterative optimization algorithm. The difference in splitting can also affect the quality of the resulting AXDA approximation, as can be observed by the difference in quality of the estimator produced. The MMSE estimator reported in \Cref{fig:results:peppers30} for \cite{Vono2019icassp} appears much smoother compared to the proposed algorithm for the same value of the regularization parameter $\kappa$.
Note that the uncertainty level of the proposed approach is slightly lower and appears more diffuse than~\cite{Vono2019icassp}.

\begin{table}
    \centering

    \resizebox{0.98\textwidth}{!}{%
    \begin{tabular}{lllrrrrrrrr} \toprule
        & Dataset & Algo. & $\SNR(\xb_{\text{MMSE}})$ & $\SNR(\xb_{\text{MAP}})$ & $\SSIM(\xb_{\text{MMSE}})$ & $\SSIM(\xb_{\text{MAP}})$ & Time per iter.  & Runtime \\
        &&&&&&& ($\times 10^{-2}$ s) & ($\times 10^{2}$ s) \\ \midrule
        \multirow{9}{*}{\rotatebox{90}{$\overline{x}_{\max} = 20$}} & House & \cite{Vono2019icassp} & 18.37 & 15.08 & 0.60 & 0.20 & 5.59 (0.12) & 1.68 \\
        & ($L=3^2$) & Proposed & 20.21 & 16.18 & 0.60 & 0.26 & \textbf{2.49 (0.09)} & \textbf{0.75} \\ \cmidrule{2-9}
        %
        & House & \cite{Vono2019icassp} & 18.00 & 14.92 & 0.60 & 0.20 & 8.52 (0.22) & 2.56 \\
        & ($L=7^2$) & Proposed & 19.86 & 15.98 & 0.59 & 0.24 & \textbf{4.88 (0.13)} & \textbf{1.47} \\ \cmidrule{2-9}
        %
        & Peppers & \cite{Vono2019icassp} & 18.98 & 15.08 & 0.66 & 0.25 & 16.71 (0.38) & 5.01 \\
        & ($L=7^2$) & Proposed & 20.52 & 16.07 & 0.67 & 0.30 & \textbf{7.73 (0.33)} & \textbf{2.32} \\ \cmidrule{2-9}
        %
        & Peppers & \cite{Vono2019icassp} & 18.90 & 15.07 & 0.66 & 0.25 & 36.68 (0.33) & 11.01 \\
        & ($L=15^2$) & Proposed & 20.52 & 16.04 & 0.66 & 0.30 & \textbf{23.32 (0.40)} & \textbf{7.00} \\

        \midrule\midrule

        \multirow{9}{*}{\rotatebox{90}{$\overline{x}_{\max} = 30$}} & House & \cite{Vono2019icassp} & 18.23 & 16.46 & 0.64 & 0.29 & 4.71 (0.10) & 1.41 \\
        & ($L=3^2$) & Proposed & 20.18 & 17.84 & 0.66 & 0.34 & \textbf{2.03 (0.05)} & \textbf{0.61} \\ \cmidrule{2-9}
        %
        & House & \cite{Vono2019icassp} & 17.89 & 16.25 & 0.64 & 0.28 & 7.01 (0.10) & 2.10 \\
        & ($L=7^2$) & Proposed & 19.80 & 17.55 & 0.65 & 0.33 & \textbf{3.93 (0.10)} & \textbf{1.18} \\ \cmidrule{2-9}
        %
        & Peppers & \cite{Vono2019icassp} & 19.03 & 16.77 & 0.69 & 0.35 & 17.04 (0.45) & 5.11 \\
        & ($L=7^2$) & Proposed & 20.71 & 17.96 & 0.71 & 0.40 & \textbf{10.26 (1.22)} & \textbf{3.08} \\ \cmidrule{2-9}
        %
        & Peppers & \cite{Vono2019icassp} & 19.00 & 16.72 & 0.69 & 0.35 & 36.81 (0.44) & 11.04 \\
        & ($L=15^2$) & Proposed & 20.74 & 17.98 & 0.71 & 0.41 & \textbf{23.22 (0.34)} & \textbf{6.97} \\
        \bottomrule
    \end{tabular}
    }
    \caption{\small
    Comparison between~\cite{Vono2019icassp} and the proposed approach with $K = 1$. Datasets have been generated from ground truth images with different maximum intensity $\overline{x}_{\max} \in \{20, 30\}$ and convolution kernel sizes $L \in \{3, 7, 15\}$. Results are reported in terms of estimation quality, runtime per iteration and total runtime.}
    \label{tab:results_serial}
\end{table}

\begin{figure}[htbp]
    \centering\small
    \begin{subfigure}[t]{0.31\textwidth}
         \centering
         \includegraphics[keepaspectratio, width=0.95\textwidth]{figures/results/beta_fixed/peppers_M30_ks15_obs.pdf}
         \includegraphics[keepaspectratio, width=0.95\textwidth]{figures/results/beta_fixed/peppers_M30_ks15_ground_truth.pdf}
         \caption{Observations and \\ground truth}
     \end{subfigure}
     %
     \begin{subfigure}[t]{0.31\textwidth}
         \centering
         \includegraphics[keepaspectratio, width=0.95\textwidth]{figures/results/beta_fixed/peppers_M30_ks15_pmyula_mmse.pdf}
         \includegraphics[keepaspectratio, width=0.95\textwidth]{figures/results/beta_fixed/peppers_M30_ks15_pmyula_95p_ci.pdf}
         \caption{MMSE estimator and 95\% credible interval \cite{Vono2019icassp}}
     \end{subfigure}
     %
     \begin{subfigure}[t]{0.31\textwidth}
         \centering
         \includegraphics[keepaspectratio, width=0.95\textwidth]{figures/results/beta_fixed/peppers_M30_ks15_psgla_mmse.pdf}
         \includegraphics[keepaspectratio, width=0.94\textwidth]{figures/results/beta_fixed/peppers_M30_ks15_psgla_95p_ci.pdf}
         \caption{MMSE estimator and 95\% credible interval (prop., $K=1$)}
     \end{subfigure}

    \vspace{-0.5cm}

    \caption{\small
    Estimators and 95\% credibility intervals for \texttt{peppers}, with $\overline{x}_{\max} = 30$ and $L=15^2$.}
    \label{fig:results:peppers30}
\end{figure}

\paragraph*{Strong scaling experiment}
For this experiment, the behaviour of the proposed method is investigated with a varying number of workers $K \in \{1, 2, 4, 8, 16, 32\}$, using the maximum intensity level $\overline{x}_{\max} = 30$ for the \texttt{peppers} dataset.

\Cref{tab:results_strong_scaling_small_kernel} and \Cref{tab:results_strong_scaling_large_kernel} show the speedup for the proposed approach when $L=7^2$ and $L=15^2$, respectively.
In \Cref{tab:results_strong_scaling_large_kernel}, a close to ideal speedup (\emph{i.e.}, close to the number of workers $K$) is observed as $L=15^2$. Note that speedup factors larger than the number of cores may occur, depending on the cache state of the machine at the time the experiments have been run. In comparison, the ideal runtime per iteration (i.e., neglecting communication time) for~\cite{Vono2019icassp} for a client-server architecture with $K = 4$ cores is $1.48 \times 10^{-1}$ s when $L=7^2$ ($1.04 \times 10^{-1}$ s for TV-related terms), and $3.41 \times 10^{-1}$ s when $L=15^2$ ($1.02 \times 10^{-1}$ s for TV-related terms). In both cases, the time required to update the other splitting variables is almost 10 times lower. For this experiment, a client server approach thus leads to a limited runtime performance, given the heterogeneity in the complexity of the tasks assigned to workers. All these results illustrate the ability of the proposed sampler to provide estimators at a fraction of the runtime of the serial implementation by using an increasing number of cores $K$.

\paragraph*{Weak scaling experiment}
For this experiment, the behaviour of the proposed method is investigated when both the size of the problem and the number of workers considered simultaneously increase.
A fixed problem size per worker is considered, using $K \in \{1, 4, 16\}$. Datasets derived from upsampled versions of the \texttt{house} image are considered, with a maximum intensity level $\overline{x}_{\max} = 30$, using $(N, L) \in \{(256^2, 3^2),$ $(512^2, 7^2), (1022^2, 11^2)\}$. These numbers ensure that both the problem size $M$ and the associated number of workers $K$ evolve in the same proportions over the configurations tested.

\Cref{tab:results_weak_scaling} shows scaled speedup factors (that is, normalized by the factor of increase for $K$ and $M$) close to the number of cores used. Differences from a linear scaling may result from fixed communication costs, representing a larger cost per iteration as $L$ increases. Overall, the results illustrate the runtime stability of the approach for a fixed problem size per worker. This experiment efficiently processes a 1 million pixel image in about 2 minutes to obtain an estimator with the associated credibility intervals.

\begin{table*}[!]
    \centering

    \resizebox{0.98\textwidth}{!}{%
    \begin{tabular}{crrrrrrr} \toprule
        $K$ & $\SNR(\xb_{\text{MMSE}})$ & $\SNR(\xb_{\text{MAP}})$ & $\SSIM(\xb_{\text{MMSE}})$ & $\SSIM(\xb_{\text{MAP}})$ & Time per iter.  & Speedup & Runtime  \\
        &&&&& ($\times 10^{-2}$ s) && ($\times 10^{2}$ s) \\ \midrule
        1 & 20.71 & 17.96 & 0.71 & 0.40 & 10.26 (1.22) & \textbf{1.00} & 3.08 \\
        2 & 20.70 & 17.94 & 0.71 & 0.40 & 5.33 (0.10) & \textbf{1.93} & 1.60 \\
        4 & 20.72 & 17.93 & 0.71 & 0.40 & 4.35 (0.12) & \textbf{2.36} & 1.30 \\
        8 & 20.72 & 17.95 & 0.71 & 0.40 & 2.51 (0.20) & \textbf{4.08} & 0.75 \\
        16 & 20.73 & 17.97 & 0.71 & 0.41 & 1.23 (0.03) & \textbf{8.32} & 0.37 \\
        32 & 20.71 & 17.91 & 0.71 & 0.40 & \textbf{0.60 (0.06)} & \textbf{17.04} & \textbf{0.18} \\
      \bottomrule
    \end{tabular}
    }
    \caption{\small
    Results of the strong scaling experiment using a dataset with $\overline{x}_{\max} = 30$ and kernel size $L=7^2$. Performance is reported in terms of estimation quality, time per iteration and speedup.}
    \label{tab:results_strong_scaling_small_kernel}
\end{table*}

\begin{table*}[!]
    \centering

    \resizebox{0.98\textwidth}{!}{%
    \begin{tabular}{crrrrrrr} \toprule
        $K$ & $\SNR(\xb_{\text{MMSE}})$ & $\SNR(\xb_{\text{MAP}})$ & $\SSIM(\xb_{\text{MMSE}})$ & $\SSIM(\xb_{\text{MAP}})$ & Time per iter.  & Speedup & Runtime  \\
        &&&&& ($\times 10^{-2}$ s) && ($\times 10^{2}$ s) \\ \midrule
        1 & 20.74 & 17.98 & 0.71 & 0.41 & 23.22 (0.34) & \textbf{1.00} & 6.97 \\
        2 & 20.72 & 17.94 & 0.71 & 0.40 & 12.42 (0.11) & \textbf{1.87} & 3.73 \\
        4 & 20.74 & 17.98 & 0.71 & 0.41 & 3.22 (0.10) & \textbf{7.22} & 0.97 \\
        8 & 20.74 & 17.97 & 0.71 & 0.41 & 1.91 (0.06) & \textbf{12.17} & 0.57 \\
        16 & 20.72 & 17.95 & 0.71 & 0.40 & 1.31 (0.11) & \textbf{17.73} & 0.39 \\
        32 & 20.73 & 17.97 & 0.71 & 0.40 & \textbf{0.70 (0.07)} & \textbf{33.35} & \textbf{0.21} \\
      \bottomrule
    \end{tabular}
    }
    \caption{\small
    Results of the strong scaling experiment using a dataset with $\overline{x}_{\max} = 30$ and kernel size $L = 15^2$. Performance is reported in terms of estimation quality, time per iteration and speedup.
    }
    \label{tab:results_strong_scaling_large_kernel}
\end{table*}

\begin{table*}[!]
    \centering

    \resizebox{0.98\textwidth}{!}{%
    \begin{tabular}{@{}c@{}rrrrrrrr@{}} \toprule
        $(N, L,K)$  & $\SNR(\xb_{\text{MMSE}})$ & $\SNR(\xb_{\text{MAP}})$ & $\SSIM(\xb_{\text{MMSE}})$ & $\SSIM(\xb_{\text{MAP}})$ & Time per iter.  & Sc. speedup & Runtime  \\
        &&&&& ($\times 10^{-2}$ s) && ($\times 10^{2}$ s) \\ \midrule
        $(256^2, 3^2, \mathbf{1})$  & 20.18 & 17.84 & 0.66 & 0.34 & 2.03 (0.05) & \textbf{1.00} & \textbf{0.61} \\
        $(512^2, 7^2, \mathbf{4})$  & 23.86 & 19.56 & 0.74 & 0.34 & 2.13 (0.08) & \textbf{3.81} & \textbf{0.64} \\
        $(1022^2, 11^2, \mathbf{16})$ & 27.55 & 20.65 & 0.81 & 0.33 & 3.56 (0.18) & \textbf{9.13} & \textbf{1.07} \\
      \bottomrule
    \end{tabular}
    }

    \caption{\small
    Results of the weak scaling experiment using a dataset with $\overline{x}_{\max} = 30$. The reconstruction quality of the estimators is reported with the time per iteration, the scaled speedup and the runtime.}
    \label{tab:results_weak_scaling}
\end{table*}