Piecewise Sparse Approximation via Threshold P_OMP with Application to Surface Reconstruction

Sparse signal representations have gained extensive study in recent years. In applications, there are large amounts of signals that are structured. Motivated by signal decomposition and scattered data reconstruction applications, we consider a particular type of structured signals which can be represented by a union of several sparse vectors. We deﬁne this type of signal as piecewise sparse signals. To ﬁnd a piecewise sparse representation of a signal, we propose a thresholding version of piecewise orthogonal matching pursuit(TP OMP), which aims to overcome the disadvantages of P OMP. We also establish the connection of piecewise sparsity and sampling over a union of subspaces. We evaluate the performance of the proposed greedy programs through simulations on surface reconstruction.


Introduction
Sparse recovery or sparse representation from a given signal(or image) b with linear noisy measurement is of interest in many different applications. The matrix A ∈ R m×n with m ≪ n and e is an additive noise. Large amounts of algorithms have been proposed to determine x in a stable and efficient manner, for detail refer to [1,2,3,4,5,6,7]. In standard sparse recovery problem, the undermined vector x is assumed to be sparse, i.e, contains at most s nonzero entries(s ≪ m) without any further information on the structure of x. However, in real applications the nonzero elements of x are always correlative, thus x possess certain type of structure. Research in recent years study the sparsity properties of x together with its structure, which result in structured sparsity [8,9,10]. The structured sparsity has the merit of reducing the degree of freedom of the solution thereby better recovery capacity. Furthermore, one can use the structured sparsity of the signal to reduce the sampling rate, thus attain the goal of reduction of storage in applications. Particularly, we use "piecewise sparse" in [11,12] with the following definition . , x T N ) T ∈ R n is partitioned into N components and it is assumed that every x T i ∈ R n i containing nonzero entries is sparse, s i = x i 0 , i = 1, . . . , N. We call the vector x is (s 1 , . . . , s N )-piecewise sparse.
Piecewise sparsity describes a type of structured sparsity which emerges in the application of reconstructing a surface from scattered data in PSI space [13,14,15], data separation/signal decomposition [16,17,18,19], etc. In these applications, data(including image, signal) are represented in different bases and frames, and each frame exhibits its particular geometric tendency or attribute, i.e. signal b = P + C = Φx + Ψy, where Φ and Ψ are different frames and data separation aims at obtaining sparse coefficients x and y.
Instead of solving data separation problem, we discuss the approximation problem based on the special sparsity-piecewise sparse summarized from data separation: the dictionary is a union of basis A = [A 1 , . . . , A N ] and we describe the piecewise sparse approximation as the following problem For each frame A i , x i is a sparse vector, thus x = (x T 1 , . . . , x T N ) T is a piecewise sparse vector. It is shown in [12,11] that better recovery and estimate is possible when A is in an union manner, which provide milder conditions for estimating the sparse recovery results. Sparse representations in unions of orthogonal bases are studied in [20,21,22]. We generalize the results in orthogonal cases to the cases of general bases in [12], improved bounds of x * 0 recovered by both l 0 and l 1 optimizations. Furthermore, piecewise inverse scale space method and piecewise sparse OMP(P OMP) are provided in [11,23] for piecewise sparse recovery based on l 1 and l 0 minimization,respectively, both aim at preserving the smallest magnitudes. P OMP(Piecewise Orthogonal Matching Pursuit) is a type of greedy algorithm proposed by the motivation of a combination of the Compressive Sampling Matching Pursuit(CoSaMP) [5] and Generalized OMP(GOMP) [24]. However, lots of prior information such as piecewise sparsity should be provided before the algorithm started. The focus of the present paper is on developing a generalized greedy algorithm for piecewise sparse recovery which can be both applied to large amounts of signal decomposition applications and do not require prior information on the signal. Our algorithm is based on the P OMP algorithm and the main idea of threshold greedy algorithms.
Organization The remaining of the paper is organized as follows. The P OMP algorithm is described in Section 2. A thresholding version of P OMP which do not require piecewise sparsity in prior named TP OMP algorithm is also proposed in Section 2. We then show the connection of piecewise sparsity and sampling over a union of subspaces. In Section 3 we show a practical application of union sampling-surface reconstruction with scattered data and apply the results to surface reconstruction.

P OMP Algorithm
The main iteration of P OMP algorithm in [23] are as follows: It is obvious the major difference between P OMP and CoSaMP in [5] lie in the Identification and Pruning steps. It is noticed that the P OMP locate s entries of a vector in step Identification similar to the first step of GOMP. However, the P OMP locate the largest s i entries of x i (i = 1, . . . , N) and form all the components together as the s = (s 1 + . . . + s N ) entries, which is different from both CoSaMP and GOMP. Instead of selecting s-largest components of least-square solution, the P OMP sort the components of the piecewise vectorx = (x 1 , . . . ,x N ) T solved by least-square by magnitudes of each piece. The piecewise pruning process make progress in preserving small-scaled entries in each piece from false chosen cause of noise perturbation.
Though the P OMP is designed especially for piecewise sparse recovery, it requires prior informationpiecewise sparsity level s i for i = 1, . . . , N which is often unknown in applications. Actually, the majority of the greedy algorithms for sparse recovery require the global sparsity s before algorithms started. However, the purpose of signal decomposition applications is recover sparse representation for each components(subvector) without known sparsity level. Thus it is necessary to develop a generalized piecewise greedy algorithm for piecewise sparse recovery.

Threshold P OMP Algorithm
In order to avoid the disadvantage P OMP, we develop a threshold P OMP algorithm which use threshold to adaptively select atoms in iterations as follows: Notice that the length of x i for each piece are known in prior in signal decomposition problem, we denote the length of x i as d i .
Thresholding Strategy In practice, hard thresholding leads to better results. Furthermore, Bobin etc. showed that the use of hard thresholding is likely to provide the l 0 -sparsest solution in [25].
In algorithm 2 the T H γ (u) denotes component-wise thresholding with threshold γ: hard thresholding T H γ (u) = u if |u| > γ and zero elsewhere.

Algorithm 2 Threshold P OMP(TP OMP) Recovery Algorithm
Thresholding Parameters The thresholding parameters in Algorithm 2 are set as follows: By using the thresholding strategy, the TP OMP algorithm do not require the prior information on piecewise sparsity or global sparsity, thus it is more practical.

Sampling Over a Union of Subspaces
In this section we describe the connection of piecewise sparsity and sampling over a union of subspaces. Therefore we show the application to sparse approximation of fitting surface to scattered points, which is an application of sampling over a union of subspaces.
Large amounts of reference such as [26,27] have illustrate in detail on the traditional sampling theory which recover a unknown signal x ∈ H from m samples. It is common assumed that x lie in a particular subspace A of H . The general linear sampling model is: b l = f l , x , 1 ≤ l ≤ m with f l is a series of function in H . When H = R n the sampling system can be rewritten as: where the matrix F is formed by columns of f l . It is generally assumed that in [26,27,28] x lie in a given subspace A of H . However, in applications x often do not lie in particular subspace but in a union of subspaces, such as piecewise polynomials, sparse representation, etc. [29] with the following description: where V l are subspaces.
Connection of Block sparsity and Sampling Over a Union of Subspaces Eldar and Mishali considered a particular case of (4) in [30], i.e., V i = l=s A l where {A l , 1 ≤ l ≤ N} is a series of disjoint subspaces, |l| = s represents the sum of s indices which corresponds to the block sparsity in [30].
is block sparse represented since c is s-block sparse.
Connection with Piecewise Sparse Recovery We consider another type of sampling: with a piecewise sparse vector c. It is obvious the system is equivalent to piecewise sparse recovery problem in Definition 1 by setting A = F T D.

Experiments
In science and engineering fields, scattered data reconstruction refers to the problem of fitting a smooth surface through a scattered, or nonuniform, distribution of data samples. In a typical scattered data reconstruction problem, there is a set of scattered data sites X = {x 1 , x 2 , . . . , x n } ⊆ Ω ⊆ R d and corresponding function values f | X = { f 1 , f 2 , . . . , f n }, and we seek a function g belonging to space H which fits the given data {(x k , f k )} n k=1 . A classical approach to scattered data approximation is the smoothing spline for multidimensional data (i.e. d ≥ 2), which solves the minimization: where g belongs to the Beppo−Levi space H m [15]. However, this method can be computationally expensive as the size of data grows large. Johnson et al. [15] proposed a regularized least square framework in an approximation space spanned by the shifts and dilates of a single compactly supported function (PSI space) using the uniform B-spline tool. By using the space H = N i=1 H i as the approximation space, where H i (⊆ H i+1 ) is a principle shift invariant (PSI) space generated by a single B-spline function φ i . The fitting surface g is represented as: where c i j are coefficients components, and H i = span{φ i j , j = 1, ..., m}. In accordance to the piecewise sparse model, the matrix where h > 0 is a scaling parameter that control the refinement of the PSI space, and K is the index set It is obvious this problem is a typical scheme of sampling in a union of spaces. In another word, reconstruction of a fitting surface g(x, y) with a sparse representation for a given set of scattered points {(x k , y k )} and the corresponding data { f (x k , y k )} with noise. The approximate function is represented by the multi-level B-spline bases φ i j (x, y), i = 1, . . . , N; j = 1, · · · , n i , where N means the number of levels of B-splines, and n = n 1 + · · · + n N is the total number of the bases. That is where the coefficients c = (c 1 1 , . . . , c 1 n 1 , c 2 1 , . . . , c 2 n 2 , . . . , c N 1 , . . . , c N n N ) T is the N piecewise sparse vector to be determined. Thus the problem is also a piecewise sparse recovery probelm. We use the 2D tensor product quadratic B-spline with uniform knots as the function φ in this paper. For a given set of scattered points (x i , y i ) n i=1 ⊆ Ω ⊆ R 2 and the corresponding noisy data { f i } n i=1 . Define the 2D tensor product quadratic B-spline function with uniform knots is: where P i j are the control points and N i,2 (x), N j,2 (y) are the B-spline basis function corresponding to the uniform grid on the x and y axis. Given a grid U : . ., the B-spline basis function is defined as: and

Numerical Examples
In this part, we apply TP OMP to reconstruct a fitting surface g(x, y) with a sparse representation for a given set of 30 × 30 scattered points {(x k , y k )} 900 k=1 ⊆ Ω = [−1, 1] × [−1, 1] and the corresponding data { f (x k , y k )} 900 k=1 with noises (ε ∈ N(0, σ 2 I n )), compared with the approximation surface reconstructed by Multilevel Least Absolute Shrinkage and Selection Operator(MLASSO) in [31] and CoSaMP in [5]. The approximate function is represented by the multi-level B-spline bases φ i j (x, y), i = 1, . . . , N; j = 1, · · · , n i , where N means the number of levels of B-splines, and n = n 1 + · · · + n N is the total number of the bases. That is where the coefficients c = (c 1 1 , . . . , c 1 n 1 , c 2 1 , . . . , c 2 n 2 , . . . , c N 1 , . . . , c N n N ) T is the N block vector to be determined. Since all bases φ i j (x, y) are linear dependent, so we want to find a sparse solution among all solutions which minimize the fitting error ∑ 900 k=1 | f (x k , y k ) − g(x k , y k )| 2 . Furthermore, due to the different scaled supports of the multi-level B-spline bases, the different scaled features of g(x, y) can be expressed by the corresponding multiple coefficients c.
We test the following functions, with N = 4 levels of B-spline bases.   The difference between f (x, y) and g(x, y) is measured by the normalized root mean square (RMS) error [13]: Parameters Setting In this part, the parameters in TP OMP are set as C 1 = 50 and C 2 = 1.2. We use the regularization parameter selection rule in [31] to set λ 1 = 0.1, λ 2 = 0.02, λ 3 = 0.01, and λ 4 = 0.1 for the MLASSO model min Since the global sparsity s is required before the CoSaMP started(which is the disadvantage of CoSaMP), we set s = 50 in general for CoSaMP, which is equivalent to the maximum iteration number of TP OMP.

Results and Discussion
Remark 1 We claim that we do not compare the performance of P OMP with TP OMP since the latter is an thresholding version of P OMP with the merit of less prior information, the numerical performance of TP OMP is similar to P OMP.    Figure 1-3 shows the approximation performance with respect to function f 1 , f 2 and f 3 . It is observed that the approximation surface recovered by MLASSO misses some details in Figure 1, and the approximation surface recovered by CoSaMP is more sensitive to noise than that of MLASSO and TP OMP. Both the approximation surfaces recovered by MLASSO and CoSaMP fluctuate due to the noise in Figure 2. The peak of MLASSO in Figure 3 can not reach the original peak. For three examples, the approximation surface recovered by TP OMP approximate the scattered data well.
It is observed from Tab.1 that TP OMP's superiority in terms of running time. We use piecewise sparsity in order to show the piecewise sparse structure recovery results by different algorithms. Smaller but not all zeros piecewise sparsity implies better recovery quality. Due to the highly dependency on λ i , MLASSO always fail to result good piecewise sparsity. Experimental results demonstrate that the TP OMP algorithm exhibits competitive performance in the surface reconstruction, while with fast running time and better piecewise sparsity.

Conclusion
In this paper, we exploit the piecewise structure of signal in sparse recovery problem, based on which we propose a threshold greedy algorithm motivated by the P OMP algorithm. The proposed thresholding algorithm enjoys the merits of less required prior information and less computational cost. We then present numerical experiments in application to reconstructing a surface with scattered data in order to demonstrate the efficiency of our algorithm.