Solving 2D Fredholm Integral from Incomplete Measurements Using Compressive Sensing^*

1. Introduction

We present a method of solving the two-dimensional (2D) Fredholm integral of the first kind from a limited number of measurements. This is particularly useful in the field of nuclear magnetic resonance (NMR), in which making a sufficient number of measurements takes several hours. Our work is an extension of the algorithm in [56] based on the new idea of matrix completion; cf. [10, 29, 52].

A 2D Fredholm integral of the first kind is written as

where k₁ and k₂ are continuous Hilbert–Schmidt kernel functions and f, g ∈ L²(ℝ²); cf. [27]. 2D Fourier, Laplace, and Hankel transforms are all common examples of Fredholm integral equations. Applications of these transformations arise in any number of fields, including methods for solving PDEs [31], image deblurring [4, 38], and moment generating functions [40]. This paper specifically focuses on Laplace-type transforms, where the kernel singular values decay quickly to zero.

To present the main idea of the problem, the data M is measured over sampling times τ₁ and τ₂ and is related to the object of interest ℱ(x, y) by a 2D Fredholm integral of the first kind with a tensor product kernel,

where ε(τ₁, τ₂) is assumed to be Gaussian white noise. In most applications, including NMR, the kernels k₁ and k₂ are explicit functions that are known to be smooth and continuous a priori. Solving a Fredholm integral with smooth kernels is an ill-conditioned problem, since the kernel’s singular values decay quickly to zero [34]. This makes the problem particularly interesting, as small variations in the data can lead to large fluctuations in the solution.

For our purposes, ℱ(x, y) represents the joint probability density function of the variables x and y. Specifically in NMR, x and y can be the measurements of the two combination of the longitudinal relaxation time T1, transverse relaxation time T2, diffusion D, and other dynamic properties. Knowledge of the correlation of these properties of a sample is used to identify its microstructure properties and dynamics [6].

This paper focuses on the discretized version of the 2D Fredholm integral,

where our data is the matrix M ∈ ℝ^N₁×N₂, matrices K₁ ∈ ℝ^N₁×N_x and K₂ ∈ ℝ^N₂×N_y are discretized versions of the smooth kernels k₁ and k₂, and the matrix F ∈ ℝ^N_x×N_y is the discretized version of the probability density function ℱ(x, y) which we are interested in recovering. We also assume that each element of the Gaussian noise matrix E has zero mean and constant variance. And since we have assumed that ℱ(x, y) is a joint probability density function, each element of F is nonnegative.

Venkataramanan, Song, and Hürlimann [56] laid out an efficient strategy for solving this problem given complete knowledge of the data matrix M. The approach centers around finding an intelligent way to solve the Tikhonov regularization problem,

where || · ||_F is the Frobenius norm.

There are three steps to the algorithm in [56] for solving (1.2).

1. Compress the data. Let the SVD of K_i be

   $$K_i=U_i{S}_i{V}_i^{\prime },i\in \{1,2\}.$$

   Because K₁ and K₂ are sampled from smooth functions k₁ and k₂, the singular values decay quickly to 0. Let s₁ be the number of nonzero singular values of K₁, and let s₂ be the number of nonzero singular values of K₂. Then U_i ∈ ℝ^N_i×s_i and S_i ∈ ℝ^s_i×s_i for i = 1, 2, as well as V₁ ∈ ℝ^N_x×s₁ and V₂ ∈ ℝ^N_y×s₂.

   The data matrix M can be projected onto the column space of K₁ and the row space of K₂ by $U_1{U}_1^{\prime }M{U}_2{U}_2^{\prime }$. We denote this as $\stackrel{\sim }{M}=U_1^{\prime }M{U}_2$. The Tikhonov regularization problem (1.2) is now rewritten as

   $$\widehat{F}=\underset{F\ge 0}{argmin}{\Vert U_1\stackrel{\sim }{M}{U}_2^{\prime }-U_1{U}_1^{\prime }{K}_{1}F{K}_2^{\prime }{U}_2{U}_2^{\prime }\Vert }_F^2+{\Vert M\Vert }_F^2-{\Vert U_1\stackrel{\sim }{M}{U}_2^{\prime }\Vert }_F^2+\alpha {\Vert F\Vert }_F^2$$

   (1.3)

   $$=\underset{F\ge 0}{argmin}{\Vert \stackrel{\sim }{M}-(S_1{V}_1^{\prime })F{(S_2{V}_2^{\prime })}^{\prime }\Vert }_F^2+\alpha {\Vert F\Vert }_F^2,$$

   (1.4)

   where (1.4) comes from U₁ and U₂ having orthogonal columns, and the second and third terms in (1.3) being independent of F. The key note here is that M̃ ∈ ℝ^s₁×s₂, which significantly reduces the complexity of the computations.

2. Optimization. For a given value of α, (1.4) has a unique solution due to the second term being quadratic. We shall detail the method of finding this solution in section 4.

3. Choosing α. Once (1.4) has been solved for a specific α, an update for α is chosen based on the characteristics of the solution in step 2. Repeat steps 2 and 3 until convergence. Again, this is detailed in section 4.

The approach in [56] assumes complete knowledge of the data matrix M. However, in applications with NMR, there is a cost associated with collecting all the elements of M, namely, time. With the microstructure-related information contained in the multidimensional diffusion-relaxation correlation spectrum of the biological sample [49, 22, 25, 20, 35, 54] and high-resolution spatial information that magnetic resonance imaging (MRI) techniques can provide, there is a need to combine the multidimensional correlation spectra NMR with 2D/3D MRI for preclinical and clinical applications [21]. Without any acceleration, however, it could take several days to acquire this data.

In practice, the potential pulse sequences for the combined multidimensional diffusion-relaxation MRI would be single spin echo (90°–180° acquisition and spatial localization) with saturation, inversion recovery, driven-equilibrium preparation to measure T1-T2 correlation, and diffusion weighting preparation for D-T2 measurements. With these MRI pulse sequences, a single point in the 2D T1-T2 or D-T2 space is acquired for each “shot,” and the total time for the sampling of the T1-T2 or D-T2 space is determined directly by the number of measurements required to recover F from (1.2). A vastly reduced number of sample points in M, together with rapid MRI acquisition techniques, which can include, e.g., parallel imaging [50], echo planar imaging (EPI) [24], gradient-recalled echo [36], and sparse sampling with compressed sensing [45], could reduce the total experiment time sufficiently to make this promising technique practicable for preclinical and clinical in vivo studies.

Notice that, despite collecting all N₁ × N₂ data points in M, step 1 of the algorithm immediately throws away a large amount of that information, reducing the number of data points to a matrix of size s₁ × s₂. M̃ is effectively a compressed version of the original M, containing the same information in a smaller number of entries. But this raises the question of why all of M must be collected when a large amount of information is immediately thrown away, since we are interested only in M̃.

This question is what motivates the introduction of a compressive sensing–type approach. The task is to undersample signals that are “compressible,” meaning that the signal is sparse in some basis representation [12, 11, 23]. The problem of recovering M falls into a subset of this field known as low-rank matrix completion; see [10, 29, 14].

An n × n matrix X that is rank r requires approximately nr parameters to be completely specified. If r ≪ n, then X is seen as being compressible, as the number of parameters needed to specify it is much less than its n² entries. It is less clear how to recover X from a limited number of coefficients efficiently. But the results of [10] showed that it is possible to recover X from, up to a constant, nr log(n) measurements by employing a simple optimization problem. These findings were inspired, at least in part, by [12, 11]. Also, the types of measurements we utilize in this paper, operator bases with bounded norm, originated via quantum state tomography [30].

Compressive sensing has been used in various forms of medical imaging for several years. The authors of [45] originally proposed speeding up MRI acquisition. The authors of [47] introduced group sparsity into consideration for accelerating T2-weighted MR. Finally, the authors of [55, 39] both introduced basic ideas of compressive sensing into the NMR framework and attained promising results, but the results did not utilize the matrix completion aspect of the physical problem and did not introduce any theoretical guarantees for reconstruction.

This paper develops an alternative to [56] which incorporates matrix completion in order to recover M̃ from significantly fewer measurements. It is organized as follows. Section 2 examines how recovery of M̃ fits into existing theory and shows that data from the 2D Fredholm integral can be recovered from 10% of the measurements. Section 3 covers the practical considerations of the problem and discusses the error created by our reconstruction. Section 4 covers the algorithms used to solve the low-rank minimization problem and invert the 2D Fredholm integral to obtain F. Section 5 shows the effectiveness of this reconstruction on simulated data. Appendix A contains a detailed proof of the central theorem of this paper.

2. Data recovery using matrix completion

3. Numerical considerations

Section 2 gives theoretical guarantees about the error of estimating M̃₀ with the recovered M̃. We shall address several issues related to practical applications in this section. We shall let M̃₀ be the original compressed data matrix we are hoping to recover, and let M̃ be the approximation obtained by solving (2.8) for the sampling operator ℛ_Ω. We consider the guarantee given in (2.4) term by term.

For the rest of this paper, we take the kernels K₁ and K₂ to be Laplace-type kernels with quickly decaying singular values. For our purposes, we shall use the kernels k₁(τ₁, x) = 1 − e^−τ₁/x and k₂(τ₂, y) = e^−τ₂/y to represent the general data structure of most multiexponential NMR spectroscopy measurements. The same kernels shall be used in section 5 for simulations and experiments. Also, τ₁ is logarithmically sampled between 0.0005 and 4, and τ₂ is linearly sampled between 0.0002 and 0.4, as these are typical values in practice. Also for this section, F is taken to be a two-peak distribution, namely Model 3 from section 5.

When needed, we set s₁ = s₂ = 20. This choice is determined by the discrete Picard condition (DPC) [33]. For ill-conditioned kernel problems Kf = g, with {u_i} denoting left singular vectors of K and {σ_i} the corresponding singular values, the DPC guarantees that the best reconstruction of f is given by keeping all σ_i ≠ 0 such that $\frac{\mid u_1^{\ast }g\mid }{{\sigma }_i}$ on average decays to zero as σ_i decrease. For our kernels with tensor product structure in (1.1), Figure 1 shows the relevant singular values and vectors to keep. The s₁ = s₂ = 20 rectangle provides a close estimate for what fits inside this curve, implying that at a minimum we could set s₁ = s₂ = 20 to satisfy the DPC. The DPC provides a stronger condition than simply keeping the largest singular values or attempting to preserve some large percentage of the energy [32].

4. Algorithm

The algorithm for solving for F in (1.1) from partial data consists of three steps. An overview of the original algorithm in [56] is in section 1. Our modification and the specifics of each step are detailed below.

1. Construct M̃ from given measurements. Let y = ℛ_Ω(M̃₀) + e be the set of observed measurements, as in (2.7). Even though section 2 makes guarantees for solving (2.8), we can instead solve the relaxed Lagrangian form

   $$min\mu {\Vert X\Vert }_{\ast }+\frac{1}{2}{\Vert {\mathcal{R}}_{\mathrm{\Omega }}(X)-y\Vert }_2^2.$$

   (4.1)

   To solve (4.1), we use the singular value thresholding algorithm from [5, 46]. To do this, we need some notation. Let the matrix derivative of the L₂ norm term be written as

   $$\begin{array}{l}g(X)={\mathcal{R}}_{\mathrm{\Omega }}^{\ast }({\mathcal{R}}_{\mathrm{\Omega }}(X)-y)\\ =U_1^{\prime }({\mathcal{A}}_{\mathrm{\Omega }}^{\ast }({\mathcal{A}}_{\mathrm{\Omega }}(U_{1}X{U}_2^{\prime })-y))U_2.\end{array}$$

   We also need the singular value thresholding operator 𝒮_ν that reduces each singular value of some matrix X by ν. In other words, if the SVD of X = UΣV′, then

   $${\mathcal{S}}_{\nu }(X)=U\mathrm{\sum }^{\sim }{V}^{\prime },{\mathrm{\sum }^{\sim }}_{i,j}=\{\begin{array}{ll}max({\mathrm{\sum }}_{i,i}-\nu ,0),\hfill & \hfill i=j,\\ 0\hfill & \hfill \text{otherwise.}\end{array}$$

   Using this notation, the algorithm can then be written as a simple, two-step iterative process. Choose a τ > 0. Then, for any initial condition, solve the iterative process

   $$\{\begin{array}{l}{Y}^k=X^k-\tau g(X^k),\hfill \\ X^{k+1}={\mathcal{S}}_{\tau \mu }(Y^k).\hfill \end{array}$$

   (4.2)

   The choices of τ and μ are detailed in [46], along with adaptations of this method that speed up convergence. However, this method is guaranteed to converge to the correct solution.

   This means that, given partial observations y, the iteration scheme in (4.2) converges to a matrix M̃, which is a good approximation of M̃ +0. Once M̃ has been generated, we recover F by solving

   $$\underset{F\ge 0}{argmin}{\Vert \stackrel{\sim }{M}-(S_1{V}_1^{\prime })F{(S_2{V}_2^{\prime })}^{\prime }\Vert }_F^2+\alpha {\Vert F\Vert }_F^2.$$

   (4.3)

2. Optimization. For a given value of α, (4.3) has a unique solution due to the second term being quadratic. This constrained optimization problem is then mapped onto an unconstrained optimization problem for estimating a vector c.

   Let f be the vectorized version of F and m be a vectorized version of M̃. Then we define the vector c from f implicitly by

   $$f=max(0,K^{\prime }c),\text{where}K=(S_1{V}_1^{\prime })\otimes (S_2{V}_2^{\prime }).$$

   Here, ⊗ denotes the Kronecker product of two matrices. This definition of c comes from the constraint that F ≥ 0 in (4.3), which can now be reformed as the unconstrained minimization problem

   $$min\left(\frac{1}{2}{c}^{\prime }[G(c)+\alpha I]c-c^{\prime }m\right),$$

   (4.4)

   where

   $$G(c)=K\left[\begin{array}{cccc}H(K_{1,·}^{\prime }c)& 0& \cdots & 0\\ 0& H(K_{2,·}^{\prime }c)& \cdots & 0\\ ⋮& ⋮& & ⋮\\ 0& 0& \cdots & H(K_{N_x\times N_y,·}^{\prime }c)\end{array}\right]K^{\prime }$$

   and H(x) is the Heaviside function. Also, K_i,· denotes the ith row of K. The optimization problem (4.4) is solved using a simple gradient descent algorithm.

3. Choosing α. There are several methods for choosing the optimal value of α.

   • Butler–Reeds–Dawson (BRD) method. Once an iteration of step 2 has been completed, it is shown in [56] that a better value of α can be calculated by

     $${\alpha }_{\mathit{opt}}=\frac{\sqrt{s_1{s}_2}}{\Vert c\Vert }.$$

     If one iterates between step 2 and the BRD method, the value of α converges to an optimal value. This method is very fast; however, it can have convergence issues in the presence of large amounts of noise, as well as on real data [53].

   • S-curve. Let F_α be the value returned from step 2 for a fixed α. The choice of α should be large enough that F_α is not being overfitted and unstable to noise, yet small enough that F_α actually matches reality. This is done by examining the “fit error”

     $$\chi (\alpha )={\Vert M-K_1{F}_{\alpha }{K}_2^{\prime }\Vert }_F.$$

     This is effectively calculating the standard deviation of the resulting reconstruction. Plotting χ(α) for various values of α generates an S-curve, as shown in Figure 7. The interesting value of α occurs at the bottom “heel” of the curve (i.e., $\frac{dlog\chi (\alpha )}{dlog\alpha }\approx .1$). This is because, at α_heel, the fit error is no longer demonstrating overfitting as it is to the left of α_heel, yet is still matching the original data, unlike to the right of α_heel. This method is slower than the BRD method; however, it is usually more stable in the presence of noise.

For the rest of this paper, we use the S-curve method of choosing α.

5. Simulation results

In our simulations, we shall use the kernels k₁(τ₁, x) = 1 − e^−τ₁/x and k₂(τ₂, y) = e^−τ₂/y and sample τ₁ logarithmically and τ₂ linearly, as was done in section 3. Our simulations revolve around inverting subsampled simulated data to recover the density function F(x, y). We shall test three models of F(x, y). In Model 1, F(x, y) is a small variance Gaussian. In Model 2, F(x, y) is a positively correlated density function. In Model 3, F(x, y) is a two-peak density, one peak being a small circular Gaussian and the other being a ridge with positive correlation.

The data is generated for a model of F(x, y) by discretizing F and computing

where E is Gaussian noise. That data is then randomly subsampled by keeping only a λ fraction of the entries.

Each true model density F(x, y) is sampled logarithmically in x and y. τ₁ is logarithmically sampled N₁ = 30 times, and τ₂ is linearly sampled N₂ = 4000 times. Each model is examined for various SNRs and values of λ, and α is chosen using the S-curve approach for each trial.

Let us also define the signal-to-noise ratio (SNR) for our data to be

Note that [56] has extensively examined steps 2 and 3 of this algorithm, including the effects of α and the SNR on the reconstruction of F. Our examination focuses on the differences between the F generated from full knowledge of the data and the F generated from subsampled data. For this reason, F_full refers to the correlation spectrum generated from full knowledge of the data using the algorithm from [56]. F_λ refers to the correlation spectrum generated from only a λ fraction of the measurements using our algorithm.

6. Conclusion

In this paper, we introduce a matrix completion framework for solving 2D Fredholm integrals. This method allows us to invert the discretized transformation via Tikhonov regularization using far fewer measurements than previous algorithms. We proved that the nuclear norm minimization reconstruction of the measurements is stable and computationally efficient, and demonstrated that the resulting estimate of ℱ(x, y) is consistent with using the full set of measurements. This allows us in application to reduce the measurements conducted by a factor of 5 or more.

While the theoretical framework of this paper applies to 2D NMR spectroscopy, the approach is easily generalized to larger-dimensional measurements. This allows for accelerated acquisition of 3D correlation maps [3] that would otherwise take days to collect. This shall be a subject of forthcoming work.