Recovery of binary sparse signals from compressed linear measurements via polynomial optimization

TL;DR

This paper introduces a novel polynomial optimization approach for recovering binary sparse signals from limited linear measurements, offering a new perspective and comparison to existing methods in binary compressed sensing.

Contribution

It presents a new polynomial optimization formulation specifically for binary sparse signal recovery and analyzes its performance relative to current state-of-the-art methods.

Findings

The polynomial optimization approach effectively recovers binary sparse signals.

Compared to existing methods, it offers competitive or improved recovery performance.

The approach provides a new mathematical framework for binary compressed sensing.

Abstract

The recovery of signals with finite-valued components from few linear measurements is a problem with widespread applications and interesting mathematical characteristics. In the compressed sensing framework, tailored methods have been recently proposed to deal with the case of finite-valued sparse signals. In this work, we focus on binary sparse signals and we propose a novel formulation, based on polynomial optimization. This approach is analyzed and compared to the state-of-the-art binary compressed sensing methods.

Figures

Equations

$$y=Ax,\text{leavevmode}\text{leavevmode}\text{leavevmode}x\in \{0,1{\}}^n,\text{leavevmode}\text{leavevmode}y\in \text{mathds}{R}^m,\text{leavevmode}A\in \text{mathds}{R}^{m,n},\text{leavevmode}\text{leavevmode}m<n$$

$$\begin{array}{rl}& \underset{x\in [0,1{]}^n}{\min }\sum _{i=1}^n\left(x_i-x_i^2\right)\text{leavevmode}\text{leavevmode}\text{ s. t. }\text{leavevmode}y=Ax.\end{array}$$

$$\begin{array}{rl}\underset{x\in [0,1{]}^n,z\in \text{mathds}{R}^{m,n-1}}{\min }\sum _{i=1}^n\left(x_i-x_i^2\right)& \text{ s. t., for each }i=1\dots ,m,\\ y_i-A_{i,1}{x}_1+z_{i,1}& =0\\ -A_{i,2}{x}_2-z_{i,1}+z_{i,2}& =0\\ & ⋮\\ -A_{i,n}{x}_n-z_{i,n-1}& =0.\end{array}$$

$$\begin{array}{rl}& \underset{x\in [0,1{]}^n}{\min }\sum _{i=1}^n\left(x_i-x_i^2\right)\text{ s. t. }\Vert Ax-y{\Vert }_{\infty }\le \eta .\end{array}$$

$$\begin{array}{rl}\underset{x\in [0,1{]}^n,z,w\in \text{mathds}{R}^{m,n-1}}{\min }\sum _{i=1}^n\left(x_i-x_i^2\right)& \text{ s. t. for each }i=1\dots ,m,\\ \eta +y_i-A_{i,1}{x}_1+z_{i,1}\ge 0,\text{leavevmode}\text{leavevmode}\text{leavevmode}& \eta -y_i+A_{i,1}{x}_1+w_{i,1}\ge 0\\ -A_{i,2}{x}_2-z_{i,1}+z_{i,2}\ge 0,\text{leavevmode}\text{leavevmode}\text{leavevmode}& A_{i,2}{x}_2-w_{i,1}+w_{i,2}\ge 0\\ & ⋮\\ -A_{i,n}{x}_n-z_{i,n-1}\ge 0,\text{leavevmode}\text{leavevmode}\text{leavevmode}& A_{i,n}{x}_n-w_{i,n-1}\ge 0\end{array}$$

Full text

1 Introduction

Solving underdetermined linear systems with finite-valued solutions has recently attracted a lot of attention in signal processing. This mathematical problem is encountered whenever a digital signal has to be recovered from compressed or incomplete linear measurements, which occurs in a number of applications, ranging from digital communications and digital image processing, see, e.g., [8, 31, 29, 5, 1, 16, 17], to localization, source separation, and spectrum sensing, see, e.g., [2, 22, 18, 23, 11, 14].

Formally, the problem is defined as follows: given a finite alphabet $\mathcal{A}$ (that is, a finite set of symbols), one aims to recover $\widetilde{x}\in\mathcal{A}^{n}$ from compressed measurements $y=A\widetilde{x}$ (possibly corrupted by noise), where $A\in\mathds{R}^{m,n}$, $m<n$, is a known matrix. The solution is assumed to be unique within $\mathcal{A}^{n}$. The problem is intrinsically combinatorial: if $c$ is the cardinality of $\mathcal{A}$, the solution could be found by listing all the $c^{n}$ possible vectors, but clearly this is unfeasible for large $n$.

If $\widetilde{x}$ is sparse, i.e., it has many null components, sparsity-promoting procedures can be applied, in particular the ones developed in compressed sensing (CS, [7, 15]). We recall that the CS paradigm states that (real-valued) sparse vectors can be recovered from few linear measurements, under suitable conditions on the matrix $A$ (for example, random matrices are fitting). Standard sparsity-promoting and CS techniques are conceived for real-valued signals, and do not envisage a possible prior information on the discrete nature of the signal. To fill this gap, in the last years, novel CS strategies have been proposed, which exploit the knowledge of $\mathcal{A}$. These strategies are shown to significantly improve the performance with respect to algorithms unaware of the discrete structure, see, e.g., [18, 11, 14] and the references therein.

In the literature on the recovery of finite-valued sparse signals from linear compressed measurements, a lot of work is devoted to the binary case $\mathcal{A}=\{0,1\}$. This line of research is also known as binary compressed sensing (BCS, [25, 28, 13]). CS strategies tailored for BCS generally require less measurements than classical CS techniques. Different approaches are proposed, e.g., Bayesian models in [31], bi-partite graph models in [25], $\ell_{0}$ minimization in [23], $\ell_{\infty}$ minimization in [1], $\ell_{1}$ minimization in [29, 2], greedy methods in [11], and non-convex models in [13].

The interest on binary sparse signals is motivated by several applications. For example, in the framework of indoor localization, this setting is frequently considered [3]: an area is subdivided into $n$ cells, and $m$ sensors are deployed to receive signals from targets on the area. In the training phase, a target transmits in turn from each cell, in order to build a dictionary $A\in\mathds{R}^{m,n}$. Each entry $A_{i,j}$ represents the received signal strength measured by sensor $i$ when the target is in the cell $j$. In the running phase, the cells which contain transmitting targets are detected as solution of $y=Ax$, where $x\in\{0,1\}^{n}$. $A$ has then a mixing effect on signals transmitted by targets. Since usually the number of targets is much smaller than the number of cells, $x$ is sparse and CS can be applied. Beyond localization, the literature on the recovery of binary vectors from linear measurements is extensive in digital data transmission, jump systems, source separation, fault detection, see, e.g., [30, 12, 9, 4] and references therein. Interestingly, we notice that the problem $Ax=y$ with $x\in\{0,1\}^{n}$ is also a variant of the popular 0-1 knapsack problem [10]. Each equation of the system is a 0-1 knapsack problem whose goal is to fill a knapsack of capacity $y_{i}$ choosing among $n$ items with weights $A_{i,1},\dots,A_{i,n}$.

In this work, we propose a novel polynomial optimization approach to recover binary sparse signals from compressed linear measurements. Specifically, we show that the problem can be formulated as a non-convex polynomial optimization problem (POP, [21]), whose solution can be estimated efficiently via semidefinite programming. Our contribution can be summarized as follows. First, we propose a non-convex POP whose global minimum is exactly the desired solution, in the noise-free case. Second, we show that this POP can be solved via semidefinite programming; moreover, it has a sparse structure as described in [20], which reduces the computational complexity. Finally, we consider the noisy case, and we perform numerical tests to illustrate the recovery performance. In particular, we compare the proposed approach to the state-of-the-art BCS methods proposed in [11, 13].

The paper is organized as follows. In Section 2, we introduce and analyze the proposed optimization problem. In Section 3, we extend the model to the noisy setting. In Section 4 is devoted to numerical simulations. Finally, we draw some conclusions.

2 Problem statement and analysis

Let us consider the problem:

[TABLE]

where $A$ and $y$ are known. We assume that (1) is well-posed.

Assumption 1**.**

The solution $\widetilde{x}\in\mathds{R}^{n}$ of problem (1) is unique.

This assumption holds for many classes of $A$, for example, when the entries of $A$ are in general position (see [33] for the definition), which is guaranteed when the entries of $A$ are generated at random according to any continuous distribution.

The key idea of this work is to tackle problem (1) by reformulating it as follows:

[TABLE]

The optimization problem in (2) is well-posed, as proven in the following proposition.

Proposition 1**.**

Under Assumption 1, the unique solution of problem (2) is $\widetilde{x}$, for any $m\geq 1$.

Proof.

The global minima of $f(x):=\sum_{i=1}^{n}\left(x_{i}-x_{i}^{2}\right)$ over $[0,1]^{n}$ are in $\{0,1\}^{n}$. In particular, $f(x)\geq 0$ for any $x\in[0,1]^{n}$, and $f(x)=0$ only when $x\in\{0,1\}$. This is consistent with the constraint $y=Ax$, which has a unique binary solution $\widetilde{x}$. Therefore, $\widetilde{x}$ is the global minimum of problem (2). ∎

Problem (2) is non-convex, therefore its solution is not straightforward. However, since it is a POP, results proven in [19, 26, 21] can be exploited to the purpose. In a nutshell, these results state that the solution of a POP can be achieved by solving a sequence of relaxed semidefinite programming (SDP) problems, whose global minima converge to the global minimum of the POP. Thus, an estimate of the POP solution can be obtained by solving an SDP problem of sufficiently large relaxation order. On the other hand, the dimensions of the SDP problems increase with the relaxation order, therefore it is important to keep the order as low as possible. Interestingly, for Problem (2), we can prove that relaxation order 1 is sufficient to achieve the exact global minimum.

Proposition 2**.**

The global minimum of the SDP relaxation of order 1 of Problem (2) corresponds to the global minimum of Problem (2). The global minimizer of Problem (2) can be extracted if the solution of the SDP relaxation has rank 1.

Proof.

The first part follows from Theorem 4.2 in [19], which, among others, states the following result. Let $p(x):\mathds{R}^{n}\to\mathds{R}$ be a polynomial and $K$ be the compact set described by $\{g_{i}(x)\geq 0,i=1,\dots,m\}$, where $g_{i}$’s are polynomials. Let $p^{\star}=\min_{x\in K}p(x)$. If $p(x)-p^{\star}=q_{0}+\sum_{i=1}^{m}g_{i}(x)q_{i}$ where $q_{j}$, $j=0,\dots,m$, are non-negative constants, then the solution of the relaxed SDP of order $1$ is equal to $p^{\star}$.

For Problem (2), we know that $p^{\star}=0$, and we can replace $x_{i}\in[0,1]$ with $x_{i}-x_{i}^{2}\geq 0$, thus the condition of [19, Theorem 4.2] holds with $q_{0}=0$, $q_{i}=1$ for $g_{i}(x)=x_{i}-x_{i}^{2}$, and $q_{i}=0$ for the constraints $Ax=b$. Finally, the possibility of extracting the minimizer when the rank is 1 is a direct application of [21, Theorem 6.6]. ∎

In conclusion, we can reformulate Problem 2 as a single SDP problem. More precisely, this SDP problem corresponds to the SDP relaxation for quadratic problems, which is written explicitly, e.g., in [34, Section 4.5].

Remark 1**.**

Problem (2) and the related analysis can be extended to any binary alphabet $\{\alpha,\beta\}$ by using $(x-\alpha)(x-\beta)$ as cost functional. Similarly, extensions to any finite alphabet $\{\alpha_{1},\dots,\alpha_{c}\}$, with $\alpha_{1}<\dots<\alpha_{c}$ are possible by using the cost functional $\sum_{i=1}^{c-1}[(x-\alpha_{i})(x-\alpha_{i+1})]_{+}$, where $[\cdot]_{+}$ denotes the positive part. In this last case, the problem is semi-algebraic, and the theory in [21] still applies. Such extensions will be object of future work.

2.1 Complexity reduction via chordal sparsity

Proposition 2 recasts the non-convex problem (2) into SDP, thus its solution can be computed via convex optimization. A drawback of this approach might be the dimension of the SDP problem, which increases with $n$. However, the chordal sparsity property reduces the complexity, as described in [20, 34, 36]. Briefly, chordal sparsity is defined as follows: if the matrices involved in an SDP problem are sparse, the sparsity pattern can be represented by an undirected graph, see, e.g., [36]. If this graph is chordal, i.e., each cycle composed by at least four nodes has a chord, then the SDP problem can be decomposed into smaller sub-problems (see Theorems 1 and 2 in [36]), which reduces the complexity.

As illustrated in [20], chordal sparsity can be evaluated directly on the original POP as follows. Given a constrained POP with $n$ variables, such as (2), let us build $p$ sets $I_{l}\subset\{1,\dots,n\}$, $l=1,\dots,p$, with the following properties: their union is $\{1,\dots,n\}$; the variables $x_{i}$ contained in each constraint equation are concerned with a single $I_{l}$; the objective function is sum of monomials such that each monomial is concerned with a single $I_{l}$. We then define the running intersection property as in [20]: for every $l=1,\dots,p-1$, $I_{l+1}\cap\left(\cup_{j=1}^{l}I_{j}\right)\subseteq I_{s}$ for some $s\leq l$. If the considered POP satisfies the running intersection property, then the graph associated with its SDP relaxation is chordal, and Theorems 1 and 2 in [36] hold. In particular, the $I_{l}$’s represent the maximal cliques of such graph.

Problem (2) is sparse, as there are no mixed products between variables; this implies that the matrices in the SDP relaxation are null except for the main diagonal, the first row, and the first column, see [34, Section 4.5] for details. Moreover, Problem (2) fulfills the running intersection property if we rewrite the constraints $Ax=y$ by adding suitable slack variables $z_{i,j}$, $i=1\dots,m$, $j=1\dots,n-1$. Specifically,

[TABLE]

The following result holds (the proof is straightforward and omitted for brevity).

Proposition 3**.**

Problem (3) satisfies the running intersection property with the following sets: $I_{1}=\{x_{1},z_{1,1},\dots,z_{m,1}\};$ $I_{l}=\{x_{l},z_{1,l-1},z_{1,l},\dots,z_{m,l-1},z_{m,l}\}$ for $l=2,\dots,n-1$; $I_{n}=\{x_{n},z_{1,n-1},\dots,z_{m,n-1}\}$.

According to [20, 36], the running intersection property allows us to recast our SDP problem of dimension $n^{2}$ into $p=n$ sub-problems of dimension $m^{2}$, which is favorable when $m\ll n$.

2.2 Relation to prior literature

We review the most important BCS approaches to highlight the differences with respect to the proposed Problem (2). In [29], the authors propose to tackle BCS by solving a basis pursuit restricted to $[0,1]^{n}$, namely, $\min_{x\in[0,1]^{n}}f(x)=\sum_{i=1}^{n}x_{i}$ such that $y=Ax$, $A\in\mathds{R}^{m,n}$, $m<n$. This model, later revisited in [2, 18], is shown to require less measurements than classical basis pursuit, and, in particular, when $m>n/2$, it provides the exact solution with high probability for specific random $A$, see [29, 2]. A similar result is obtained by using linear programming in [24]. In [16, 17], the authors propose a sum-of-absolute-value optimization by assuming to know the probability of $x_{i}=0$. In [13], a non-convex model is proposed: $\min_{x\in[0,1]^{n}}\frac{1}{2}\left\|y-Ax\right\|_{2}^{2}+\lambda\sum_{i=1}^{n}\left(x_{i}-\frac{x_{i}^{2}}{2}\right)$, $\lambda>0$, which can be interpreted as a Lasso [32] with concave (non-decreasing) penalization instead of $\ell_{1}$ penalization to compel sparsity. The global minimum is exactly the desired signal, under weak conditions, and $\ell_{1}$ iterative reweighting algorithms are shown to achieve it, if $m$ is sufficiently large. In [11], the PROMP algorithm, based on orthogonal matching pursuit, is proposed, which requires a larger $m$ than [13], while it is faster.

3 Extension to the noisy case

Problem (2) can be extended to the case of noisy measurements, namely $y=A\widetilde{x}+\zeta$, $\zeta\in\mathds{R}^{n}$ being an unknown disturbance. If we assume that $\|\zeta\|_{\infty}\leq\eta$, where $\eta>0$ is known, Problem (2) can be reformulated as follows:

[TABLE]

As $\|Ax-y\|_{\infty}\leq\eta$ corresponds to $|A_{i}x-y_{i}|\leq\eta$ for each $i=1,\dots,n$, we can write the sparse formulation as follows:

[TABLE]

where $z_{i,j}$ and $w_{i,j}$, $i=1\dots,m$, $j=1\dots,n-1$, are slack variables. The number of slack variables is doubled with respect to the noise-free case, which increases the whole dimension of the problem. However, since the constraints are all linear, the increase of complexity is not dramatic. Moreover, if $\eta$ is sufficiently small, we can still expect a unique solution under the same conditions of the noise-free case.

4 Numerical results

In this section, we propose numerical simulations that support the effectiveness of the proposed approach. We consider the following setting: the unknown binary signal $\widetilde{x}\in\{0,1\}^{n}$ has $n=100$ components, among which $k\in[10,50]$ components are equal to 1, with sparsity ratio $\frac{k}{n}\in\left[\frac{1}{10},\frac{1}{2}\right]$. We test $m\in[20,50]$. The support of the signal is generated uniformly at random, and the entries of $A$ are Gaussian $\mathcal{N}(0,\frac{1}{m})$. To recover $\widetilde{x}$, we solve the POP’s (2) and (4) by solving the corresponding SDP relaxations with the solver CDCS, see [36, 35], on Matlab R2018a. CDCS solves SDP problems by implementing the alternating direction method of multipliers (ADMM, [6]), which is faster than, e.g., interior point methods. Moreover, CDCS exploits the chordal sparsity of SDP to efficiently reduce the dimensionality of the problem, which perfectly matches with problems (3) and (5).

Our approach is compared to the state-of-the-art BCS algorithms RWR [13] and PROMP [11]. RWR is based on repeating a local minimization iterative algorithm starting from different random initial points, and is shown to outperform basis pursuit based methods, see [13, Section IV]. The following recovery accuracy metrics are considered (we denote by $\widehat{x}$ the estimate of $\widetilde{x}$): exact recovery rate, i.e., how many times the whole signal is exactly recovered; false positive rate, i.e., the number of occurrences $\widehat{x}_{i}\neq 0$ when $\widetilde{x}_{i}=0$ divided by the true number of zeros $n-k$; false negative rate, i.e., the number of occurrences $\widehat{x}_{i}=0$ when $\widetilde{x}_{i}\neq 0$ divided by the true number of ones $k$. The run time is illustrated as well. The results are averaged over 100 runs. Estimations are quantized to $\{0,1\}$ (by imposing a threshold at $\frac{1}{2}$) when algorithms do not provide binary solutions.

4.1 Noise-free experiment

In Figures 1 and 2, we show the results in the noise-free case. In Figure 1, we can see that the proposed POP approach provides the best recovery accuracy, while PROMP is the least accurate. In particular, in contrast to RWR and PROMP, for POP $m=45$ is always sufficient to recover exactly a signal with $k\leq 15$ (the light yellow denotes $100\%$ of success). We remark that the ratio $\frac{m}{k}=3$ is favorable with respect to CS standards [15]. When $k\geq 30$, RWR and PROMP always fail ($0\%$ of success), while POP is still successful in more than $50\%$ of runs with $m=50$. POP and RWR have similar run time, of order $10^{-1}$ seconds. PROMP is faster ($10^{-4}$ seconds); nevertheless, as already said, its performance accuracy is scarce.

From Proposition 1, one might expect to find the exact solution even when $m=1$, using the POP approach. However, the corresponding SDP problem might have several global minima. To obtain the right minimizer the SDP solution must have rank 1. A way to induce low rank could the penalization of the nuclear norm, see [27]; this strategy will be investigated in the future. In this work, we observe that, in practice, a sufficient small $k$ and sufficient large $m$ provide the exact solution.

In Figure 2, we further examine the estimation accuracy by considering the correct recovery of single symbols, in terms of false negative and false positive rate. POP has a more favorable false negative rate with respect to RWR and PROMP, as it does not not induce sparsity. On the other hand, the false positive rate is similar for the three methods, except for large $k$.

4.2 Noisy experiment

We finally revisit the previous experiment by adding measurement noise. We set $\eta=5\times 10^{-2}$, which corresponds to a signal-to-noise ratio $\frac{k}{m\eta^{2}}$ between $25$ and $40$ dB.

In Figures 3 and 4, we can see that, as expected, the general performance is slightly worse when noise is present. However, the mutual behaviors between POP, RWR, and PROMP are similar to those observed in the noise-free case.

5 Conclusion

In this paper, we develop a non-convex quadratic approach to recover binary signals from few linear measurements. The effectiveness of the proposed method leverages on recent results on polynomial optimization, and has better recovery accuracy than binary compressed sensing algorithms, in particular when signals are not very sparse. Future work will envisage the extension of the proposed approach to larger, non-binary, alphabets, and to non-sparse problems, e.g., combinatorial problems, like the 0-1 knapsack problem.