Resolution and uncertainty of Laplace inversion spectrum

Magnetic Resonance Imaging 25 (2007) 445 – 448

Resolution and uncertainty of Laplace inversion spectrum Yi-Qiao Song Schlumberger-Doll Research, Ridgefield, CT 06877, USA

Abstract Laplace inversion has been widely used in analyzing nuclear magnetic resonance data to obtain a spectrum of decay time constant, such as the T 2 spectrum. However, due to the ill-conditioned nature of such inversion, it is difficult to determine the reliability and error of the resulting spectrum. This article describes methods used to analyze resolution and to obtain bounds of integral quantities of the spectrum. D 2007 Elsevier Inc. All rights reserved. Keywords: Laplace inversion; Relaxation; Diffusion; Resolution; Error

1. Introduction In the study of heterogeneous materials (both natural substances and manmade products), nuclear magnetic resonance (NMR) spin relaxation spectra and diffusion constants are often used as fingerprints of molecular species, structure, and dynamics. For example, water and crude oil present in oil reservoirs can be distinguished by diffusion and relaxation experiments [1,2]. Typically, spin relaxation and diffusion are manifested as decaying signals. Data analysis often involves Laplace inversion to obtain a spectrum (or distribution) of relaxation times or diffusion constants. The inversion is ill-conditioned in the sense that a small noise in the data can cause large changes in the spectrum. Despite the ill-conditioned nature, there are wellestablished mathematical methods (e.g., Tikhonov regularization [3] and maximum entropy method [4,5]) used to handle this and other types of inversion [6,7]. Numerical recipes [8] and NMR literature [9–12] are a good introduction to this problem. One of the critical issues regarding such inversion is the resolution of resulting spectra. For example, given a spectrum with a single peak, one may ask whether spectral width is determined by finite signal-to-noise ratio (SNR) or whether it reflects the true width of an underlying phenomenon. We will analyze the salient features of numerical Laplace inversion and provide a simple numerical method to estimate the resolution of the spectrum.

E-mail address: [email protected]. 0730-725X/$ – see front matter D 2007 Elsevier Inc. All rights reserved. doi:10.1016/j.mri.2006.11.023

Second, an inversion spectrum is often used to further calculate other quantities, such as porosity, partial porosity or other quantities averaged over the spectrum. We develop a method to determine the upper and lower bounds of such quantities. This method does not limit itself to only regularized solutions and does not require prior assumptions about the spectrum.

2. Spectral resolution For two spectra with a narrow peak and a broader peak, respectively, time domain data can be quite similar and difficult to distinguish with a finite amount of noise. We analyze this behavior systematically to obtain the sharpest spectral basic function to describe the data. The wavelength of such basic function determines the best possible resolution for the given data and noise. 2.1. General mathematical treatment For simplicity, we demonstrate the idea of estimating resolution on T 2 relaxation data, written in the discretized form: M ¼ KF þ :

ð1Þ

Here M is the data vector of length n, e is the noise and F is an unknown spectrum vector discretized at N points. The kernel matrix K is of dimensions n N, K ij uexp(s i /T j ). Typically, the relaxation time T is expected to vary from milliseconds to seconds. Thus, in Eq. (1), the variable T is often logarithmically spaced, and the estimated density F

446

Y.-Q. Song / Magnetic Resonance Imaging 25 (2007) 445–448

corresponds to the spectrum on a logarithmic scale. We define kulog10T. Using singular value decomposition (SVD), the kernel matrix K can be given as: K ¼ U RV

T

ð2Þ

where R is a diagonal matrix with singular values (s i ) in descending order, and U and V are unitary matrices. The dimensions of U and V are n n s and N n s, respectively, where n s is the number of nonzero singular values. V T is the matrix transpose of V. The column vectors of U denoted by u i are referred to as left-hand singular vectors and are orthogonal to each other, forming a subregion within a data space (a vector space of dimension n). Similarly, the column vectors of V, denoted by v i , are referred to as right-hand singular vectors and form a subspace in a spectrum space. When some of the singular values are zero or close to zero, the matrix K is singular. In this case, there exist nonzero vectors that satisfy KF = 0, and the amplitude of these vectors cannot be determined from the measured data. The least squares solution to this problem (Eq. (1)) has an analytical form F = VA1U TM. The spectrum is a linear combination of the singular vectors v i : F¼

ns X

bi vi ;

i¼1

bi ¼

uTi M : si

ð3Þ

idea of restricting contribution from small singular values to avoid noise amplification. Eq. (6) is appealing for its simplicity in concept and ease of numerical analysis. 2.2. Resolution Eq. (4) shows that F is a linear combination of v i . When higher singular values and functions contribute to the estimated spectrum, the resulting spectrum can be sharper. Thus, a natural definition of the resolution (Dk) of a reconstructed spectrum can be given by the bwavelengthQ of v r of the smallest singular value chosen to estimate the spectrum. Singular vectors are oscillatory in nature. Thus, wavelength can be loosely defined in terms of the number of zero crossings of v r : Dkf

k2  k1 : r

ð7Þ

Here, E1 and E2 are the lower and upper limits of log10T, respectively. Fig. 1 shows the effects of such truncation and the resulting resolution. For a test spectrum (e.g., solid line, Fig. 1A), the corresponding time domain data M is calculated. Then, Eq. (6) was used to determine r, and Eq. (4) was used to reconstruct the spectrum, shown as a dotted line. The reconstructed spectrum is significantly wider than the original test spectrum due to truncated SVD. This

Eq. (3) shows that singular vectors corresponding to smaller singular values contribute heavily to the estimated spectrum. Thus, any noise recorded in the data M is amplified by small singular values. A natural approach to resolving the noise amplification problem is to limit SVD to only the first r singular values: F¼

r X

bi vi :

ð4Þ

i¼1

The maximum r can be determined using the time domain data and the known noise. Given the singular vectors u i : M

r X

ci ui þ ;

ð5Þ

i¼1

where c i is the weight of u i , c i =uTi M. When one multiplies the above equation with u rT, one obtains u rTM6c r +u rTe. Thus, the minimum singular value (the largest r) is reached when noise contribution is comparable with c r,|u rTM|~|u rTe|=r, where r is the noise variance. In other words, we select all singular values that satisfy: uTr M z r:

ð6Þ

In other more sophisticated approaches such as Tikhonov regularization [3] and Maximum entropy method [4], the truncation of SVD is gradual, instead of the sharp cutoff of the current method. However, they are similar in the basic

Fig. 1. Plots of test spectra and the corresponding reconstructed spectrum with truncated SVD. All simulated data use 1000 echoes and an echo spacing of 1 ms. (A) One test spectrum (solid line) and its corresponding reconstruction (dashed line) are shown. (B) Five test spectra (solid lines) and their corresponding reconstruction (dashed lines) are shown, demonstrating the broadening effect for different T 2 values. The SNR is 1000. (C) The same types of plot as (B), with an SNR of 100.

Y.-Q. Song / Magnetic Resonance Imaging 25 (2007) 445–448

447

both undesirable qualities. The two properties, misfit and complexity, are in conflict and cannot be made arbitrarily small simultaneously; thus, by a suitable choice of a, a satisfactory compromise is sought to provide a plausible density function. In fact, there are many alternative solutions that satisfy Eq. (1).

Fig. 2. Data misfit as a function of x. The shaded region corresponds to all functions F obeying F N 0 and w TF = x. The hatched area corresponds to those solutions that further satisfy data equation Eq. (1).

broadening (i.e., resolution) is determined by SNR and the structure of kernel function. Since the singular vector (v i ) may not have a uniform pitch along the T 2 axis, the width of reconstructed spectra may be different. Fig. 1B shows the resolution at different T 2 values. Reconstructed spectra show similar widths for test spectra below T 2 b 101 s; the width is noticeably broader at T 2~100 s and much broader at T 2~4 s. This behavior is intuitively reasonable since the echo train only extends to 1 s, and there are insufficient data to define a sharp spectrum beyond 1 s. Fig. 1C further demonstrates the effect of reducing SNR on the reconstructed spectrum as compared with those in Fig. 1B. As SNR is reduced, the maximum number of SVD is reduced (Eq. (6)); thus, the resolution is further reduced. Such reduction of resolution is often observed in the results of regularized Laplace inversion algorithm. Since these algorithms employ the same principle for SVD truncation (albeit more sophisticated) as the one demonstrated here, the width of the resulting spectrum of such inversion is consistent with the resolution derived in this work. Such consistency has been demonstrated in numerical simulations [13]. It is important to point out that regularized solutions only represent an exceedingly small fraction of all solutions F satisfying Eq. (1). The concept of resolution described in this work means that, for a given data set, it is impossible to resolve fine spectral features within the resolution limit. It does not imply that the reconstructed spectrum (such as the dashed lines in Fig. 1) are true solutions. 3. Uncertainty The standard solution of Eq. (1) uses the technique of regularization, which seeks to minimize: jjM  KFjj2 þ ajjFjj2 ;

ð8Þ

with some suitable choice of a. The rationale for the method is that in Eq. (8), the first term represents a measure of model misfit and the second term represents complexity—

Fig. 3. Calculation of bounds on porosity for given T 2 decay data with SNR = 20. (A) The X 1(/) curve. The best fit is shown as an open square. (B) Plot of several spectra: the original model spectrum with two d function peaks, the spectrum of the best-fit solution and the spectra at the upper and lower bounds. (C) The corresponding time domain data for the best-fit, upper and lower bounds.

448

Y.-Q. Song / Magnetic Resonance Imaging 25 (2007) 445–448

An alternative to regularization that does provide uncertainty information is the Bayesian formulation [7], but this approach uses statistical characterization of density function and assumption of a prior distribution that we wish to avoid. The associated computational effort can also be quite extensive. We investigate the properties of the set of all models matching the measurements to a specified degree, a class much wider than that of traditional regularized models. For example, since the resulting spectrum is often used to further calculate other quantities for various applications, we will consider the uncertainty of such quantities for all possible solutions to be constrained only by the data. To illustrate the idea, we will consider a simple linear function of F such as: Iw ða; bÞ ¼

Z

b

wðT ÞF ðT ÞdT

ð9Þ

bounds are shown in Fig. 3B together with the original model spectrum, illustrating the potential range of the spectrum at the extrema of porosity estimates.

4. Conclusion We have outlined easily computed numerical methods to estimate the spectral resolution of a Laplace inversion spectrum and the uncertainty of the integral properties of the spectrum. They could be useful in determining the reliability of some of the features in a spectrum, such as the presence of double peaks, in estimating the error bars of quantities derived from such experiments. These methods can also be used as a guide to designing diffusion and relaxation experiments aiming to maximize information. Although we use NMR data as examples, the method and discussion are certainly applicable to data analysis beyond NMR.

a

for finite integration limits (a,b) and bounded functions w z 0. The general properties of the solutions are illustrated in Fig. 2, where all functions (in F space) are located in the half space (shaded area) above a concave curve. The figure is defined by v 2u||MKF||2 and x = w TF. w is the discretized weighting function. For a fixed value of x, functions satisfying x =w TF exist, and they exhibit a range of v 2 with a minimum. If the noise variance of the data is known, such as the horizontal line at v 2 = X 02, functions within a hatched area are the only solutions satisfying Eq. (1). This plot illustrates that one of the properties of these solutions is that x of them will be within the bounds of x + and x . A numerical procedure was recently developed [14] to directly evaluate x + and x  using a non-negative least squares (NNLS) algorithm. We illustrate the method through an example (Fig. 3) with constraint as the total porosity /. Two thousand echo train data were simulated with t e = 1 ms with an SNR of 20. The model spectrum presents two d function peaks, as shown in Fig. 3A. Specifically, NNLS is used first to obtain the best fit corresponding to the minimum v 2 point in Fig. 3, with the corresponding spectrum F 0. Then, a series of / from 0.8 to 1.3 was used to perform NNLS fitting to obtain the minimum v 2 for a given /, forming the X 1(/) curve. The intersects of the X 1 curve with the desired variance X 02 are the bounds.These calculations are straightforward, using MATLAB with NNLS routines such as lsqlin and lsqnonneg. Fig. 3A shows the X 1 curve, with porosity as the constraint. The spectra at the best-fit lower and upper

References [1] Kleinberg RL, Vinegar HJ. NMR properties of reservoir fluids. Log Anal 1996;37(6):20. [2] Kleinberg R. Well logging. In: Grant DM, Harris RK, editors. Encyclopedia of nuclear magnetic resonance. New York7 John Wiley and Sons; 1995. [3] Tikhonov AN, Arsenin VY. Solutions of ill-posed problems. New York7 John Wiley and Sons; 1977. [4] Skilling J. Classic maximum entropy. In: Skilling J, editor. Maximum entropy and Bayesian methods. Dordrecht7 Kluwer Academic; 1989. p. 45 – 52. [5] Gull SF. Developments in maximum entropy data analysis. In: Skilling J, editor. Maximum entropy and Bayesian methods. Dordrecht7 Kluwer Academic; 1989. p. 53 – 71. [6] Lawson CL, Hanson RJ. Solving least squares problems. Englewood Cliffs (NJ)7 Prentice-Hall; 1974. [7] Tarantola A. Inverse problem theory. Amsterdam (The Netherlands)7 Elsevier; 1987. [8] Press WH, Teukolsky SA, Vetterling WT, Flannery BP. Numerical recipes in C. Cambridge7 Cambridge University Press; 1997. [9] Fordham EJ, Sezginer A, Hall LD. Imaging multiexponential relaxation in the ( y, logT 1) plane, with application to clay filtration in rock cores. J Magn Reson A 1995;113:139. [10] Kroeker RM, Henkelman RM. Analysis of biological NMR relaxation data with continuous distributions of relaxation times. J Magn Reson 1986;69:218. [11] Borgia GC, Brown RJS, Fantazzini P. Uniformpenalty inversion of multiexponential decay data. J Magn Reson 1998;132:65. [12] Borgia GC, Brown RJS, Fantazzini P. Uniform penalty inversion of multiexponential decay data: II. Data spacing, T 2 data, systematic data errors, and diagnostics. J Magn Reson 2000;147:273. [13] Song Y-Q, Venkataramanan L, Burcaw L. Determining the resolution of Laplace inversion spectrum. J Chem Phys 2005;122:104104. [14] Parker R, Song Y-Q. Assigning uncertainties in the inversion of NMR relaxation data. J Magn Reson 2005;174(2):314.