Modified OMP algorithm for exponentially decaying signals.

A group of signal reconstruction methods, referred to as compressed sensing (CS), has recently found a variety of applications in numerous branches of science and technology. However, the condition of the applicability of standard CS algorithms (e.g., orthogonal matching pursuit, OMP), i.e., the existence of the strictly sparse representation of a signal, is rarely met. Thus, dedicated algorithms for solving particular problems have to be developed. In this paper, we introduce a modification of OMP motivated by nuclear magnetic resonance (NMR) application of CS. The algorithm is based on the fact that the NMR spectrum consists of Lorentzian peaks and matches a single Lorentzian peak in each of its iterations. Thus, we propose the name Lorentzian peak matching pursuit (LPMP). We also consider certain modification of the algorithm by introducing the allowed positions of the Lorentzian peaks' centers. Our results show that the LPMP algorithm outperforms other CS algorithms when applied to exponentially decaying signals.

Introduction

Compressed sensing (CS) theory provides the methods to solve an underdetermined linear equation = within a set of sparse vectors x ∈ ℂ. The unconstrained linear equation for an m × n matrix of a maximum rank may have infinitely many solutions for m < n. However, the number of sparse solutions happens to be finite. Vectors in ℂ having at most l ∈ ℕ non-zero coordinates are referred to as the l-sparse vectors.

All l-sparse vectors fulfilling the equation = may be found by considering linear sub-problems of solving = with the constraint ∈ ℂ, where ℂ ⊂ ℂn denotes the subspace of vectors supported by I ⊂ [1,n]. This provides all l-sparse solutions, but the overall problem was proven to be NP-hard. Remarkably, the sparsest vector solving = may be found as the solution of the l1-minimization problem: where ‖ ‖ 1 denotes the l-norm of . The l-norm convexity and the linear constraint = allow the application of linear programming algorithms, as long as has a restricted isometry property; see [1].

A number of CS algorithms has been proposed, including: iterative soft thresholding (IST) [2], orthogonal matching pursuit (OMP) [3], compressive sampling matching (COSAMP) [4], stage-wise orthogonal matching pursuit (StOMP) [5] and iterative re-weighted least squares (IRLS) [6]. Inspired by the greedy CS algorithm, OMP, and motivated by a case of nuclear magnetic resonance (NMR) spectroscopy, we propose an approach that is designed specifically to recover a signal being a combination of exponentially decaying components. The spectrum of such a signal consists of Lorentzian peaks. The algorithm matches a single peak in each of its iterations, and thus, the name Lorentzian peak matching pursuit (LPMP) is proposed. The peaks' centers and heights are found as in the OMP algorithm, whereas the widths of the peaks are determined in a novel way.

NMR Spectroscopy and CS

The sources of the free induction decay (FID) signal measured in NMR spectroscopy are magnetic moments of atomic nuclei with non-zero spin, e.g., 1H,13C,15N,19F, that undergo coherent precession when excited by radio-frequency irradiation in the presence of a high static magnetic field [7]. The time dependence of an ideal FID signal f(t) is given by the following linear combination: of the exponentially decaying oscillations: where N is a number of nuclei groups differing in resonance frequency; a are amplitudes (the (x0221D) number of nuclei in the k-th group); γ is a relaxation rate for the k-th group; and ω is an oscillation frequency of the k-th group.

The precession frequencies (and thus, the frequencies of the emitted FID signal) are in the order of MHz, but their band is very narrow, typically a few kHz. In such a narrow range, tens or even hundreds of peaks are present. Often, the spectral resolution (peak width compared to the difference in resonance frequencies) is too low to resolve peaks and, thus, makes an analysis difficult. Fortunately, the peak dispersion can be improved by the acquisition of the multidimensional NMR signal [8] where additional (“indirect”) time dimensions are introduced. Peaks in the multidimensional NMR spectrum correspond to groups of nuclei coupled by some physical interaction. Using the above notation, the formula of the two-dimensional FID is: where t1 is indirectly sampled time and t2 directly sampled time.

The reason for the application of CS in multidimensional NMR spectroscopy is a time-consuming sampling of the indirect dimensions of the FID signal (t1 in Equation (2)). There are three reasons for the lengthy acquisition. Firstly, each measurement point in the indirect dimension takes a few seconds of “real” time. Secondly, the sampling density has to fulfil the Nyquist theorem [9]. Finally, the sampling has to reach far into the indirect time domain to provide sufficient spectral resolution [10]. The latter, together with the Nyquist condition for the sampling density, means that, often, many thousands of sampling points have to be acquired. This makes an experiment unacceptably long (days in the case of three or four dimensions) or completely impossible in the case of higher dimensionality (5D and more). Thus, there is a need for alternative sampling approaches that allow one to save expensive spectrometer time preserving the spectral information. Most of them use the non-uniform sampling, where the major part of the points is not measured, but reconstructed after the experiment based on certain assumptions [11,12]. The CS, where the assumption is spectral sparsity, has been successfully applied in NMR exploiting the algorithms based on convex l1-norm or non-convex l-norm minimization [13,14]. The history of the application of the “greedy” algorithm, OMP, in NMR is long; in fact, the algorithm was transferred from radio astronomy, where it existed under the name, CLEAN [15]. Direct application of CLEAN in NMR is possible [16], but far from being perfect, because of the aforementioned approximate sparsity of the signal. Thus, various modifications of CLEAN have been proposed [17-19]. Importantly, none of these works refer to OMP formalism and do not explain the modifications by the sparsity criteria. We propose an approach that directly exploits the fact that the NMR spectrum consists of a low number of Lorentzian peaks. Thus, LPMP is based on the sparsity concept, which is more relevant in the case of NMR than classical CS methods.

OMP Algorithm

Let us fix the notation. For m,n ∈ ℤ and m < n, we write [m, n] = {m, m + 1,…, n}. Let m, n ∈ ℕ. For matrix ∈ M(ℂ), its range and kernel are denoted ran and ker , respectively. The Hermitian conjugate of is denoted by *. The number of non-zero coordinates of ∈ C is denoted by ‖‖0, and supp (the support of ) denotes the set of non-zero coordinates. Note that ‖‖0 = # supp (the cardinality of a set I is denoted by # I). We define χ ⊂ ℂ:

OMP algorithm arose in the context of the sparse approximation problem; see [3]. In order to formulate the problem, let us fix a matrix ⊂ M(ℂ) satisfying ran = ℂ. The matrix is in this context referred to as the dictionary matrix.

Definition 1

Adopting the above notation, we say that ∈ χ is an optimal l-sparse approximate representation of a vector ∈ ℂ if:

We say that y is l-representable, if there exists a vector ∈ χ, such that = .

For a subset I ⊂ [1, n] of cardinality l, we define a vector subspace ℂ ⊂ ℂ:

The restriction of to ℂ is denoted by : ℂ → ℂ.

Definition 2

We say that ⊂ M(ℂ) is l-injective if ker = {0} for all I ⊂ [1, n], # I = l.

For an l-injective matrix and I ⊂ [1, n], # I = l, we define:

Since x is a vector satisfying A = , where is an orthogonal projection of onto ran , we may see the uniqueness of for l-injective . Let us note that ∈ χ is an optimal l-sparse approximate representation of ∈ ℂ if:

Although l-injectivity of the matrix does not exclude a vector ′ ∈ χ, ′ ≠ , such that ‖ — ‖2 = ‖ — ‖2, we shall usually choose an optimal l-sparse approximate and denote it by .

The OMP algorithm provides a certain approximation of x ∈ ℂn, which we denote by Remarkably, is found on the basis of Denoting the support by , the index i ∈ [1, n] is given by:

The OMP vector is defined as ; see Equation (3) for the definition of . For the coherence, restricted isometry property (RIP) conditions that imply the equality , we refer to [3] and [20], respectively. For the recovery limits of the OMP, see [21]. The OMP algorithm (Algorithm 1) is schematically described below.

Input:

measurement matrix A ∈ M(ℂ)

measurement vector ∈ ℂ

accuracy parameter ε > 0

Output:

x ∈ ℂ

Initialization:

I = ∅, = 0

The main loop:

while ‖ — ‖2 ≥ ε do

end while

Free Induction Decay Signals in NMR

The free induction decay signal in NMR consists of the combination of exponentially decaying oscillations (see Equations (1) and (2)). The signal can be multidimensional, but for the sake of simplicity, we will limit ourselves to the one-dimensional indirect part of the two-dimensional signal, corresponding to f(t) from Equation (1). Given the damping factor γ ≥ 0 and the oscillation frequency ω ∈ (x0211D), a single decaying oscillation is given by:

Its Fourier transform: has the form of complex Lorentzian peak:

The FID signal is sampled at the rate given by the Nyquist theorem for the range cut off ν ∈ [–Ω, Ω]. In order to perform the NMR signal processing, the frequency resolution ΔΩ is introduced, where and N ∈ ℕ. For any l ∈ [–N,N], we get the basic complex oscillation e, where . The NMR signals can be approximated by the finite combination of the basic complex oscillations:

The approximation being -periodic works only for . Assuming that the series evaluated at coincides with sampling points of f, i.e., f(t) for k ∈ [0,2N], we get:

Introducing the time unit and the frequency unit , we may enumerate the time and frequency by the integers [0, 2N] and [–N, N], respectively. Then, the relation between the peak intensities and the time samples is given by = , where is the Fourier transform matrix:

l ∈ [–N, N] and k ∈ [0, 2N]. In order to model the Lorentzian peaks within our framework, let us introduce the width peak unit Γ and a damping matrix ∈ M(2+1)×(2+1)(ℂ): for k, k′ ∈ [0, 2N]. Then, the Lorentzian peak of width Γ centered at l is represented by l : where is the l-th column of . Similarly, the Lorentzian peak of the width jΓ centered at l is represented by:

Using the above notation, the FID signal is given by: where ∈ ℂ.

Lorentzian Peak Matching Pursuit

The spectrum corresponding to FID signal (Equation (5)) is given by:

We shall assume that j ∈ [0, J], which gives the maximal width JΓ and the width resolution Γ. Let K = {k1, k2,…, k} ⊂ [0, 2N] be the sampling scheme, ∈ ℂ the measured signal = (t) and = F the matrix consisting of rows of the Fourier matrix enumerated by the elements of K.

The Lorentzian peak matching pursuit algorithm matches a Lorentzian peak in each of its iterations. In order to match a peak LPMP establishes first the peak's center l. The center l1 ∈ [–N, N] of the first peak is determined by:

In order to determine the width of the first peak j1 ∈ [0, J], we find l1, ∈ ℂ, such that:

Then, j1 ∈ [0, J] is a number satisfying: and we put .

In the second step of the algorithm, we determine the second peak's center l2 by:

In order to determine the second peak's width j2 ∈ [0, J], we find 2, 1,1 ∈ ℂ minimizing: where z1,z2 ∈ ℂ Defining j2 E [0, J] as a number satisfying: we put .

We begin the third step of the algorithm by determining the third peak's center l3:

Then, j3 and l3 are found as in the previous steps.

Note that in the k-th step of the algorithm, the peak's width jΓ is found and not changed afterwards. The b, coordinate may, in turn, vary in the n-th step of the algorithm for n > k. The algorithm is performed until the accuracy threshold ‖ – ‖ < ε is achieved; for a different stopping criterion, see Section 5.1.

Stopping Criterion

If a number of Lorentzian peaks in the spectrum cannot be assumed a priori, then the stopping criterion for the LPMP becomes a crucial issue.

Notably, the dimension of a measurement vector y provides a mathematical bound for the number of iterations in LPMP. The implementation of the LPMP algorithm uses an inversion of a matrix, which ceases to be invertible when the number of peaks exceeds the number of measurements. In what follows, we shall introduce the stopping criteria that further restrict the number of iterations. One commonly used criterion is the threshold ε of accuracy used in Algorithm 2.

Let us denote the result of the k-th LPMP iteration by f̂. The k + 1 iteration is performed if: where α ∈ [0,1] is an user-defined number. The criterion reflects the idea that the LPMP is continued as long as its actual iteration increases the explanation of the measured data by the factor α with respect to the previous iteration. In our simulations described in the next section, we put α = 0.98.

Alternatively, the noise amplitude σ can be used for a stopping criterion. In this case, LPMP stops on the k-th iteration if |^| ≤ σ.

Input:

measurement matrix ∈ M×2+1(ℂ)

measurement vector ∈ ℂ

accuracy ε > 0

Lorentzian peaks where j ∈ [0, J], l ∈ [−N, N]

Output:

^∈ ℂ2+1

Initialization:

^= 0, = 0, = 0, i = 0

The main loop:

while ‖ − ^‖2 ≥ ε do

i = i + 1

end while

Mask

NMR experiments are often performed in series, with some or all peaks preserving their positions in at least some of the dimensions. This fact has been exploited to restrict the allowed peak frequencies and to improve the results of NUS (non-uniform sampling) reconstruction, e.g., in hyperdimensional spectroscopy [22], the SIFT method [23] or high-dimensional spectra [24]. It is easy to use this information in the LPMP algorithm by introducing the subset Mask ⊂ [−N, N] of admissible peak centers in the frequency domain. In what follows, we give a detailed description of the masked LPMP algorithm (Algorithm 3), taking into an account the stopping criterion described in Section 5.1.

Input:

measurement matrix ∈ M×2+1(ℂ)

measurement vector ∈ ℂ

stopping parameter α > 0

Lorentzian peaks where j ∈ [0, J], l ∈ [−N, N]

Mask ⊂ [−N, N]

Output:

^ ∈ ℂ2+1

Initialization:

^ = 0, = 0, = 0,i = 0,δ = ‖‖2

The main loop:

while ‖ – ^‖2 ≤ αδ or i = 0 do

δ = ‖ – ^‖2

i = i + 1

end while

Results and Discussion

In this section, we present the results of the LPMP performance tests. Extensive simulations on different levels of sampling sparsity provided a quantitative characterization of the method. As an example of an application, we have chosen the challenging NOESY (nuclear Overhauser effect spectroscopy) experiment [25] with a high dynamic range of peak intensities.

The parameters of the synthetic spectrum used in the simulations are given in Table 1. The spectral parameters were chosen to simulate various experimental scenarios, i.e., (I) peak widths differ from peak to peak, which corresponds to the differences in relaxation rates in the case of the NMR experiment; (II) the range of intensities is high; (III) two of the peaks overlap; and (IV) the peaks' centers and widths are off-grid.

Table 1.

Peak coordinates for the simulation presented in Figure 1.

Center	5.6	22.7	67.8	70.4	194.5	239.25
width	5.2	5.4	2.6	2.8	4	1.5

An example of the LPMP recovery of an NMR spectrum based on 120 sampling points drawn at random from a 256-point full grid is presented in Figure 1. The signal was contaminated with 5% of noise. We performed systemic simulations with a varying sampling level. For each level, 50 different sampling schedules were generated. Denoting the frequency spectrum vector obtained from the full sampling by and the LPMP recovery by , we compute the recovery error ‖–‖, and then, we average it over 50 sampling distributions. The LPMP performance is depicted in Figure 2 (left), where it is also compared with other popular CS algorithms: IST, OMP, COSAMP and StOMP. For all sampling levels, LPMP is revealed to be superior to other CS algorithms.

Figure 1.

Example of the Lorentzian peak matching pursuit (LPMP) recovery. The size of the spectrum is 256. The number of randomly chosen samples is 120. The result of LPMP recovery (blue line) is compared with the original spectrum (red line).

Figure 2.

Quality of the reconstructions of the spectrum from Figure 1. (Left) The comparison of the performance of various algorithms: LPMP (solid line), iterative soft thresholding (IST) (dashed line), stage-wise orthogonal matching pursuit (StOMP) (dotted line), compressive sampling matching (COSAMP) (stared line) and orthogonal matching pursuit (OMP) (squared line). (Right) The performance of the masked LPMP algorithm (solid line) compared with the performance of non-masked LPMP algorithm (dashed line). The absolute recovery errors (in intensity units) are plotted vs. the number of sampling points.

In order to test the masked version of LPMP algorithm, we introduced the mask of ±5 frequency points around the peaks' centers. The simulation results are presented in Figure 2 (right). As expected, masking has a more pronounced effect on the low sampling levels.

The qualitative aspects of LPMP reconstruction can be evaluated on the data from the 2D NOESY experiment. Each reconstructed 1D slice of the 2D NOESY spectrum contains one dominating peak accompanied by a group of smaller peaks. Such a case of large intensity range is very challenging for CS methods. In Figure 3, we present the results of the LPMP spectrum recovery obtained on the basis of 200, 225 and 250 random measurements of the corresponding FID signal of length 512. The 2D NOESY spectrum was recorded on a 600-MHz Varian UNITY Inova spectrometer at 25 °C using a sample of human ubiquitin protein (1 mM concentration).

Figure 3.

1D indirect dimension slice from 2D NOESY spectrum: LPMP-reconstruction from 200, 225 and 250 samples out of a 512-point full grid compared with a fully sampled spectrum.

In order to explain the effectiveness of LPMP algorithm in the NMR context, let us emphasize the key difference between the OMP and LPMP methods. OMP tries to obtain the frequency-domain representation having the minimal number of non-zeros, which cannot be exact, even for a single Lorentzian peak. LPMP improves OMP by replacing every determined spectral non-zero with the Lorentzian peak of the appropriate width. Roughly speaking, this corresponds to the concept of sparsity in the sense of a number of Lorentzian peaks.

Conclusions

Our results show that LPMP is superior to the generic CS methods when applied to exponentially decaying signals. Comparison with OMP, StOMP, IST and COSAMP revealed that the assumption of the Lorentzian spectral line shape improves the quality of the reconstruction with LPMP. Moreover, its performance can be further improved by limiting (“masking”) the support of a signal. We believe that this method will find applications in NMR spectroscopy. Other fields of possible application may include different kinds of spectroscopy, wireless communications, sonar and radar. In fact, reconstruction of any non-uniformly sampled exponentially-damped signal can benefit from the application of LPMP.