Sparse recovery: from vectors to tensors

Recent advances in various fields such as telecommunications, biomedicine and economics, among others, have created enormous amount of data that are often characterized by their huge size and high dimensionality. It has become evident, from research in the past couple of decades, that sparsity is a flexible and powerful notion when dealing with these data, both from empirical and theoretical viewpoints. In this survey, we review some of the most popular techniques to exploit sparsity, for analyzing high-dimensional vectors, matrices and higher-order tensors.


INTRODUCTION
The problem of sparse recovery is ubiquitous in modern science and engineering applications. In these applications, we are interested in inferring a high-dimensional object, namely a vector, a matrix or a higher-order tensor, from very few observations. Notable examples include identifying key genes driving a complex disease, and reconstructing high-quality images or videos from compressive measurements, among many others.
More specifically, consider linear measurements of an n-dimensional object x of the form: where ·, · stands for the usual inner product in R n , and a k are a set of prespecified n-dimensional vectors. The number of measurements m is typically much smaller than n so that the linear system (1) is underdetermined whenever m < n. Thus it is impossible to recover x from y k in the absence of any additional assumption. The idea behind sparse recovery is to assume that x actually resides in a subspace whose dimensionality is much smaller than the ambient dimension n.
A canonical example of sparse recovery is the socalled compressive sensing for vectors, where x is assumed to have only a small number of, albeit unknown, nonzero coordinates. More generally, we call a vector x ∈ R n k-sparse if it can be represented by up to k elements from a predetermined dictio-nary. Another common example is the recovery of low-rank matrices where x ∈ R n 1 ×n 2 is assumed to have a rank much smaller than min {n 1 , n 2 }. In many practical situations, we are also interested in signals of higher-order multilinear structure. For example, it is natural to represent multispectral images by a third-order multilinear array, or tensor, with the third index corresponding to different bandwidths. Clearly, vectors and matrices can be viewed as firstorder and second-order tensors as well. Despite the connection, moving from vectors and matrices to higher-order tensors could present significant new challenges. A common way to address these challenges is to unfold tensors to matrices; see e.g. [1][2][3][4]. However, as recently pointed out in [5], the multilinear structure is lost in such matricization and, as a result, methods based on these techniques often lead to suboptimal results.
A general approach to sparse recovery is through solving the following constrained optimization problem: min z S(z) subject to y k = a k , z , k = 1, · · · , m, (2) where S(·) is an objective function that encourages sparse solutions. The success of this approach hinges upon several crucial aspects including, among others, r how to choose the object function S(·); r how to solve the optimization problem (2); There is, by now, an impressive literature addressing these issues when x is a vector or matrix; see e.g. [6][7][8][9][10][11][12][13][14][15][16][17][18], among numerous others. In this review, we aim to survey some of the key developments, with a particular focus on applications to image and video analysis.
The rest of the paper is organized as follows. In the sections entitled 'Recovery of sparse vectors', 'Recovery of low-rank matrices' and 'Recovery of low-rank higher-order tensors', we discuss the recovery of sparse vector, low-rank matrix and low-rank tensor signals, respectively. Finally, a couple of illustrative examples in image and video processing are given in the section entitled 'Applications'.

RECOVERY OF SPARSE VECTORS
We first consider recovering a sparse vector signal, a problem more commonly known as compressive sensing [10][11][12][13].

Compressive sensing of sparse signals
With slight abuse of notation, let y be an mdimensional vector whose coordinates are the measurement y k , and A be an m × n matrix whose rows are given by a k . It is then not hard to see that (1) can be more compactly written as y = Ax, where x is an n-dimensional vector. Following the jargon in compressive sensing, hereafter we shall refer to A as the sensing matrix.
Obviously, when m < n, there may be infinitely many z that agree with the measurements in that y = Az. Since x is known a priori to be sparse, it is then naturally to seek among all these solutions the one that is sparsest. As mentioned before, an obvious way to measure sparsity of x is its 0 norm: where |·| stands for the cardinality of a set, leading to the following recovering x by a solution to min z∈R n z 0 subject to y = Az. ( Under mild regularity conditions on the sensing matrix A, it can be shown that the solution to (3) is indeed well defined and unique, and thus correctly recovers x [9]. However, it is also well known [19] that solving (3) is NP-hard in general and thus infeasible to compute for even moderate-size problems. The most popular way to overcome this challenge is the 1 relaxation, which minimizes the 1 norm instead, leading to min z∈R n z 1 subject to y = Az. (4) Assuming that x is k-sparse and the unique solution to (3), a key question pertaining to the 1 relaxation (4) is: under what condition is it also the unique solution to (4)? The answer can be characterized by various properties of the sensing matrix including the mutual incoherence property (MIP [20]), null space property (NSP [21]) and restricted isometry property (RIP [9]), among others. We shall focus primarily on RIP here. Interested readers are referred to [22,23] and references therein for further discussions on MIP and NSP. Definition 2.1 [9] A sensing matrix A ∈ R m×n is said to satisfy the RIP of order k if there exists a constant δ k ∈ [0, 1) such that, for every k-sparse vector z ∈ R n , Similarly, A is said to obey the restricted orthogonality property (ROP) of order (k, k ) if there exists a constant θ k,k such that, for every k-sparse vector z and k -sparse vector z with non-overlapping support sets, The constants δ k and θ k,k are called the k-restricted isometry constant (RIC) and (k, k )-restricted orthogonality constant (ROC), respectively. The concept was first introduced by Candès and Tao [9], who showed that, if δ k + θ k, k + θ k, 2k < 1, then x is the unique solution to (4). This condition has since been weakened. For example, Candès and Tao [13] showed that δ 2k + 3θ k, 2k < 2 suffices, and Cai et al. [24] required that δ 1.25k + θ k, 1.25k < 1. More recently, Cai and Zhang [25] further weakened the condition to δ k + θ k, k < 1 and showed that the upper bound 1 is sharp in the sense that, for any > 0, the condition δ k + θ k, k < 1 + is not sufficient to guarantee exact recovery. Sufficient conditions for exact recovery by (4) that only involve RIC have also been investigated in the literature. For example, [26] argued that δ 2k < √ 2 − 1 implies that x is the unique solution to (4). This was improved to δ 2k < 0.472 in [27]. More recently, [28] showed that, for any given constant t ≥ 4/3, the condition δ tk < (t − 1)/t guarantees the exact recovery of all k-sparse signals by 1 minimization. Moreover, for any > 0, δ tk < (t − 1)/t + is not sufficient to ensure exact recovery of all k-sparse signals for large k. Wang et al. 3 An immediate question following these results is how to design a sensing matrix that satisfies these conditions so that we can use (4) to recover x. It is now well understood that, for many random ensembles, that is, where each entry of A is independently sampled from a common distribution such as a Gaussian, Rademacher or other sub-Gaussian distribution, δ k < with overwhelming probability, provided that m ≥ C −2 k log (n/k), for some constant C > 0. There has also been some recent progress in constructing deterministic sensing matrices that satisfy these RIP conditions; see e.g. [29][30][31].

Compressive sensing of block-sparse signals
In many applications, the signal of interest may have more structured sparsity patterns. The most common example is so-called block-sparsity where sparsity occurs in a blockwise fashion rather than at the individual coordinate level. More specifically, let x ∈ R n be the concatenation of b signal 'blocks': x [1] x n 1 +1 · · · x n 1 +n 2 x [2] · · · x n−n b +1 · · · x n x [b] ] , (7) where each signal 'block' x[i] is of length n i . We assume that x is block k-sparse in that there are at most k nonzero blocks among x [i]. As before, we are interested in the most block-sparse signal that satisfies y = Az. To circumvent the potential computational challenge, the following relaxation is often employed: where the mixed 2 / 1 norm is defined as It is not hard to see that, when each block has size n i = 1, (8) reduces to the 1 minimization given by (4). More generally, the optimization problem in (8) is convex and can be recast as a second-order cone program, and thus can be solved efficiently. One can also extend the notion of RIP and ROP to the block-sparse setting; see e.g. [32,33]. For brevity, (9) Similarly, A is said to obey the block-ROP of order (k, k ) if there exists a constant θ k,k |I such that, for every block k-sparse vector z and block k -sparse vector z with disjoint supports, The constants δ k|I and θ k,k |I are referred to as block k-RIC and block (k, k )-ROC, respectively. Clearly, any sufficient RIP conditions of standard 1 minimization can be naturally extended to the setting of block-sparse recovery via mixed 2 / 1 minimization so that, for example, θ k,k|I + δ k|I < 1 is also a sufficient condition for x to be the unique solution to (8).

Nonconvex methods
In addition to the 1 -minimization-based approach, there is also an extensive literature on nonconvex methods for sparse recovery where instead of the 1 norm of z one minimizes a nonconvex objective function in z. The most notable example is the q (0 < q < 1) (quasi-)norm, leading to min z∈R n z q subject to y = Az.
Other notable examples of nonconvex objective functions include smoothly clipped absolute deviation (SCAD) [40] and the minimax concave plus function [41], among others.

Compressive sensing with general dictionaries
Thus far, we have focused on sparsity with respect to the canonical basis of R n . In many applications, it might be more appropriate to have sparsity with respect to more general dictionaries; see e.g. [42][43][44][45][46]. More specifically, a signal x is represented by x = Dα with respect to a dictionary D ∈ R n×n , where α ∈ R n is the coordinate in the dictionary and is known a priori to be sparse or nearly sparse. One can obviously treat A = AD as the new sensing matrix and apply any of the aforementioned methods to exploit the sparsity of α. One of the drawbacks, however, is that nice properties of the original sensing see e.g. [43,44,47,48].

Compressive phase retrieval
In the previous subsections, we have mainly discussed the problem of recovering a sparse signal from a small number of linear measurements. However, in some practical scenarios one can only observe some nonlinear measurements of the original signal. The typical example is the so-called compressive phase retrieval problem. In such a scenario, we observe the magnitude of the Fourier coefficients instead of their phase, which can be modeled as the form For recovering x, several studies [49][50][51] have considered the following 1 minimization: By introducing a strong notion of RIP, Voroninski and Xu [51] built a parallel result for compressive phase retrieval with the classical compressive sensing. Specifically, they proved that a k-sparse signal x can be recovered from m = O(k log (n/k)) random Gaussian phaseless measurements by solving (14). Unlike the standard convex 1 minimization (4), the problem (14) is a nonconvex problem, but some efficient algorithms have been developed to compute it; see e.g. [50,52].

RECOVERY OF LOW-RANK MATRICES
We now consider recovering a low-rank matrix, a problem often referred to as matrix completion.

Matrix completion via nuclear norm minimization
In many practical situations such as collaborative filtering, system identification and remote sensing, to name a few, the signal that we aim to recover oftentimes is a matrix rather than a vector. To signify this fact and distinguish from the vector case treated in the last section, we shall write the underlying signal as a capital letter, X ∈ R n 1 ×n 2 , throughout this section. In these applications, we often observe a small fraction of the entries of X. The task of matrix completion is then to 'complete' the remaining entries. Formally, let be a subset of [ The goal of matrix completion is to recover X based on {X ij : (i, j) ∈ }, particularly with the sample size | | much smaller than the total number n 1 n 2 of entries. To fix ideas, we shall assume that is a uniformly sampled subset of [n 1 ] × [n 2 ], although other sampling schemes have also been investigated in the literature, e.g. [53][54][55]. Obviously, we cannot complete an arbitrary matrix from a subset of its entries. But it is possible for low-rank matrices as their degrees of freedom are much smaller than n 1 n 2 . But low rankness alone is not sufficient. Consider a matrix with a single nonzero entry; it is of rank one but it is impossible to complete it unless the nonzero entry is observed. A formal way to characterize low-rank matrices that can be completed from {X ij : (i, j) ∈ } was first introduced in [14].
Definition 3.1. Let U be a subspace of R n of dimension r and P U be the orthogonal projection onto U. Then the coherence of U (with respect to the standard basis (e i )) is defined to be It is clear that the smallest possible value for μ(U) is 1 and the largest possible value for μ(U) is n/r. Let M be an n 1 × n 2 matrix of rank r and with column and row spaces denoted by U and V, respectively. We shall say that M satisfies the incoherence condition with parameter μ 0 if max (μ(U), μ(V)) ≤ μ 0 . Now let X be an incoherent matrix of rank r. In a similar spirit as the vector case, a natural way to reconstruct it from {X ij : (i, j) ∈ } is to seek, among all matrices whose entries indexed by agree with our observations, the one with the smallest rank: Again, to overcome the computational challenges in directly minimizing matrix ranks, the following convex program is commonly suggested: where the nuclear norm · * is the sum of all singular values. As before, we are interested in when the solution to (16) is unique and correctly REVIEW Wang et al. 5 recovers x. Candès and Recht [14] were the first to show that this is indeed the case, for almost all that are large enough. These results are probabilistic in nature due to the randomness of , that is, one can correctly recover x using (16) with high probability regarding min {n 1 , n 2 }. The sample size requirement to ensure exact recovery of X by the solution of (16) was later improved in [16] to | | ≥ Cr(n 1 + n 2 ) · polylog(n 1 + n 2 ), where C is a constant that depends on the coherence coefficients only. It is clear to see that this requirement is (nearly) optimal in that there are O(r(n 1 + n 2 )) free parameters in specifying a rank-r matrix.

Matrix completion from affine measurements
More generally, one may consider recovering a lowrank matrix based on affine measurements. More specifically, let A : R n 1 ×n 2 → R m be a linear map such that A(X ) = y . We aim to recover X based on the information that A(X ) = y . It is clear that the canonical matrix completion problem discussed in the previous subsection corresponds to the case when A(X ) = {X i j : (i, j ) ∈ }. Similarly, we can proceed to reconstruct X by the solution to min Z∈R n 1 ×n 2 Z * subject to y = A(Z).
(17) It is of interest to know under which kind of sensing operator A can X be recovered in exactly this way. An answer is given in [56], which extends the concept of RIP to general linear operator for matrices. Definition 3.2. A linear operator A : R n 1 ×n 2 → R m is said to satisfy the matrix RIP of order r if there exists a constant δ Z r such that holds for all matrices Z ∈ R n 1 ×n 2 of rank at most r.
Besides the aforementioned low-rank matrix completion problems, some recent studies, e.g. [58][59][60], have drawn attention to the so-called high-rank matrix completion, in which the columns of the matrix belong to a union of subspaces and, as a result, the rank can be high or even full.

Mixed sparsity
In some applications, the matrix that we want to recover is not necessarily of low rank but differs from a low-rank matrix only by a small number of entries; see e.g. [17,18]. In other words, we may write X = S + L, where S is a sparse matrix with only a small number of nonzero entries while L is a matrix of low rank. It is clear that, even if we observe X entirely, the decomposition of X into a sparse component S and a low-rank component L may not be uniquely defined. General conditions under which such a decomposition is indeed unique are provided in, for example, [18]. In the light of the previous discussions, it is natural to consider reconstructing X from observations {X ij : (i, j) ∈ } by the solution to min Z 1 ,Z 2 ∈R n 1 ×n 2 This strategy has been investigated extensively in the literature; see e.g. [17]. Further developments in this direction can also be found in [61][62][63][64], among others.

Nonconvex methods
Just as 1 norm is a convex relaxation of the 0 norm, the nuclear norm is also a convex relaxation of the rank for a matrix. In addition to these convex approaches, nonconvex methods have also been proposed by numerous authors. The most common example is the Schatten-q (0 < q < 1) (quasi-)norm defined by where σ 1 , σ 2 , ··· are the singular values of X. It is clear that, when q → 0, X sq → rank(X) while the nuclear norm corresponds to the case when q = 1.

RECOVERY OF LOW-RANK HIGHER-ORDER TENSORS
In an increasing number of modern applications, the object to be estimated has a higher-order tensor structure. Typical examples include video inpainting [69], scan completion [70], multichannel EEG (electroencephalogram) compression [71], traffic data analysis [72] and hyperspectral image restoration [73], among many others. Similar to matrices, in many of these applications, we are interested in recovering a low-rank tensor either from observing a subset of its entries or a collection of affine measurements.
Despite the apparent similarities between matrices and higher-order tensors, it is delicate to extend the idea behind nuclear norm minimization to the latter because matrix-style singular value decomposition does not exist for higher-order tensors. A common approach is to first unfold a higher-order tensor to a matrix and then apply a matrix-based approach to recover a tensor. Consider, for example, a thirdorder tensor X ∈ R d 1 ×d 2 ×d 3 . We can collapse its second and third indices, leading to a d 1 × (d 2 d 3 ) matrix X (1) . If X is of low rank, then so is X (1) . We can exploit the low rankness of X (1) by minimizing its nuclear norm. Clearly we can also collapse the first and third indices of X , leading to a matrix X (2) , and the first and second, leading to X (3) . If we want to complete a low-rank tensor X based on its entries X (ω) 3 ], we can then consider recovering X by the solution to the following convex program: Efficient algorithms for solving matrix nuclear norm minimization (16) such as the alternating direction method of multipliers (ADMM) and Douglas-Rachford operator splitting methods can then be readily adapted to solve (22); see e.g. [3,4].
As pointed out in [5], however, such an approach fails to exploit fully the multilinear structure of a tensor and is thus suboptimal. Instead, directly minimizing the tensor nuclear norm was suggested: (23) where the tensor nuclear norm · * is defined as the dual norm of the tensor spectral norm · . For more discussions on the tensor nuclear and spectral norms, please see [74]. Unfortunately, unlike the matrix nuclear norm, computing the tensor nuclear norm, and thereby the problem (23), is NP-hard. Hence, various relaxations and approximate algorithms have been introduced in the literature; see e.g. [75][76][77][78][79], and references therein. This is a research area in its infancy and many interesting issues need to be addressed.

APPLICATIONS
In the previous sections, we have given an overview of some basic ideas and techniques in dealing with sparsity in vectors and low rankness in matrices and tensors. We now give a couple of examples to illustrate how they can be used in action.

Background subtraction with compressive imaging
Background subtraction in image and video has attracted a lot of attention in the past couple of decades. It aims at simultaneously separating video background and extracting the moving objects from a video stream, and can provide important clues for various applications such as moving object detection [80] and object tracking in surveillance [81], among numerous others.
Conventional background subtraction techniques usually consist of four steps: video acquisition, encoding, decoding and separating the moving objects from the background. This scheme needs to fully sample the video frames with large computational and storage requirements, followed by well-designed video coding and background subtraction algorithms. To alleviate the burden of computation and storage, a newly-developed compressive imaging scheme [82][83][84] has been used for background subtraction by combining the video acquisition, coding and background subtraction into a single framework, as illustrated in Fig. 1. It simultaneously achieves background subtraction and video reconstruction. In this setting, the main objective is then to maximize the reconstruction and separation accuracies using as few compressive measurements as possible.
Several studies have been carried out for background subtraction from the perspective of compressive imaging. In the seminal work of Cevher et al. [85], the background subtraction problem is formulated as a sparse recovery problem. They showed that the moving objects can be recovered by learning a low-dimensional compressed representation of the background image. More recently, the robust principle component analysis (RPCA) approach has also been used to deal with the problem of background subtraction with compressive imaging, in which they commonly model the video as a matrix with columns of vectorized video frames and then decompose the matrix into a low-rank matrix L and a sparse matrix S; see e.g. [86][87][88][89].
Although methods based on RPCA have achieved satisfactory performance, they fail to exploit the finer structures of background and foreground after vectorizing the video frames. It appears more advantageous to model the spatio-temporal information of background and foreground using direct tensor representation of a video. To this end, a novel tensor RPCA approach has been proposed in [90] for background subtraction from compressive measurements by decomposing the video into a static background with spatiotemporal correlation and a moving foreground with spatio-temporal continuity within a tensor representation framework. More specifically, one can use 3D total variation (3D-TV) to characterize the spatio-temporal continuity underlying the video foreground, and low-rank Tucker decomposition to model the spatio-temporal correlation of the video background, which leads to the following tensor RPCA model: where the factor matrices U 1 and U 2 are orthogonal in columns for two spatial modes, the factor matrix U 3 is orthogonal in columns for the temporal mode, the core tensor G interacts with these factors and the 3D-TV term · 3D-TV is defined as Because a 3D patch in a video background is similar to many other 3D patches over the video frames, one can model the video background using several groups of similar video 3D patches, where each patch group corresponds to a fourth-order tensor. Integrating the patch-based modeling idea into (24), one can easily get a patch-group-based tensor RPCA model. It should be noted that solving the nonconvex tensor RPCA model (24) as well as its patchgroup-based form are computationally difficult. In practice, we can find a local solution using a multiblock version of ADMM. For more details, please refer to Section V of [90]. Fig. 2 gives an example based on three real videos. It is evident that the proposed tensor models enjoy a superior performance over the other popular matrix models both in terms of the quality of the reconstructed videos and in terms of the separation of the moving objects. This suggests that using a direct tensor modeling technique to deal with practical higherorder tensor data can utilize more useful structures.

Hyperspectral compressive sensing
Hyperspectral imaging employs an imaging spectrometer to collect hundreds of spectral bands  [90] and three popular matrix models (i.e. SparCS [86], ReProcs [87] and SpLR [88]) under the sampling ratio 1/30. The first column shows the original video frames from different video volumes (a)-(c); the second to sixth columns correspond to the results produced by all the compared methods, respectively. Here, for each method, the reconstruction result of the original video frame (upper panels) and the detection result of moving objects in the foreground (lower panels) is shown.
ranging from ultraviolet to infrared wavelengths for the same area on the surface of the Earth. It has a wide range of applications including environmental monitoring, military surveillance and mineral exploration, among numerous others [91,92]. Figuratively speaking, a hyperspectral image can be treated as a 3D (x, y, λ) data cube, where x and y represent two spatial dimensions of the scene, and λ represents the spectral dimension comprising a range of wavelengths. Typically, such hyperspectral cubes are collected by an airborne sensor or a satellite and sent to a ground station on Earth for subsequent processing. Noting that the dimension of λ is usually in the hundreds, hyperspectral cubes are of a fairly large size even for moderate dimensions of x and y. This makes it necessary to devise effective techniques for hyperspectral data compression, due to the limited bandwidth of the link connection between the satellite/aerospace and the ground station.
In the last few years, a popular hyperspectral compressive sensing (HCS) scheme based on the principle of compressive sensing for hyperspectral data compression, as illustrated in Fig. 3, has been ex-tensively investigated. Similar to other compressive sensing problems, the main objectives of HCS are to design easy hardware encoding implementation and develop an efficient sparse reconstruction procedure. In what follows we shall focus primarily on the latter. For hardware implementation, interested readers are referred to [93,94] and references therein.
Like other natural images, hyperspectral images can be sparsified by using certain transformations. Traditional sparse recovery methods such as 1 minimization and TV minimization are often used for such purposes; see e.g. [95,96]. To further exploit the inherent spectral correlation of hyperspectral images, a series of work based on low-rank modeling has been carried out in recent years. For example, Golbabaee and Vandergheynst proposed in [97] a joint nuclear and 2 / 1 minimization method to describe the spectral correlation and the joint-sparse spatial wavelet representations of hyperspectral images. Then, they modeled the spectral correlation together with the spatial piecewise smoothness of hyperspectral images by using a joint nuclear and   Here the sampling ratio is 1%. The last column shows the original image bands. The columns from the first to the fourth correspond to the produced results from the KCS method [96], the JNTV method [98], the ST-NCS method [100] and the tensor method (25), respectively.
TV norm minimization method [98]. Similar to the surveillance videos, however, modeling a hyperspectral cube as a matrix cannot utilize the finer spatialand-spectral information, leading to suboptimal reconstruction results under relatively low sampling ratios (e.g. 1% of the whole image size).
To further exploit the compressibility underlying a hyperspectral cube, one may consider such a single cube as a tensor with three modes (width, height and band) and then identify the hidden spatialand-spectral structures using direct tensor modeling techniques. Precisely, all the bands of a hyperspectral image have very strong correlation in the spectral domain, and each band, if considered as a matrix, has relatively strong correlation in the spatial domain; such spatial-and-spectral correlation can be modeled through low-rank Tucker decomposition. In addition, the intensity at each voxel is likely to be similar to its neighbors, which can be characterized by smoothness using the so-called 3D-TV penalty. In summary, [99] considered the following joint tensor Tucker decomposition and 3D-TV minimization model: It is clear that the above minimization problem is highly nonconvex. One often looks for good local solutions using a multi-block ADMM algorithm. Fig. 4 gives an example of the first band of four hyperspectral datasets, respectively reconstructed by different methods, with the sampling ratio at 1%. It is evident that the proposed tensor method could provide nearly perfect reconstruction. In addition, REVIEW sparse tensor and nonlinear compressive sensing (ST-NCS) performs slightly better than Kronecker compressive sensing (KCS) and joint nuclear/TV norm minimization (JNTV) in terms of reconstruction accuracy, because of the use of a direct tensor sparse representation of a hyperspectral cube. Both these findings demonstrate the power of tensor modeling techniques. It is also worth noting that, compared with ST-NCS, the images reconstructed with the method (25) are clearer and sharper.