Filtering techniques for efficient inversion of two-dimensional Nuclear Magnetic Resonance data

The inversion of two-dimensional Nuclear Magnetic Resonance (NMR) data requires the solution of a first kind Fredholm integral equation with a two-dimensional tensor product kernel and lower bound constraints. For the solution of this ill-posed inverse problem, the recently presented 2DUPEN algorithm [V. Bortolotti et al., Inverse Problems, 33(1), 2016] uses multiparameter Tikhonov regularization with automatic choice of the regularization parameters. In this work, I2DUPEN, an improved version of 2DUPEN that implements Mean Windowing and Singular Value Decomposition filters, is deeply tested. The reconstruction problem with filtered data is formulated as a compressed weighted least squares problem with multi-parameter Tikhonov regularization. Results on synthetic and real 2D NMR data are presented with the main purpose to deeper analyze the separate and combined effects of these filtering techniques on the reconstructed 2D distribution.


Introduction
Nuclear Magnetic Resonance (NMR) of 1 H nuclei of hydrogenous fluids confined in porous media provides a set of techniques exploited to study pore space structure and fluids behavior in a wide range of systems from sedimentary rocks to biological materials [1]. Parameters like the longitudinal (T 1 ) and the transverse (T 2 ) relaxation times, as well as the self-diffusion molecular coefficient, can be determined and associated to properties of the fluids and porous media. When a hydrogenous fluid saturates a porous medium, T 1 and T 2 follow distributions of values and the problem is to correlate these two-dimensional distributions of relaxation times with the two-dimensional experimental data. The resulting two dimensional distributions are shown as two-dimensional relaxation maps characterized by peaks of different size and shape over flat regions. The problem of determining the realxation times distribution from NMR data is a wellknown ill-posed inverse problem. In [2], the 2DUPEN algorithm has been presented, employing multi-parameter Tikhonov regularization and lower bound constraints, for the inversion of twodimensional NMR relaxation data. In order to compute local values of the regularization parameter, 2DUPEN uses the Uniform PENalty (UPEN) principle [3] stating that the product of the regularization parameter and the curvature value in each point of the relaxation times distribution should be constant. Therefore, 2DUPEN is an iterative algorithm where, at each iteration, the locally adapted regularization parameters are automatically updated using the UPEN principle, and an approximate distribution is obtained by solving a L 2 -regularized least squares problem subject to lower bound constraints, usually non negativity constraints. The Newton Projection method [4] is used for the solution of the constrained subproblem. The results of numerical experiments on synthetic and real data show that 2DUPEN is effective in reconstructing peaks and flat regions of the unknown distribution with the same accuracy. However, when dealing with real 2D NMR relaxation data, 2DUPEN may require a huge computational work since it involves solving a sequence of large scale non-negatively constrained least squares problems. Therefore it is necessary to reduce the computation effort both from the point of view of the storage requirement and the computation time. In this work, this aim is obtained by introducing and combining two classical techniques, namely the Singular Value Decomposition (SVD) and the Mean Windowing (MW) methods. The application of SVD for data compression has a long history in inverse problems and also in NMR data inversion. In [5] Sezginer used the SVD method to compress the echo train data and kernel matrix. In [6] echo train data and kernel matrix compression are extended to inverse problems with tensor product structure. In that work a preliminary step applies SVD to reduce problem size by projecting the data and kernels onto a lower dimensional subspace. The main difference between [6] and our method is the employing of the uniform penalty principle to compute the locally adapted regularization parameters whereas in [6] a scalar global regularization parameter is used. Hence two different constrained least squares models are solved. Moreover the projected Newton method is used in the present work to solve the constrained minimization problem where in [6] the BRD [7] method is used. The MW method has been used in [3] to compress the data blocks considerably reducing the amount of data. In [8], a preliminary use of MW and SVD filters have been successfully tested to improve the computational cost of 2DUPEN when applied to real 2D NMR data. This enhanced version of 2DUPEN was called I2DUPEN. The MW and SVD filters turned out to be effective in reducing the problem size and the computational effort of 2DUPEN with real data, thus motivating a deeper analysis of these filters on simulated data. Hence, the main purpose of this paper is to analyse the separate and combined effects of these filtering techniques on synthetic and real 2D NMR data. Our experimental analysis shows that the combined use of the MW and SVD filters has a considerable effect on reducing the computational cost of 2DUPEN while preserving the reconstruction quality. In particular I2DUPEN preserves the principal feature of UPEN algorithm, that is to reproduce the underlying distribution in presence of peaks and flat regions. The sequel is organized as follows: in Section 2 the problem of the inversion of 2D NMR data is described. Section 3 presents the filtering techniques used in the 2DUPEN algorithm and numerical results are reported in Section 4. Finally, conclusions are given in Section 5.

Problem description
If we consider 2D NMR relaxation data acquired using a conventional Inversion-Recovery (IR) or Saturation-Recovery (SR) experiment detected by a Carr-Purcell-Meiboom-Gill (CPMG) pulse train [9], then the evolution time t 1 in IR or SR and the evolution time t 2 in CPMG are two independent variables and the 2D NMR relaxation data S(t 1 , t 2 ) can be expressed as: where the unknown F (T 1 , T 2 ) is the distribution of T 1 and T 2 relaxation times, the kernels have the following expressions: and e(t 1 , t 2 ) represents Gaussian additive noise. The unknown distribution F (T 1 , T 2 ) is known to be greater or equal to a constant ρ ∈ R; in this work, ρ = 0. After considering M 1 × M 2 samples of the times t 1 and t 2 , and N 1 × N 2 samples of the relaxation times T 1 and T 2 , problem (1) can be discretized as: where is the vector reordering of the distribution and e ∈ R M is the vector with the discretized noise. In order to avoid artifact peaks in the computed distribution, 2DUPEN uses multiple-parameter Tikhonov regularization and solves the minimization problem where · is the L 2 norm, L ∈ R N ×N is the discrete Laplacian operator and λ i are the regularization parameters, i = 1, . . . , N . The 2DUPEN method is an iterative procedure where, at each iteration, suitable values for the λ i 's are determined by imposing that all the non-null products λ i (Lf ) 2 i have the same constant value (UPEN principle [2,3]) and an approximated distribution is obtained by solving a problem of the form (3) with the Newton Projection method. 2DUPEN is able to achieve high quality reconstructions [2] but our experience indicates that considerable computational effort is needed. The computational effort is given by the cost of multiplications between the matrix K 2 ⊗ K 1 (and its transpose) and vectors in the Newton Projection method (see [2] for details about the algorithm). Using the properties of the Kronecker products, the matrices K 2 ⊗K 1 and (K 2 ⊗K 1 ) T are never formed and the matrix vector products are computed as: is the vector obtained by columnwise reordering the matrix V; therefore, in our case, x = vec(X) and X ∈ R N 1 ×N 2 . The matrix vector multiplications involving (K 2 ⊗ K 1 ) T are performed using the relation ( . Therefore an effective reduction of computation cost and storage requirements is obtained, in this work, by means of techniques that allow us to decrease the problems size while preserving properties of the Kronecker product matrices. For this reason we do not propose here methods, such as parallel block iterative algorithms, which require the formation of the matrix K 2 ⊗ K 1 , with storage requirement of M 1 · M 2 · N 1 · N 2 bytes. Moreover the formation of K 2 ⊗ K 1 would require the parallelization of matrix vector products involving major changes in the software and hardware devices usually available in real NMR data analysis applications.

Filtering techniques
In order to accelerate the implementation of the 2DUPEN method we explore projection methods (SVD filtering) and data reduction methods (MW technique) sparately and then we couple them together. In the following paragraphs we describe the different techniques and deeply analyze the Improved 2DUPEN algorithm (I2DUPEN) which includes both SVD and MW filters.
The Mean Windowing-like filter. The CPMG data blocks of an IR-CPMG or SR-CPMG sequence, used to acquire the two dimensional T 1 −T 2 experimental data, usually have thousands of points; in order to reduce their number, the CPMG data blocks are averaged on windows of logarithmically increasing size, decreasing the number of points by orders of magnitude. In order to take into account the number of points averaged in each window [3], weighting factors are introduced in the data fitting term of (3). Since each IR or SR point has the same weight, the weighting matrix can be expressed as a Kronecker product of diagonal matrices, i.e. (B ⊗ I) where B contains the numbers of points of each CPMG windowed block, and I is the identity matrix. Assuming that the size of the whole 2D data is M 1 × M 2 , after the application of the windowing, the vector reordered data iss ∈ R M 1 ·M 2 with M 2 M 2 . Introducing the weighting matrix (B ⊗ I), we obtain the following weighted least squares problem: where K 2 ∈ R M 2 ×N 2 is the t 2 kernel of the reduced size problem, B ∈ R M 2 ×M 2 . We have: Using (5) and the Kronecker product property (B ⊗ I)(K 2 ⊗ K 1 ) = (BK 2 ⊗ K 1 ), problem (4) can be written as: min where K 2 = BK 2 and s = (B ⊗ I)s.
The SVD filter. While the MW filter reduces the data size by replacing the points in a window with their average value, with the SVD filter it is possible to reduce the data dimensions by projecting s onto a lower-dimensional subspace spanned by the first left singular vectors of K 1 and K 2 .
with U i , V i orthogonal matrices and Σ i diagonal matrices of the singular values, i = 1, 2. Hence, the data fitting term of (3) can be written as: We use the Truncated SVD (TSVD) and consider only the SVD components corresponding to the largest singular values. Let τ be a positive parameter and let ρ 1 , ρ 2 be the number of singular values of Σ 1 , Σ 2 greater or equal to τ , we replace problem (3) with where: Filtered 2DUPEN methods. The MW and SVD filters can be used separately or together obtaining the following modifications of the 2DUPEN algorithm. If the SVD filter is applied to problem (3) without MW data reduction, then we obtain a modified 2DUPEN method, named 2DUPENm 1 , which applies 2DUPEN to solve problem (8). If the SVD filter is applied, the data are reduced by MW but the weights are not used in the data fitting term, we have a modified 2DUPEN method, named 2DUPENm 2 , which uses 2DUPEN to solve a SVD filtered problem of the form (8) whereK 2 is obtained by TSVD of the matrix K 2 andŝ = (Û T 2 ⊗Û T 1 )s. If the data are reduced by the MW filter, the weights are used in the data fitting term and no SVD filter is used, we obtain a modified 2DUPEN method, named 2DUPENm 3 , which applies 2DUPEN to solve (6). Finally, the so-called Improved 2DUPEN (I2DUPEN) algorithm uses both the MW and SVD filters and solves a problem of the form (8) with weighting factors, whereK 2 is obtained by TSVD of the matrix K 2 andŝ = (Û T 2 ⊗Û T 1 ) s.  [2] for detailed explanation), the I2DUPEN algorithm is stated as follows: (iii) Compute the TSVD of K 1 and K 2 and setK 1 =Σ 1V by applying a few iterations of the Gradient Projection method [10].
where β 0 , β p , β c are 2DUPEN parameters [2]. We observe that the core iteration of 2DUPEN is preserved in I2DUPEN and the filters act only on the kernels and the data. The SVDs of the kernel matrices are performed only once and their computation cost do not affect significantly the global computation cost.

Numerical Results
The effectiveness of the proposed strategies is extensively evaluated on a synthetic test problem in paragraph 4.1 and its behaviour on real data is described in paragraph 4.2.

Experiments on a synthetic test problem
In this paragraph we compare, on a synthetic test problem, the proposed improvement strategies 2DUPENm 1 , 2DUPENm 2 , 2DUPENm 3 and I2DUPEN. The following parameters have been used in all the 2DUPEN modifications: β 0 = 10 −6 , β p = β c = 1, tol = 0.001. We consider here a model distribution f * ∈ R N , N = 64 · 64, and use it to synthesize the simulated IR-CPMG sequence data with 128 × 2048 non uniformly distributed T 1 − T 2 data. The test problem has two Gaussian peaks positioned along the main diagonal over a zero flat area ( figure 1(a)). The noisy data s are obtained by adding Gaussian white noise of level δ to the synthetic data: s = (K 2 ⊗K 1 )f * +e where e = δη and η is a normal random Gaussian vector such that η = 1. Several values of δ have been tested, we report here the results obtained with δ = 10 −2 which is close to the noise level in the real data. The Relative Error (Err) is used to evaluate the reconstruction quality. In table 1 we report the reconstruction errors and the computation time for the 2DUPENm 1 , 2DUPENm 2 and I2DUPEN algorithms for different threshold values τ . The 2DUPEN algorithm solves this problem in 17770.81 s with relative error Err = 1.972 10 −2 ( figure 1(b)). Observing the left picture of figure 2 we see the log-uniform decreasing distribution   of the singular values of K 1 and K 2 yielding to a severe ill conditioning of the problem. At first we introduce only the SVD filter with threshold τ ranging from 0 (i.e. no filter) to 10 −2 (2DUPENm 1 ). In this case, the whole acquisition set of 128 · 2048 data is reduced by projection to 64 · 64 samples when τ = 0 down to a minimum of 12 · 11 samples when τ = 10 −1 , as reported in column 2 of table 1. We observe that the size reduction causes a proportional decrease of the computational time and an increase of the relative error values (columns 3 and 4 of table 1). The computed distribution is shown in figure 1(d). We examine now the results obtained using data reduced by the MW filter. In our example the reduced data have 128 · 118 samples. The singular values distribution does not change significantly for the introduction of the weights (right plot of figure 2) and the problem is still badly ill-conditioned. We see that the MW filter causes an increases of the reconstruction errors (columns 6 and 9 of table 1) even if the computed 2D maps are visually indistinguishable. We observe that I2DUPEN has always the smallest computation times and gives restorations with comparable qualitative accuracy. Finally, 2DUPENm 3 obtains, in 2062.90 s, a computed 2D map with relative error Err = 5.002 · 10 −2 ( figure 1(f) From the different reconstructions and numerical results reported in figure 1 and table 1, we observe that the introduction of a filtering technique causes an increase of the reconstruction errors. We observe that when the full data s ∈ R M are available, 2DUPENm 1 is the best improvement of 2DUPEN both in terms of achieved quality and time reduction for large values of τ . However, if only the windowed datas are available, 2DUPENm 2 is the best way to improve 2DUPEN accuracy, while I2DUPEN is the best improvement of 2DUPEN computational efficiency. Concluding, if the whole data set is available it is not convenient to apply MW because 2DUPENm 1 with large threshold (for example τ 10 −3 ) efficiently computes the most accurate distribution. On the other hand, when the data are available only in MW compressed form, I2DUPEN is more suitable because it is faster than 2DUPENm 2 even if slightly less accurate.

Experiments on real data
In this paragraph we report the results obtained by testing I2DUPEN and its modifications on a real data set. A SR-CPMG test set of moderate size is used in order to be able to test 2DUPEN and all the I2DUPEN variants. A sample (5 × 5 × 2 cm 3 ) of biocalcarenite with porosity about 47%, was saturated under high vacuum with distilled water and promptly sealed into glass containers in order to avoid the evaporation of water. NMR measurements were performed at 25 by a NMR single-sided apparatus based on a Kea II NMR console and a NMR-MOUSE PM25 (Magritek, GmbH), controlled by the software Prospa (Magritek, GmbH). The 2D measurement, longitudinal-transverse relaxation curve (T 1 − T 2 ) was acquired by a SR-CPMG pulse sequence. The pulse sequence is characterized by an application of a train of π/2 pulses, in order to pre-saturate the transverse magnetization, and an encoding-reading stage The T 1 relaxation signal was acquired with 32 values of τ chosen in geometrical progression from 1 ms up to 400 ms, with N E = 1024 (number of acquired echos, echo times T E = 60µs) on each CPMG, and number of scans equal to 24. All curves were acquired using phase-cycling procedures. The following parameters have been used in all the 2DUPEN modifications: β 0 = 10 −5 , β p = 10, β c = 30, tol = 0.05. The reconstruction with N 1 × N 2 = 100 × 100 points has been computed by 2DUPEN in 72.55 s. The 3D surface and contour map are reported in figure 3: a peak hight 1.25 · 10 −5 a.u. (arbitrary unit) is located at T 2 = 53 ms and T 1 = 135 ms. Small peaks closed to the bottom border, probably are related to some small artefacts on the acquired data. It is noteworthy to highlight that the UPEN algorithm it isn't only a powerful tool to invert data, but also a tool to evaluate the quality of acquired data. In this case the data reduced by MW  peak. The ρ 1 , ρ 2 values are the same for all the methods. As for synthetic data, 2DUPENm 2 is the less performing algorithm since always gives an inaccurate localization of the peak even if the computation times are comparable to those of I2DUPEN. We remark that the correct peaks localization is a very important information to be recovered in real applications. Concerning 2DUPENm 1 it requires the highest computation times and obtains a quite precise peaks localization. Hence we can conclude that I2DUPEN is the best improvement of the 2DUPEN algorithm for accuracy and computation times.

Conclusions
In several NMR applications, users are interested in efficient methods able to reconstruct relaxation time distributions where the peaks are correctly localized and the underlying 3D distribution volume is well approximated. Our experimental analysis shows that I2DUPEN is the best modification of 2DUPEN in order to efficiently obtain correct information about the unknown relaxation times distributions. The loss of information due to the MW filter is less critical on real data because of unavoidable instrumental limitations and artifacts. This explains the better accuracy of I2DUPEN compared to I2DUPENm 1 in case of real data.