Statistical-Physics-Based Reconstruction in Compressed Sensing

Compressed sensing has triggered a major evolution in signal acquisition. It consists of sampling a sparse signal at low rate and later using computational power for the exact reconstruction of the signal, so that only the necessary information is measured. Current reconstruction techniques are limited, however, to acquisition rates larger than the true density of the signal. We design a new procedure that is able to reconstruct the signal exactly with a number of measurements that approaches the theoretical limit, i.e., the number of nonzero components of the signal, in the limit of large systems. The design is based on the joint use of three essential ingredients: a probabilistic approach to signal reconstruction, a message-passing algorithm adapted from belief propagation, and a careful design of the measurement matrix inspired by the theory of crystal nucleation. The performance of this new algorithm is analyzed by statistical-physics methods. The obtained improvement is confirmed by numerical studies of several cases.

Popular Summary

When you want an image with $100\times 100$-pixel resolution, how many imaging data points do you need? $10000$, the obvious answer seems to be. No, says compressed sensing, a branch of signal-processing science and computational statistics. This research branch develops methods of acquiring and reconstructing a signal that rely on fewer data-acquisition points or measurements than what was previously considered possible, and has enabled faster and more precise measurements in a wide range of applications, including cameras, computed tomography, magnetic-resonance imaging, and genome sequencing. Current techniques are, however, still not optimal, requiring more data points or measurements than necessary. In this interdisciplinary paper, we have designed and tested a new compressed-sensing strategy that allows us to go significantly beyond the known limits for data acquisition and signal reconstruction, by bringing methods and inspirations from statistical physics to bear on compressed sensing—a field that traditionally does not belong to physics.

Mathematically, the problem of compressed sensing is easy to formulate. The signal to be reconstructed is an $N$-component one, represented by a vector$\mathbf{s}$ ($N$ can be viewed as the total number of pixels desired), and compressed sensing deals with signals that are “sparse,” with only a fraction (${\rho }_0$) of the $N$ components being nonzero. The actual acquired data is represented by an $M$-component vector, $\mathbf{y}$, with $M$ also being only a fraction $\alpha$ of $N$. $\mathbf{y}$ and $\mathbf{s}$ are connected linearly by a known mapping. Reconstruction of signal $\mathbf{s}$ then becomes a problem of performing the inverse mapping. Current strategies use linear-programming-based optimization to determine the right inverse mapping. But, they require more data points than would be absolutely necessary, i.e., $\alpha$ must be larger than ${\rho }_0$.

The strategy we have designed, called the seeded belief-propagation algorithm, approaches the signal acquisition and reconstruction in a radically different way. It views the ultimate signal to be reconstructed as the result of a sufficient sampling of a probabilistic Boltzmann-measure-like distribution, which depends on the mapping that corresponds to data-acquisition measurements and the acquired data points. Using a class of specially designed mappings that exploit the physics intuition about crystal nucleation and a sampling technique known to statistical physics as the belief-propagation algorithm, our new compressed-sensing protocols achieve, in a computationally efficient manner, exact reconstruction of the original signal, not only when $\alpha$ is larger than ${\rho }_0$, but also as $\alpha$ approaches ${\rho }_0$, in other words, as the number of acquired data points approaches the absolute minimum. In fact, the examples we have tested show that the gains are stunning.

We would like to think that our interdisciplinary approach breaks new ground in the field of compressed sensing, and we anticipate many further fundamental developments and practical applications.

Figure 1

Two illustrative examples of compressed sensing in image processing. Top: The original image, known as the Shepp-Logan phantom, of size $N={128}^2$, is transformed via one step of Haar wavelets into a signal of density ${\rho }_0\approx 0.15$. With compressed sensing, one is thus in principle able to reconstruct the image exactly with $M\ge {\rho }_{0}N$ measurements, but practical reconstruction algorithms generally need $M$ larger than ${\rho }_{0}N$. The five columns show the reconstructed figure obtained from $M=\alpha N$ measurements, with decreasing acquisition rate $\alpha$. The first row is obtained with the ${\ell }_1$-reconstruction algorithm [3, 4] for a measurement matrix with independent identically distributed (iid) elements of zero mean and variance $1/N$. The second row is obtained with belief propagation, using exactly the same measurements as used in the ${\ell }_1$ reconstruction. The third row is the result of the seeded belief propagation, introduced in this work, which uses a measurement matrix based on a chain of coupled blocks (see Fig. 4). The running time of all three algorithms is comparable. (Asymptotically, they are all quadratic in the limit of the large signal size.) In the lower set of images, we took as the sparse signal the relevant coefficients after two-step Haar transform of the picture of Lena. Again for this signal, with density ${\rho }_0=0.24$, the s-BP procedure reconstructs exactly down to very low measurement rates. (Details are in Appendix ; the data is available online [6].)

Figure 2

Phase diagrams for compressed-sensing reconstruction for two different signal distributions. On the left, the ${\rho }_{0}N$ nonzero components of the signal are independent Gaussian random variables with zero mean and unit variance. On the right, they are independent $\pm 1$ variables. The measurement rate is $\alpha =M/N$. On both sides, we show, from top to bottom: (a) The phase transition ${\alpha }_{{\ell }_1}$ for ${\ell }_1$ reconstruction [5, 11, 12] (which does not depend on the signal distribution). (b) The phase transition ${\alpha }_{\mathrm{EM}\mathrm{\text{−}}\mathrm{BP}}$ for EM-BP reconstruction based for both sides on the probabilistic model with Gaussian $\varphi$. (c) The data points which are numerical-reconstruction thresholds obtained with the s-BP procedure with $L=20$. The point gives the value of $\alpha$ where exact reconstruction was obtained in 50% of the tested samples; the top of the error bar corresponds to a success rate of 90%; the bottom of the bar to a success of 10%. The shrinking of the error bar with increasing $N$ gives numerical support to the existence of the phase transition that we have studied analytically. These empirical reconstruction thresholds of s-BP are quite close to the $\alpha ={\rho }_0$ optimal line, and they move closer to it with increasing $N$. The parameters used in these numerical experiments are detailed in Appendix . (d) The line $\alpha ={\rho }_0$ that is the theoretical reconstruction limit for signals with continuous ${\varphi }_0$. An alternative presentation of the same data using the convention of Donoho and Tanner [5] is shown in Appendix .

Figure 3

When sampling from the probability $\widehat{P}(\mathbf{x})$, in the limit of large $N$, the probability that the reconstructed signal $x$ is at a squared distance $D=\sum _i(x_i-s_i{)}^2/N$ from the original signal $s$ is written as $c{e}^{N\Phi (D)}$, where $c$ is a constant and $\Phi (D)$ is the free entropy. Left: $\Phi (D)$ for a Gauss-Bernoulli signal with ${\rho }_0=\rho =0.4$ in the case of an unstructured measurement matrix $\mathbf{F}$ with independent random Gaussian-distributed elements. The evolution of the EM-BP algorithm is basically a steepest ascent in $\Phi (D)$ starting from large values of $D$. It goes to the correct maximum at $D=0$ for large values of $\alpha$, but it is blocked in the local maximum that appears for $\alpha <{\alpha }_{\mathrm{EM}\mathrm{\text{−}}\mathrm{BP}}({\rho }_0=0.4)\approx 0.59$. For $\alpha <{\rho }_0$, the maximum is not at $D=0$, and exact inference is impossible. The seeding matrix $\mathbf{F}$, leading to the s-BP algorithm, succeeds in eliminating this local maximum. Right: Convergence time of the EM-BP and s-BP algorithms obtained through the replica analysis for ${\rho }_0=0.4$. The EM-BP convergence time diverges as $\alpha \to {\alpha }_{\mathrm{EM}\mathrm{\text{−}}\mathrm{BP}}$ with the standard $L=1$ matrices. The s-BP strategy allows one to go beyond the threshold: Using ${\alpha }_1=0.7$ and increasing the structure of the seeding matrix (here, $L=2$, 5, 10, 20), we approach the limit $\alpha ={\rho }_0$. (Details of the parameters are given in Appendix .)

Figure 4

Construction of the measurement matrix $\mathbf{F}$ for seeded compressed sensing. The elements of the signal vector are split into $L$ (here, $L=8$) equal-sized blocks; the number of measurements in each block is $M_p={\alpha }_{p}N/L$ (here, ${\alpha }_1=1$, ${\alpha }_p=0.5$ for $p=2,\dots ,8$). The matrix elements $F_{\mu i}$ are chosen as random Gaussian variables with variance $J_{q,p}/N$ if variable $i$ is in the block $p$ and measurement $\mu$ is in the block $q$. In the s-BP algorithm, we use $J_{p,q}=0$, except for $J_{p,p}=1$, $J_{p,p-1}=J_1$, and $J_{p-1,p}=J_2$. Good performance is typically obtained with relatively large $J_1$ and small $J_2$.

Figure 5

Evolution of the mean-squared error at different times as a function of the block index for the s-BP algorithm. Exact reconstruction first appears in the left block whose rate ${\alpha }_1>{\alpha }_{\mathrm{EM}\mathrm{\text{−}}\mathrm{BP}}$ allows for seeded nucleation. It then propagates gradually block by block, driven by the free-entropy difference between the metastable and equilibrium states. After little more than 300 iterations, the whole signal of density ${\rho }_0=0.4$ and size $N=50000$ is exactly reconstructed well inside the zone forbidden for BP (see Fig. 2). Here, we used $L=20$, $J_1=20$, $J_2=0.2$, ${\alpha }_1=1.0$, $\alpha =0.5$.

Figure 6

Convergence time of the s-BP algorithm with $L=2$ as a function of $J_1$ and $J_2$ (in $\log $ scale) for Gauss-Bernoulli signal with ${\rho }_0=0.1$. The color represents the number of iterations such that the MSE is smaller than ${10}^{-8}$. The white region gives the fastest convergence. The measurement density is fixed to $\alpha =0.25$. The parameters in s-BP are chosen as $L=2$ and ${\alpha }_1=0.3$. The axes are in ${log}_{10}$ scale.

Figure 7

Same data as in Fig. 2, but using the convention of 5. The phase diagrams are now plotted as a function of the undersampling ratio ${\rho }_{\mathrm{DT}}=K/M=\rho /\alpha$ and of the oversampling ratio ${\delta }_{\mathrm{DT}}=M/N=\alpha$.

Figure 8

Mean-squared error as a function of $\alpha$ for different values of the measurement noise. Left: Gauss-Bernoulli signal with ${\rho }_0=0.2$, and measurement noise with standard deviation $\sqrt{\Delta }={10}^{-4}$. The EM-BP ($L=1$) and s-BP strategy ($L=4$, $J_1=20$, $J=0.1$, ${\alpha }_1=0.4$) are able to perform a very good reconstruction up to much lower value of $\alpha$ than the ${\ell }_1$ procedure. Below a critical value of $\alpha$, these algorithms show a first-order phase transition to a regime with much larger MSE. Right: Gauss-Bernoulli signal with ${\rho }_0=0.4$, with a noise with standard deviation $\sqrt{\Delta }={10}^{-3}$, ${10}^{-4}$, ${10}^{-5}$. s-BP, with $L=9$, $J_1=30$, $J_2=8$, ${\alpha }_1=0.8$ decodes very well for all of these values of noises. In this case, ${\ell }_1$ is unable to reconstruct for all $\alpha <0.75$, well outside the range of this plot.