I would very much like to implement the algorithm present in this paper, however, I think I am way over my head. My interest is to implement the algorithm contained in it to get 2D maps of two relaxation times, T1 and T2, used in the study of rocks with Nuclear Magnetic Resonance (NMR). The algorithm has three steps, and I want to fully understand the first one before continuing. I will describe the problem shortly, then show my work and say exactly where my problem is.

The idea of the paper is to solve Fredholm integrals of the first kind in order to get the joint probability density function $\mathcal{F}(x,y)$. The equation has this shape:

$M_r(\tau_1,\tau_2) = \int \int k_1(x,\tau_1)k_2(y,\tau_2)\mathcal{F}(x,y)dxdy+\epsilon_r(\tau_1,\tau_2)\quad r=1,...,R$

where $M$ is the measured data, $\tau_1$ and $\tau_2$ are both experimental time parameters, $x$ and $y$ are the relaxation times (T1 and T2) and $\epsilon$ is the error. $r$ is one "slice", but in my case I want only 2D data, so $r=R=1$. I will omit the subscript from now on.

For NMR, the kernels $k_1$ and $k_2$ are known, and have the form



The equation can be rewritten in a matrix form, resulting in


The first step of the algorithm involves using singular value decomposition (SVD) on $K_1$ and $K_2$ to reduce the size of the data. So, $K_1 = U_1\Sigma_1V_1'$. And here is where I am stuck. It seems obvious, as other papers treat it as such, but I don't know how to go from my measured data $M$ to the two kernels, to perform the SVD on them.

The paper provides some instructions to create a 2D dataset to visualize this step. I have reproduced it to the best of my abilities in Python. Here is the code.

import numpy as np
import matplotlib.pyplot as plt

t1s = np.array([500E-6, 0.07, 4])  # Three \tau_1 were provided
t2s = np.linspace(200E-6, 200E-6 * 4000, num=4000)  # Typical 4000 echoes (\tau_2) were provided

def kernel1(t1, x):
    return 1 - 2 * np.exp(-t1/x)

def kernel2(t2, y):
    return np.exp(-t2/y)

# I guessed two x and y values, called relax_1 (for \tau_1) and relax_2 (for \tau_2)
# I guessed them to get curves similar to the ones in the paper
relax_1 = 3e-1
relax_2 = 2e-1

# Calculate M
signal = np.zeros((3, 4000))
for ii in range(len(t1s)):
    for jj in range(len(t2s)):
        signal[ii, jj] = (
            kernel1(t1s[ii], relax_1) *
            kernel2(t2s[jj], relax_2)
signal += np.random.randn(*signal.shape) / 20  # To add some noise

# Plot the results
plt.plot(t2s, signal[0, :], label=f't1={t1s[0]}')
plt.plot(t2s, signal[1, :], label=f't1={t1s[1]}')
plt.plot(t2s, signal[2, :], label=f't1={t1s[2]}')
plt.xlabel(r'$\tau_2$, [s]')

Comparison between my dataset and theirs As you can see, I am very close to their dataset. However, I am stuck in the next part, using SVD on $K_1$, because it's unclear to me what $K_1$ is. I have only my data, how would I "decompose" it into two kernels to perform SVD?

For example, if I try

from numpy.linalg import svd
U1, S1, V1 = svd(signal, full_matrices=True)
plt.plot(range(len(S1)), S1)

$\Sigma_1$ will have three singular values. How could $\Sigma_2$ have, according to their paper, 10 singular values, if the input matrix has the shape (3,4000)? If someone could explain this part to me, and, if it's not too much to ask, provide a snippet of code of how to do it, I would be very grateful.

Thanks a lot for reading up to this point. I expect to come back here shortly to ask about the next steps in the algorithm.


1 Answer 1


So, to answer my own question, for posterity, in case anyone faces the same problem as me.

First, the data generation is wrong. I created two kernel matrices with known relaxation times ($x$ and $y$), along $\tau_1$ and $\tau_2$. Then, to get the signal, I just multiplied both, without considering the function $F$. This was the main crux of what I didn't understand.

The point is that the function $F$ will actually give me the $x$ and $y$ distributions, and the kernels only serve to "map out" the possibilities for the algorithm to consider. So the kernels are, in essence, unchanged between runs of the algorithm, or between different measured signals $M$ (unless you decide to change the possible $x$ and $y$ values, which is not that often honestly, for NMR at least).

To generate a kernel, then, is easy; just calculate every single exponential possible, but not attribute any weights to them! These weights will come later from $F$. This can be done in Python in the following manner.

N1, N2 = 30, 4000
Nx, Ny = 100, 100
tau1s = np.logspace(np.log10(5E-4), np.log10(4), num=30)
tau2s = np.arange(2E-4, (N2+1) * 2E-4, step=2E-4)
xs = np.logspace(np.log10(0.05), np.log10(0.4), num=Nx)
ys = np.logspace(np.log10(0.05), np.log10(0.4), num=Ny)
kern1 = np.zeros((N1, Nx))
kern2 = np.zeros((N2, Ny))
for ii in range(N1):
    kern1[ii, :] = kernel1(tau1s[ii], xs) 
for jj in range(N2):
    kern2[jj, :] = kernel2(tau2s[jj], ys)

Therefore, you get $K_1$ and $K_2$, with each row being one possible decay curve.

Then, to get the matrix $M$, you can use one of the examples from the article as the matrix $F$, a small Gaussian centered at $\mu_x = 0.1$ and $\mu_y = 0.2$.

X,Y = np.meshgrid(xs,ys)
mux = 0.2
muy = 0.1
# https://en.wikipedia.org/wiki/Gaussian_function
Ft = np.exp( -( (X - mux)**2 + (Y - muy)**2) / (2.0 * sigma**2 ))

At last, to get the perfect noiseless signal, just multiply these matrices, S = kern1 @ Ft @ kern2.T.

Then, you can easily perform the SVD of $K_1$ and $K_2$ and follow the article in a reasonably straightforward way.


You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .