# On Compressed Sensing of Binary Signals for the Unsourced Random Access Channel
## 2 Compressed Sensing of Binary Signals
### MCMC Algorithm for Recovery
The idea is to construct an ergodic Markov chain such that: (i) its stationary distribution is the desired (high-dimensional) distribution one would like to sample from, namely \(\operatorname{Pr}(\mathbf{x}\,|\,\mathbf{y})\); and (ii) the chain is easy to propagate in time (e.g., its update rule is local). Having constructed such a chain, and assuming that it mixes sufficiently fast (which is often difficult to ensure), one can therefore efficiently sample from the desired distribution, up to high precision. For further background and discussion on MCMC, we refer to the work of [PERSON] and [PERSON] [61] (Chapter 3). The use of MCMC methods for solving inverse problems in signal processing and for decoding/detection in communication is by no means novel (see, e.g., [62; 63; 64; 65; 66; 67; 68; 69]). While both the idea of using LDPC codes as sensing matrices and the idea of using MCMC methods for decoding are not new, our innovation here is in combining the two concepts for the compressed sensing of binary signals problem. As becomes evident below, the sparse structure of the sensing matrix constructed from an LDPC code significantly reduces the computational load from the MCMC decoder by reducing the computational cost of each iteration. We propose using the well-known _Gibbs sampling_ method, also known as _Glauber dynamics_, which is a general-purpose recipe for sampling from high-dimensional distributions. Let \(\mathbb{Q}(\mathbf{x})\) be a distribution over \(\{0,1\}^{M}\) from which one wants to sample; in our case, of course, \(\mathbb{Q}(\mathbf{x})=\Pr(\mathbf{x}\,|\,\mathbf{y})\). We construct a chain \(\mathbf{x}^{(1)},\mathbf{x}^{(2)},\ldots\in\{0,1\}^{M}\) starting from some (arbitrary) initial state \(\mathbf{x}^{(0)}\) according to the following transition rule. Suppose that the current state is \(\mathbf{x}^{(t)}\); one samples a coordinate to update at random, \(i_{t}\sim\text{Uniform}(\{1,\ldots,M\})\), so that \(\mathbf{x}_{j}^{(t+1)}=\mathbf{x}_{j}^{(t)}\) for all \(j\
eq i_{t}\). As for coordinate \(i_{t}\), it is sampled according to the conditional distribution of \(\mathbf{x}_{i_{t}}\), with all other coordinates fixed and given by \(\mathbf{x}^{(t)}\), that is: \(\mathbf{x}_{i_{t}}^{(t+1)}\sim\mathbb{Q}(x_{i_{t}}\,|\,\mathbf{x}_{\sim i_{t}}=\mathbf{x}_ {\sim i_{t}}^{(t)})\) (we denote the vector of all coordinates, except for \(i_{t}\), by \(\mathbf{x}_{\sim i_{t}}\)). Applied to the posterior in (5), Glauber dynamics reads as follows:  
It is easy to see that the process \(\mathbf{x}^{(1)},\mathbf{x}^{(2)},\ldots\in\{0,1\}^{M}\) is an ergodic Markov chain, and therefore has a unique stationary distribution. Furthermore, it is easy to verify that \(\mathbb{Q}(\mathbf{x})\) is a stationary distribution of this chain. Thus, for \(T\) sufficiently large, we have that indeed \(\mathbf{x}^{(T)}\) is distributed as a random sample from \(\mathbb{Q}(\mathbf{x})\). Note that, when \(A\) is a sparse LDPC matrix, each iteration of Algorithm 1 is computationally very cheap. One can easily keep track of \(\mathbf{y}^{(t)}=Ax^{(t)}\) and \(\|\mathbf{y}-\mathbf{y}^{(t)}\|^{2}\) across iterations, noting that an update to a coordinate of \(\mathbf{x}^{(t)}\) requires updating only \(\
u\) coordinates of \(\mathbf{y}^{(t)}\), where \(\
u\) is the degree of a variable in \(A\). Thus, the computational complexity of Algorithm 1 is \(O(T\
u)\), where a typical choice of \(T\) should be \(T=O(M\log M)\) (see Lemma 2 below). ```
Input:\(T=\text{number of steps to run};\,\mathbf{x}^{(0)}\in\{0,1\}^{M}=\text{initial state};\,\mathbf{y}\in\mathbb{R}^{n}=\text{measurements};\) parameters \(\sigma^{2}>0\), \(\lambda\in\mathbb{R}\). for\(t=1,\ldots,T\)do\(i_{t}\sim\text{Uniform}(\{1,\ldots,M\})\); new coordinate to update \(q_{0}^{(t)}=\exp\Bigl{\{}-\frac{1}{2\sigma^{2}}\|\mathbf{y}-A\mathbf{x}_{i_{t}=0}^{(t) }\|^{2}\Bigr{\}}\), \(q_{1}^{(t)}=\exp\Bigl{\{}-\frac{1}{2\sigma^{2}}\|\mathbf{y}-A\mathbf{x}_{i_{t}=1}^{(t) }\|^{2}+\lambda\Bigr{\}}\) Update: \(\mathbf{x}_{i_{t}}^{(t+1)}=1\) w.p.