# Feasibility Pump Algorithm for Sparse Representation under Laplacian Noise

Authors: Florin Ilarion Miertoiu; Bogdan Dumitrescu
Journal: Mathematical Problems in Engineering
 (2019)

---

## Abstract

The Feasibility Pump is an effective heuristic method for solving mixed integer optimization programs. In this paper the algorithm is adapted for finding the sparse representation of signals affected by Laplacian noise. Two adaptations of the algorithm, regularized and nonregularized, are proposed, tested, and compared against the regularized least absolute deviation (RLAD) model. The obtained results show that the addition of the regularization factor always improves the algorithm. The regularized version of the algorithm also offers better results than the RLAD model in all cases. The Feasibility Pump recovers the sparse representation with good accuracy while using a very small computation time when compared with other mixed integer methods.


---

## Body

## 1. Problem Formulation

Sparse representations are among the most important linear representation methods. They are directly related to the compressed sensing problem, but can also be used for machine learning, signal or image processing. As presented in [1] different types of algorithms have been proposed. Sparse representations allow the identification of only a few relevant features or atoms, which can be used to represent the signal of interest with high accuracy. They can also help to well separate different types of signals, which makes them appropriate for denoising and classification problems.The sparse representation of a signaly∈Rm, using a dictionary D∈Rm×n, consists in finding a solution x∈Rn having the smallest possible number of nonzero coefficients to the system y=Dx. In practice, the equality constraint is relaxed and replaced by the minimization of the data misfit measure y-Dx. We adopt here the l1 norm, although usually the norm used for the error minimization is l2. The l1 norm is more robust against outliers and is optimal for Laplacian noise [2]. It is also a norm that promotes sparsity (the error between the original signal and the obtained representation has only few significant elements). Bounding the number of nonzero coefficients (denoted x0) by an imposed sparsity level K, the sparse representation is found by solving(1)minimizex∈Rny-Dx1subject  tox0≤KAn alternative formulation of the problem is obtained by minimizing the number of coefficients and requiring the data misfit measure to be under a selected thresholdu:(2)minimizex∈Rnx0subject  toy-Dx1≤uBecause the norm used to bound the number of nonzero coefficients in (1) is l0, Mixed Integer Programming (MIP) algorithms such as Branch and Bound and Feasibility Pump can be used [3].Mixed Integer Programs are optimization problems that contain both integer and continuous variables. The most used MIP algorithms are Branch and Bound [4, 5] and Cutting Planes [6]. Another popular variant is the Branch and Cut [7] algorithm, which combines both algorithms. These algorithms can offer good solutions for problems with integer constraints, but often they are very time consuming. The Feasibility Pump, initially proposed in [8, 9], is a MIP heuristic method that generates two sequences of solutions, one satisfying the linear constraints and one satisfying the integer constraint. These sequences are generated until a solution that satisfies both conditions is found. The Feasibility Pump starts from an initial solution and then proceeds through several iterations to minimize the distance between the two solutions, first by solving a linear optimization problem with the integer constraints removed in order to obtain the solution that satisfies the linear constraint, then by a rounding step in order to obtain the solution that satisfies the integer constraint. The algorithm is prone to cycling and loops and several approaches have been proposed to avoid this [10–14]. The Feasibility Pump offers much better execution time (several orders of magnitude compared to Branch and Bound or Branch and Cut), but the quality of the solution is not as good; there is no guarantee that the optimum is attained. The Feasibility Pump has been proposed for both linear and nonlinear problems [15, 16].In order for the Feasibility Pump algorithm to be used, we introduce a binary variableb∈0,1n, which shows whether an atom of the dictionary is used to represent the signal y.Problem (1) becomes(3)minimizex∈Rn,b∈0,1ny-Dx1subject  to1nTb≤K-Mb≤x≤MbThe big-M trick is used as in [17], where M is a preset parameter bounding the elements of x, chosen as M=1.1DTy∞/D22.In this paper, the Feasibility Pump is adapted in order to obtain the sparse solution of problem (3). A regularization term is also added to the objective in order to improve the convergence of the algorithm. Section 2 presents our Feasibility Pump algorithm and the implementation details. Section 3 is dedicated to experimental results showing the behavior of our algorithm and comparisons with the regularized least absolute deviation (RLAD) model, which consists of a linear programming problem. We show that our regularized Feasibility Pump algorithm is consistently better than its nonregularized version and than RLAD. Section 4 presents the conclusions and future ideas of research.


## 2. Algorithm

Feasibility Pump [8] is a MIP algorithm that tries to minimize the difference between the solution of the relaxed linear programming problem (where the integer condition is relaxed) and a solution that satisfies the integer conditions.Algorithm1 presents our approach for solving problem (1) via its reformulation (3). In what follows, we explain its steps and the refinements that we have made to standard Feasibility Pump.Algorithm 1:Sparse feasibility pump: SFP uses (4) and (5); SFPreg uses (9) and (10).
Data: Dictionary D∈Rm×n, signal to represent y∈Rm, sparsity level K∈N, maximum number of iterations Iter, weight parameters α, γ, λ (the latter for SFPreg only) Result:a feasible solution x∈Rn Solve (4)/(9) for x,b∈Rn using a linear programming solver Use rounding procedure to obtain vectorb~ while number of iterations ≤ Iterdo Solve the problem (5)/(10) for the vectors x and b if b is integer then return x end Use rounding procedure to obtain vector b~ if cycle is detected then Use perturbation (7) on b~ end Update the value of α using (6) endAs usual, the starting point for the Feasibility Pump algorithm is the solution of problem (3) when b is a vector with values in the 0,1n box, instead of a binary one. We also introduce an auxiliary variable w∈Rm to bound the data misfit measure. The resulting problem is(4)minimizex∈Rn,w∈Rm,b∈0,1n1mTwsubject  to1nTb≤K-w≤y-Dx≤w-Mb≤x≤MbThis is a linear programming problem that can be easily solved with one of the many existing algorithms. We have used the MATLAB function linprog, based on the dual simplex algorithm.This realb is then rounded to obtain a binary vector b~, by setting to 1 the K values of b corresponding to the K largest absolute values of x and the others to 0. This is in contrast with the usual rounding to the nearest integer used by Feasibility Pump, but is justified by the fact that a smaller representation error is obtained in the long run by taking the maximum allowed number of nonzero coefficients.In the ideal case where the vectorb is already binary, the algorithm stops. Note that the first constraint of (4) is satisfied; hence the number of coefficients with a value of 1 cannot exceed K. The other stopping condition is when the maximum number of iterations is reached. In this case, b is rounded as described above and the coefficients of x that correspond to zero coefficients in b are also set to zero.In each iteration of Feasibility Pump, the vectorb and the tentative solution x are improved by solving(5)minimizex∈Rn,w∈Rm,b∈0,1n1-α△b,b~+α1mTwsubject  to1nTb≤K-w≤y-Dx≤w-Mb≤x≤Mbwhere △(b,b~)=b-b~1 and α∈(0,1) is a weight. This linear program is also solved with linprog. The iteration step has an objective that combines the representation error with a term that enforces the new solution b to be near from the current integer vector b~, with the aim of making the solution b nearer from a binary vector. This kind of modification is proposed in [18, 19]. After each iteration, the parameter α is reduced and the integer condition will weigh more than the error objective. The reduction of α is done by multiplication with a value γ∈(0,1):(6)α←γα.Typically, a large γ gives the smallest error, but at the expense of a larger execution time, while a smaller γ offers faster results, but with a larger error.Although this strategy usually leads towards a binaryb, it can also lead to cycling, which we consider to occur at iteration k+1 if(i)
the pointbk+1 has already been visited in previous iterations;(ii)
the distance△(bk+1,bk) is under a certain threshold.To break the cycle, a weak perturbationv∈0,1n is used for the coefficients of b~:(7)b~k+1←12b~k+1+v.The resulting b~ has real values and is used in the linear program for the next iteration. This is done to induce a switching of the coefficients in the next vector b~, which will have binary values.We name Sparse Feasibility Pump (SFP) the algorithm as described above.Because the matrixD can be ill conditioned, a second implementation, named Sparse Feasibility Pump with regularization (SFPreg), is proposed, where the objective has an l1 regularization term added. Regularization terms in the Feasibility Pump are also considered in [20]. For the initial step, instead of (3) we aim to solve(8)minimizex∈Rn,b∈0,1ny-Dx1+λx1subject  to1nTb≤K-Mb≤x≤Mbwhere λ>0 is a given constant. Similar to the transformation of (3) into (4), we implement (8) by introducing a new variable w to bound the error and a variable xr to bound x and by relaxing the binary constraint on b, thus obtaining(9)minimizex,xr∈Rn,w∈Rm,b∈0,1n1mTw+λ1nTxrsubject  to1nTb≤K-w≤y-Dx≤wxr≤Mb-xr≤x≤xrRegularization also changes the form of the iteration step. The parameterα is now used as a weight value for the regularized error. The regularized version of (5) is(10)minimizex,xr∈Rn,w∈Rm,b∈0,1n1-α△b,b~+α1mTw+λ1nTxrsubject  to1nTb≤K-w≤y-Dx≤wxr≤Mb-xr≤x≤xrAdding the regularization term does not change the nature of the problem and the same algorithms as for SFP are used.We note that, unlike the Branch and Bound algorithms from [17], Algorithm 1 is not guaranteed to provide the exact solution of (1). However, it is significantly faster.


## 3. Results

For our tests, we use randomly generated dictionaries of size80×200 with condition numbers of 1000 and 100000; like in [17], we focus on not so well conditioned matrices because they are more difficult to handle by sparse representation algorithms. The test signals are obtained with y=Dxtrue+u, where the solutions xtrue have K∈5,7,9,11 nonzero coefficients generated randomly following a Gaussian distribution with zero mean and unit variance, in random positions; the noise u is Laplacian and its variance is chosen such that the signal to noise ratios (SNR) have values 10, 20, 30, and ∞. For each condition number, sparsity level K, and SNR, we generate 10 distinct dictionaries and solutions xtrue.For the computation of the representation error, the relative error(11)e=Dx-y1y1is used (where now x is the computed solution), in accordance with the formulation of the initial problem (1). The name of the algorithm may be added as subscript, like in eRLAD.We have implemented SFP and SFPreg as shown in Algorithm1. The initial weight α is set to 0.5 and is multiplied by an update factor γ=0.95 at each iteration; in [18, 21] the choices α=1 and γ=0.9 are used; by taking α=0.5 we aim to start with equal emphasis on the number of nonzero coefficients and on the representation error; this choice appears also to speed up the convergence. The number of iterations Iter is set to 1000. We run SFPreg with 50 equally spaced values of the regularization parameter λ from 0 to 1. The value for which the error is the smallest is considered the best.Our two algorithms are compared with the regularized least absolute deviation (RLAD) model [2](12)minimizex∈RnDx-y1+λ~∣x∥1solved with linprog. Note that this is a lasso-style relaxed problem, adapted to the l1 error norm. Problem (12) is solved for different values of λ~ until the solution x has a number of nonzero coefficients equal or close to K. After that, problem (12) is solved without the regularization term, keeping fixed the support already found. Hence, the solution is optimal for that support. This linear programming approach is generically named RLAD.The algorithms are implemented in MATLAB and tested on a computer with a 6-core 3.4 GHz processor and 32 GB of RAM.In Figure1 we show the effect of the regularization factor λ in SFPreg on error (11). The values for SFP (which corresponds to λ=0) are not shown, since they are much larger (with the exception of the noiseless case); in increasing SNR order, they are 0.51, 0.232, 0.108, and 2.7·10-12. It can be seen that as the SNR is smaller the error reduction is more significant. Although it is hard to give a general recipe on how λ should be chosen, it is visible that the error is small for a fairly large range of values. Hence, it is enough to try only few of them and take the one with best error.Figure 1
Mean representation error depending onλ for K=11, for matrix conditioning of 1000.Figures2 and 3 display the difference in relative errors between RLAD and SFPreg for different matrix condition numbers. It can be seen that the differences are significant. In order to obtain the same number of coefficients, RLAD is forced by the regularization to focus more on sparsity than on the error reduction. Both algorithms recover the solution exactly in the noiseless case.Figure 2
Difference in relative representation errors between RLAD and SFPreg for matrix conditioning of 1000. A positive value means that SFPreg is better.Figure 3
Difference in relative representation errors between RLAD and SFPreg for matrix conditioning of 100000.The relative errors (11) for the SFPreg algorithm are displayed in Figures 4 and 5. Obviously, the errors are dependent on the SNR. The error values are consistently at about the noise level, showing that SFPreg has a good behavior.Figure 4
Relative representation errors for SFPreg for matrix conditioning of 1000.Figure 5
Relative representation errors for SFPreg for matrix conditioning of 100000.The relative recovery errors computed with(13)erec=x-xtrue2xtrue2for the SFPreg algorithm are displayed in Figures 6 and 7. The recovery errors are also dependent on the SNR. For low perturbation the signal is recovered with very good accuracy.Figure 6
Relative recovery errors for SFPreg for matrix conditioning of 1000.Figure 7
Relative recovery errors for SFPreg for matrix conditioning of 100000.The means of the relative errors (11) obtained by running the tests are displayed in Tables 1 and 2. Table 1 shows the results for matrices with a condition number of 1000 and Table 2 for matrices with condition number 100000. The first value in each cell is the relative error for SFP, the second (with larger font) is for SFPreg, and the last is for RLAD. It can be seen that SFPreg has always smaller mean errors compared to RLAD. The SFPreg approach gives an important improvement, in some cases reducing the error by almost an order of magnitude compared to SFP. Also, SFP has many cases in which the representation error is much larger than for RLAD.Table 1
Mean relative errors (11) for SFP (top in each cell), SFPreg (middle), and RLAD (bottom) for matrix conditioning of 1000.
K 5 7 9 11 SNR = 10 0.262 0.379 0.427 0.471 0.255 0.255 0.248 0.241 0.447 0.450 0.500 0.471  SNR = 20 0.084 0.143 0.378 0.338 0.084 0.084 0.082 0.082 0.201 0.276 0.410 0.415  SNR = 30 0.027 0.027 0.095 0.157 0.027 0.026 0.025 0.025 0.068 0.231 0.270 0.464  SNR =∞ 1.814 ⋅ 1 0 - 11 4.556 ⋅ 1 0 - 12 1.037 ⋅ 1 0 - 12 1.161 ⋅ 1 0 - 12 5.801 ⋅ 1 0 - 14 2.646 ⋅ 1 0 - 14 1.985 ⋅ 1 0 - 14 1.132 ⋅ 1 0 - 14 5.603 ⋅ 1 0 - 10 6.772 ⋅ 1 0 - 10 9.165 ⋅ 1 0 - 10 8.621 ⋅ 1 0 - 10Table 2
Mean relative errors (11) for SFP (top in each cell), SFPreg (middle), and RLAD (bottom) for matrix conditioning of 100000.
K 5 7 9 11 SNR = 10 0.257 0.327 0.492 0.505 0.253 0.252 0.263 0.256 0.364 0.494 0.543 0.493  SNR = 20 0.086 0.086 0.185 0.445 0.085 0.082 0.081 0.079 0.143 0.224 0.438 0.406  SNR = 30 0.027 0.028 0.153 0.422 0.027 0.027 0.026 0.025 0.101 0.143 0.280 0.395  SNR =∞ 2.223 ⋅ 1 0 - 12 1.805 ⋅ 1 0 - 12 0.023 7.815 ⋅ 1 0 - 13 3.993 ⋅ 1 0 - 14 2.088 ⋅ 1 0 - 14 1.938 ⋅ 1 0 - 14 1.448 ⋅ 1 0 - 14 6.526 ⋅ 1 0 - 10 7.906 ⋅ 1 0 - 10 1.244 ⋅ 1 0 - 9 8.061 ⋅ 1 0 - 10In fact, the SFPreg algorithm has smaller representation errors in all cases with noise compared to the RLAD algorithm, for both conditioning numbers. It can be seen that the conditioning of the matrix does not affect the methods based on regularization.Turning now to the complexity of the algorithms, the mean number of iterations is 3.77 for SFP and 3.14 for SFPreg. So, regularization not only helps the algorithm to give better results but also leads to slightly faster convergence.The number of iterations depends on the selection ofα. A small α makes SFP and SFPreg converge faster but can cause a higher error. Also the update factor γ affects the speed and accuracy of the algorithms. The values reported above seem to provide a good compromise.The mean execution time of SFPreg is 0.416 seconds, while SFP requires 0.348 seconds. This is an important improvement compared to the formulation (4) solved with the Branch and Cut algorithm. As reported in [17], Branch and Cut require between 3 to 10000 seconds to find the solution for the problem. So, SFPreg is significantly faster and still obtains results that are nearly optimal.


## 4. Conclusions

In this paper an adaptation of the Feasibility Pump has been proposed to address the sparse representation problem withl1 norm of the error. As seen from the experimental results, the addition of the regularization term offers a considerable improvement to the algorithm. The Feasibility Pump gives a good solution in a very small number of iterations. Thus, the Feasibility Pump proves to be a good algorithm for finding sparse representations of signals. Due to regularization and the big-M trick, it limits the magnitude of the values of the coefficients and is thus appropriate also for ill-conditioned problems. Future research will focus on improving the randomization step, using different regularization terms, and adapting the algorithm for other sparse problems reformulations.



---

*Source: 5615243-2019-02-24.xml*