A typical reconstruction limit of compressed sensing based on Lp-norm minimization

Introduction

Compressed (or compressive) sensing is a technique for reconstructing a high dimensional signal from lower dimensional data, the components of which represent partial information about the signal, utilizing prior knowledge on the sparsity of the signal. The research history of this technique is rather long ; but the horizon of the research field is now expanding rapidly after recent publication of a series of influential papers .

We also assume that $F$ is known and that $x$ is sparse in the sense that the number of non-zero elements of $x$ is limited to $\rho N$, where $0\leq\rho\leq 1$. Then, under what conditions can the original signal $x$ be correctly reconstructed from the compressed expression $y$?

It is obvious that eq. (1) in itself cannot determine a unique solution of $x$ because the dimension of $y$, $P$, is smaller than that of $x$, $N$. However, the assumption on the sparsity of $x$ may allow correct reconstruction. In the research on compressed sensing, minimization of a cost function with respect to the $L_{p}$-normEq. (5) does not define a norm in the mathematical sense because it violates the triangle inequality.

subject to the constraints of eq. (1) has been actively studied toward designing efficient reconstruction schemes exploiting such sparsity .

Results of indicate that the following proposition holds. Let us suppose that $x$ is an arbitrary continuous real vector the number of non-zero elements of which is bounded above by $S$. When each entry of $F$ is an independently and identically distributed (i.i.d.) Gaussian random number, the probability of failure in reconstructing $x$ based on the $L_{1}$-norm minimization becomes arbitrarily small as $N$ tends to infinity if the inequalities

hold simultaneously. These inequalities constitute a sufficient condition for arbitrarily reducing the probability of failure for the $L_{1}$-based reconstruction of the arbitrary vector $x$. However, earlier studies on several other problems in information theory indicate that critical conditions of such worst cases are, in general, considerably different from those of typical cases , and are not necessarily relevant in practical situations.

This Letter is written from such a perspective. More precisely, we will herein assess a critical condition for successfully reconstructing $x$ in typical cases in the limit $N,P\to\infty$, but keeping $\alpha=P/N$ finite, utilizing methods of statistical mechanics. Results of numerical experiments reported in indicate that a critical condition of the reconstruction success for typical cases is far from that of eqs. (6) and (7). Our result is in excellent agreement with this indication.

Problem setting

For generality, we formally consider a general reconstruction scheme

utilizing a cost function with respect to the $L_{p}$-norm. We will refer to eq. (8) as $L_{p}$-reconstruction. In the following, we will generally examine the typical reconstruction performance for the cases of $p=0,1$ and $2$ in the limit $N,P\to\infty$, but keeping compression rate $\alpha=P/N$ finite.

In a recent work, utility of the $L_{p}$-norm cost function in estimating \mbox{\boldmath{x}}^{0} from \mbox{\boldmath{y}}+\mbox{\boldmath{n}} is examined, where $n$ is a zero mean Gaussian noise vector . In such problem setting, however, correct reconstruction of \mbox{\boldmath{x}}^{0}, which we will focus on hereinafter, is not possible as long as the variance per element of $n$ is finite.

Analysis

To directly assess the typical performance of the $L_{p}$-reconstruction, we have to solve eq. (8) and examine whether the solution that is obtained is identical to \mbox{\boldmath{x}}^{0} or not for each sample of randomly generated $F$ and \mbox{\boldmath{y}}(=F\mbox{\boldmath{x}}^{0}). Carrying this out analytically is, unfortunately, difficult in practice. To avoid this difficulty, we convert the constrained minimization problem of eq. (8) to a posterior distribution of the inverse temperature $\beta$, thus:

where Z(\beta;\mbox{\boldmath{y}})=\int d\mbox{\boldmath{x}}e^{-\beta||\mbox{\boldmath{x}}||_{p}}\delta\left(F\mbox{\boldmath{x}}-\mbox{\boldmath{y}}\right) plays the role of a partition function. In the limit $\beta\to\infty$, eq. (9) generally converges to a uniform distribution over the solutions of eq. (8). Therefore, one can evaluate the performance of the $L_{p}$-reconstruction scheme by examining the macroscopic behavior of eq. (9) as $\beta\to\infty$, for which one can utilize methods of statistical mechanics.

A distinctive feature of the current problem is that eq. (9) depends on the predetermined (quenched) random variables $F$ and \mbox{\boldmath{x}}^{0} (through \mbox{\boldmath{y}}=F\mbox{\boldmath{x}}^{0}), which naturally leads us to applying the replica method . Under the replica symmetric (RS) ansatz, this yields an expression of the typical free energy density as $\beta\to\infty$ as

where $\left[\cdots\right]$ represents the operation of averaging with respect to $F$ and \mbox{\boldmath{x}}^{0}, and $\mathop{\rm extr}_{X}\{{\cal G}(X)\}$ denotes extremization of a function ${\cal G}(X)$ with respect to $X$, $\Theta=\{Q,\chi,m,\widehat{Q},\widehat{\chi},\widehat{m}\}$, $Dz=dz\exp(-z^{2}/2)/\sqrt{2\pi}$ is a Gaussian measure and

The term $\mathop{\rm min}_{X}\left\{{\cal G}(X)\right\}$ denotes minimization of ${\cal G}(X)$ with respect to $X$. A sketch of the derivation is shown in A.

Three issues are noteworthy here. The first issue concerns the physical meanings of the variables introduced in eq. (12). For example, at the extremum, values of $Q$ and $m$ in eq. (12) correspond to N^{-1}\left[\left\langle|\mbox{\boldmath{x}}\right|^{2}\rangle\right] and N^{-1}\left[\mbox{\boldmath{x}}^{0}\cdot\left\langle\mbox{\boldmath{x}}\right\rangle\right], respectively, where $\left\langle\cdots\right\rangle$ denotes averaging with respect to eq. (9) as $\beta\to\infty$ and |\mbox{\boldmath{a}}| denotes the ordinary Euclidean norm $\sqrt{\sum_{i}|a_{i}|^{2}}$ for a vector \mbox{\boldmath{a}}=(a_{i}). This indicates that the typical value of the mean square error per component {\rm MSE}=N^{-1}\left[\left\langle|\mbox{\boldmath{x}}-\mbox{\boldmath{x}}^{0}|^{2}\right\rangle\right] can be assessed as

holds, which represents the de Almeida-Thouless (AT) instability condition for the present problem . When eq. (16) holds for the extremum solution of eq. (12), the RS treatment is not valid and one has to explore more general solutions taking the effect of replica symmetry breaking (RSB) into account to accurately assess the performance of the $L_{p}$-reconstruction.

Results

For $p=0,1$ and $2$, we numerically solved the RS extremization problem of eq. (12) for various pairs of $\alpha$ and $\rho$. In all cases, only a single stable solution was found. Given $\rho$ and $p$, the solution found for sufficiently large $\alpha$ was always characterized by $Q=m=\rho$ indicating successful reconstruction. However, as $\alpha$ was lowered, the success solution lost its local stability (against the RS disturbance) and a transition to a failure solution of $Q\neq m\neq\rho$ occurred.

For the success solution, conjugate variables $\widehat{Q}$ and $\widehat{m}$ were always infinitely large whereas the remaining variables $\chi$ and $\widehat{\chi}$ did not necessarily diverge. Investigating local stability of the success solution yielded a limit $\alpha_{c}(\rho)$, which represented the possibility of $L_{p}$-reconstruction in typical cases. For each of $p=0,1$ and $2$, this is summarized as follows.

The success solution, for which $\chi=0$ and $\widehat{\chi}\to\infty$, is stable if and only if $\alpha>\rho$, which indicates $\alpha_{c}(\rho)=\rho$. The condition $\alpha>\rho$ is necessary to ensure that eq. (1) has a unique solution, even in the situation that all sites of non-zero elements of $x$ are known. This means that the limit for $p=0$ achieves the best possible performance. However, due to discontinuity in the profile of $x_{0}^{*}(h;\widehat{Q})$ (Fig. 1 (a)), eq. (16) always holds for the success solution, indicating that the current RS analysis is not valid. Therefore, further exploration based on various RSB ansätze is necessary for accurately assessing the reconstruction performance, which is, however, beyond the scope of the present Letter.

2 p=1𝑝1p=1

$\widehat{\chi}$ of the success solution is determined by

where $H(x)=\int_{x}^{\infty}Dt$. Utilizing the solution of this equation, the stability condition of the success solution is expressed as

This indicates that the limit for $p=1$ can be expressed as $\alpha_{c}(\rho)=2(1-\rho)H(\widehat{\chi}^{-1/2})+\rho$. $\alpha_{c}(\rho)$ also corresponds to the criticality of eq. (16) and the RS success solution is locally stable against perturbations that break the replica symmetry as long as eq. (18) holds. Therefore, our RS analysis is valid.

3 p=2𝑝2p=2

The success solution is stable if and only if $\alpha\geq 1$, implying $\alpha_{c}(\rho)=1$. For $\alpha>\alpha_{c}(\rho)=1$, eq. (16) does not hold and the RS analysis is valid. Since $\alpha\geq 1$ makes the constraints of eq. (1) sufficient to reconstruct \mbox{\boldmath{x}}^{0} perfectly, this result means that the $L_{2}$-norm minimization is not capable of reconstructing any compressed expressions.

Plots of the results obtained are shown in Fig. 2 (a). We also depict a curve of the worst case critical condition for the $L_{1}$-reconstruction (inset), which is assessed utilizing eqs. (6) and (7) in the limit $N,P\to\infty$, keeping $\alpha=P/N$ and $\rho=S/N$ finite, for comparison. The $L_{1}$-reconstruction can be carried out in practice by interior point methods , the necessary computational cost of which grows as $O(N^{3})$ in the present large system limit. On the other hand, performing the $L_{0}$-reconstruction is, in general, NP hard, although its potential might be superior to that of the $L_{1}$-reconstruction. Fig. 2 (a), in conjunction with these, implies that the $L_{1}$-based scheme is a practically preferable method which balances computational feasibility and relatively high reconstruction capability. This figure also indicates that discrepancy of the values of critical compression rate is huge between the worst and typical case analyses. This implies that there may be much room for improvement of the worst case assessment although we must keep in mind that the criterion of reconstruction success in the present analysis, which permits reconstruction errors of asymptotically negligible size as $N\to\infty$, is different from that of the worst case analysis, in which no errors are allowed.

To justify our assessment, we performed extensive numerical experiments of the $L_{1}$-reconstruction for $\rho=0.5$, the results of which are summarized in Fig. 2 (b). In an experimental trial, an original signal \mbox{\boldmath{x}}^{0} was randomly generated so as to have exactly $S=\rho N=N/2$ non-zero elements, to which i.i.d. Gaussian random numbers of zero mean and unit variance were assigned. For numerically assessing the criticality, the number of constraints $P$ was lowered from $P=N$ one-by-one until the solution of the $L_{1}$-reconstruction, \widehat{\mbox{\boldmath{x}}}, satisfied the condition of ||\widehat{\mbox{\boldmath{x}}}-\mbox{\boldmath{x}}^{0}||_{1}>10^{-4}, and $P_{c}=P+1$ was recorded when the condition was first satisfied. For searching for \widehat{\mbox{\boldmath{x}}}, we used CVX, a package for specifying and solving convex programs . The trials were carried out $10^{6}$ times for a fixed system size $N$ and the experimental critical rate was defined as $\alpha_{c}(\rho=0.5,N)=\overline{P_{c}}/N$, where $\overline{\cdots}$ denotes the arithmetic average over the trials. Quadratic extrapolation from data for $N=10,12,\ldots,30$ yielded an experimental estimate of the critical ratio $\alpha_{c}(0.5)=\lim_{N\to\infty}\alpha_{c}(0.5,N)\simeq 0.83165$, which is in good accordance with the theoretical value $\alpha_{c}(0.5)=0.83129\ldots$ (Fig. 2 (b)). In , experiments for evaluating the critical density $\rho_{c}$ for $\alpha=0.5$ were performed for relatively large systems of $N=512$ and $1024$. Judging from comparison by eye, plots of the results are also consistent with our theoretical estimate $\rho_{c}(\alpha=0.5)=0.19284\ldots$. These indicate that our assessment is at least capable of explaining the experimental results to a high accuracy although mathematical justification of the replica method, in general, has not yet been established .

Summary and discussion

In summary, we have assessed the typical performance of compressed sensing based on minimization with respect to the $L_{p}$-norm for $p=0,1$ and $2$, utilizing the replica method under the replica symmetric (RS) ansatz. Analysis of the stability condition of a solution which represents successful reconstruction yields a critical relation between the compression rate and the signal density that represents the frequency of non-zero elements in the original signal. We have shown that the RS solution of the $L_{0}$-reconstruction achieves the best possible performance, which is, unfortunately, not stable against perturbations that break the replica symmetry. The $L_{2}$-reconstruction has no capability of compressed sensing. On the other hand, our RS analysis has clarified that the $L_{1}$-based scheme does have a considerably high reconstruction ability. Moreover, it has been recognized that the $L_{1}$-reconstruction can be solved via linear programming with a feasible computational cost. These properties are advantageous from the viewpoint of practical utility.

In this Letter, we have assumed that each entry of the compression matrix $F$ is an i.i.d. random variable with zero mean and a fixed variance. Utilizing a technique offered in , the analysis can be extended to cases in which $F$ is randomly generated so as to be characterized as

where $D$ is a diagonal matrix, whose eigenvalue spectrum asymptotically converges to a fixed distribution, and $O$ is a sample from the uniform distribution of $N\times N$ orthogonal matrices, independent of $D$. However, as long as $P\times P$ matrix $FF^{\rm T}$ is typically of full rank, which is the case when entries of $F$ are i.i.d. random numbers of zero mean and a fixed variance, the result is identical to that obtained here (see B). One can also show that the values of $\alpha_{c}(\rho)$ do not depend on details of the distribution of the non-zero elements of $x$ as long as the mean and variance are finite. This implies that the findings of this Letter generally hold for relatively wide classes of compression matrices and signals.

Performance assessment of the $L_{0}$-reconstruction based on a replica symmetry breaking ansatz and development of mean field algorithms for approximately solving the reconstruction problems with a lower computational cost are currently under way.

Note added – After submitting this Letter, the authors noticed that the typical criticality of the $L_{1}$-reconstruction was explored in for compression matrices consisting of i.i.d. zero mean Gaussian random column vectors utilizing techniques of combinatorial geometry. It turns out that their weak threshold corresponds to our result for $\alpha_{c}(\rho)$ with $p=1$. In view of their criterion of reconstruction success, in which no errors are allowed, our result implies that the criticality of $L_{1}$-reconstruction is “tight” in the sense that it does not change irrespective of whether or not we allow small errors which are vanishing asymptotically as $N\to\infty$. The connection between our analysis and theirs further suggests a possibility of wide application of statistical-mechanics tools to problems in large-dimensional random combinatorial geometry.

Appendix A Derivation of eq. (12)

play a key role in deriving eq. (12), where \mbox{\boldmath{y}}=F\mbox{\boldmath{x}}^{0} is used and $u_{\mu}^{a}=\sum_{i=1}^{N}F_{\mu i}x_{i}^{a}$ $(\mu=1,2,\ldots,P;\;a=0,1,2,\ldots,n)$. When $F_{\mu i}$ are i.i.d. Gaussian random variables of mean zero and variance $1/N$, which is mainly assumed in this Letter, the central limit theorem guarantees that $u_{\mu}^{a}$ can be handled as zero mean multivariate Gaussian random variables which are characterized by the covariances $\left[u_{\mu}^{a}u_{\nu}^{b}\right]_{F}=Q_{ab}\delta_{\mu\nu}$ for a fixed set of \mbox{\boldmath{x}}^{0},\mbox{\boldmath{x}}^{1},\mbox{\boldmath{x}}^{2},\ldots,\mbox{\boldmath{x}}^{n}, where $\left[\cdots\right]_{F}$ denotes the operation of averaging with respect to $F$, and Q_{ab}=Q_{ba}=N^{-1}\mbox{\boldmath{x}}^{a}\cdot\mbox{\boldmath{x}}^{b}. $\delta_{\mu\nu}$ is unity for $\mu=\nu$ and vanishes otherwise. Under the RS ansatz

this indicates that $u_{\mu}^{a}$ can be expressed as $u_{\mu}^{a}=\sqrt{Q-q}s_{\mu}^{a}+\sqrt{q}t_{\mu}$ $(a=1,2,\ldots,n)$ and $u_{\mu}^{0}=\sqrt{\rho-m^{2}/q}s_{\mu}^{0}+m/\sqrt{q}t_{\mu}$, where $s_{\mu}^{a}$ and $t_{\mu}$ are i.i.d. Gaussian random variables of zero mean and unit variance. Employing these expressions to eq. (21) yields

On the other hand, the saddle point method offers an expression for the volume of the subshell corresponding to the RS order parameters (26) as

Appendix B Treatment of rotationally-invariant matrix ensembles

For $\alpha=P/N\leq 1$, let us suppose that compression matrix $F$ is characterized as eq. (19), where eigenvalues of diagonal matrix $D$ asymptotically follow a fixed distribution $r(\lambda)=(1-\alpha)\delta(\lambda)+\alpha\widetilde{r}(\lambda)$ as $N,P\to\infty$ with keeping $\alpha=P/N\sim O(1)$ and $O$ is sampled from the uniform distribution of $N\times N$ orthogonal matrices. This matrix ensemble is invariant under any rotation of coordinates. For simplicity, we assume that the support of $\widetilde{r}(\lambda)$ is a certain finite range away from the origin $\lambda=0$, which implies that $P$ row vectors of a typical sample of $F$ are linearly independent. Under the RS ansatz, employment of a formula developed in to eq. (20) offers an expression

for $x>0$. For $x\gg 1$, extremization in eq. (34) provides $\Lambda\simeq(1-\alpha)/x$, which leads to an asymptotic form $G(-x)\simeq-(\alpha/2)\ln x-(\alpha/2)(1+\int d\lambda\widetilde{r}(\lambda)\ln\lambda)-((1-\alpha)/2)\ln(1-\alpha)$. Utilizing this expression in assessment of eq. (33), where contributions of $O(\ln\tau)$ which arise from the $G$-functions are canceled with the term of $-(n\alpha/2)\ln\tau$ avoiding divergence as $\tau\to+0$, indicates that the difference between eqs. (33) and (28) is only a constant independently of $\beta$. This means that in the vanishing temperature limit as $\beta\to\infty$ free energy for the rotationally-invariant ensemble exactly accords with eq. (12). Therefore, the result is identical to that obtained in the main part of this Letter.

References

References