Reading the 2007 Sparse MRI paper by Lustig, Donoho and Pauly

The problem the paper sets out to solve

Scan time is proportional to the number of k-space samples, and the number Nyquist sampling demands is what makes conventional scans slow. The paper starts from the observation that MR images are often sparse or compressible in a suitable transform: most coefficients are close to zero. If three requirements are met, fewer samples can be taken and the image computed back.

The three requirements: (a) the image has a sparse representation in a known transform domain; (b) the aliasing caused by k-space undersampling is incoherent, noise-like, in that domain; (c) a nonlinear reconstruction enforces both sparsity of the representation and consistency with the acquired data. Figure 1 of the paper shows how the image, k-space and transform domains are related.

Source: M. Lustig · MIT 6.3000 Signal Processing

Incoherence: the point spread function and Eq. [1]

Figure 2 shows the idea on a simple example: a sparse signal undersampled 8-fold. Equispaced undersampling with zero filling gives a set of shifted replicas, equally likely as the original and impossible to tell apart. Random undersampling gives artifacts that look like random noise but are really the energy of each nonzero coefficient leaking elsewhere. Strong components stand above the interference and can be found by thresholding; the interference they cause is computed and subtracted, and then the weaker components are found.

To quantify this, the paper defines the point spread function, PSF, and the transform point spread function, TPSF (Eq. [2]): how much of one transform coefficient leaks onto coefficient $j$ after undersampling. When $N$ of $D$ grid points are chosen at random, the real and imaginary parts of the off-centre values each behave much like zero-mean white Gaussian noise, with the standard deviation of Eq. [1], where $p=D/N$ is the reduction factor.

\sigma_{\mathrm{SPR}}=\sqrt{\frac{p-1}{D}}

Eq. [1] of the paper. The text says the derivation is in Appendix II; it is in Appendix B.

Source: M. Lustig · MIT 6.3000 Signal Processing

Three Cartesian cases and the sampling design

The paper compares three kinds of Cartesian sampling. Single-slice 2D imaging can only undersample the phase encodes, so the interference spreads along one direction and, by Eq. [1], is $D^{1/4}$ times larger than for pure 2D random sampling; only modest accelerations can be expected. Multislice 2D imaging can choose different phase-encode lines in each slice. The authors prefer 3D imaging, where random undersampling of the $k_y$ – $k_z$ plane reaches the high incoherence of the 2D case.

The sampling is designed by a Monte Carlo procedure: set the grid from the desired resolution and field of view, build a probability density that falls with a power of the distance from the origin, and draw samples from it. Because the draw is random it can produce a bad pattern, so the procedure is repeated many times and the set with the smallest peak TPSF sidelobe is kept. The reason for variable density is that image energy is concentrated near the k-space origin.

Source: M. Lustig · MIT 6.3000 Signal Processing

Reconstruction: constrained optimization and conjugate gradients

The reconstruction is Eq. [3]: among images that differ from the measured data by less than $\varepsilon$ , find the one with the smallest $\|\Psi m\|_1$ , with $\varepsilon$ usually set below the expected noise. Appendix A writes it in Lagrangian form, [A1], chooses $\lambda$ so the residual is about $\varepsilon$ , and solves it by nonlinear conjugate gradients with a backtracking line search; the examples need 80 to 200 iterations. Their Matlab implementation takes about 150 iterations, roughly 30 s, to reconstruct a $480\times92$ angiogram at 5-fold acceleration.

The sparsifying transform depends on the application: wavelets, finite differences (total variation) or both. The brain example uses Daubechies 4 wavelets with a total-variation penalty.

\min_m\ \|F_u m-y\|_2^2+\lambda\|\Psi m\|_1

Eq. [A1] of the paper.

Source: M. Lustig · MIT 6.3000 Signal Processing

Results and the limits the authors admit

Simulation: a $100\times100$ phantom with 18 randomly placed features of 6 sizes and 3 intensities, 575 nonzero pixels. At 8-fold and 12-fold acceleration (12-fold with variable density only) the compressed-sensing reconstruction is perfect; as the acceleration grows further, what is lost is the low-contrast features rather than the usual resolution.

Measurements: a multislice 2D fast-spin-echo brain scan, with 80 of 192 phase encodes chosen at random with quadratic variable density for each slice, a 2.4-fold acceleration; and a 3D abdominal angiogram undersampled 5-fold in retrospect, equivalent to cutting the scan from 22 s to 4.4 s.

The discussion is direct: as the acceleration grows, the most distinct artifacts of compressed sensing are not the usual loss of resolution or increase in aliasing interference but the loss of low-contrast features. The method is therefore particularly attractive where images have high-resolution, high-contrast features and imaging must be fast.

Source: M. Lustig · MIT 6.3000 Signal Processing

Worked example: recompute Eq. [1] and the factor of 16

The argument in Appendix B takes a few lines. Write the total energy of the point spread function as $N/D$ ; the main lobe holds $(N/D)^2$ , and the off-centre energy is the difference. Dividing by the number of off-centre pixels and by the main-lobe energy, and setting $p=D/N$ , gives Eq. [1]. Counting the off-centre pixels exactly as $D-1$ gives the right-hand expression below.

With $D=256$ and $N=32$ , Eq. [1] gives 0.165359 and the exact expression 0.165683. The course benchmark draws 400 random sample sets and measures 0.165683, the exact value. The excess kurtosis of the real part is −0.057, close to the Gaussian value of 0, in line with the paper’s statement that it behaves much like white Gaussian noise.

Now the factor of 16. A $256\times256$ image has $D=65\,536$ pixels. Pure 2D random sampling spreads the interference over all of them; single-slice 2D imaging can spread it only over the 256 pixels along one direction. The ratio is $\sqrt{65\,536/256}=16=D^{1/4}$ , as the paper states.

\sigma^2=\frac{N/D-(N/D)^2}{(D-1)(N/D)^2}=\frac{D-N}{N(D-1)}\;\approx\;\frac{p-1}{D}

For

D\gg1

this returns Eq. [1].

Source: M. Lustig · MIT 6.3000 Signal Processing

How this course’s experiment differs from the paper

The k-space experiment follows the paper’s approach: random phase-encode lines, the centre 8 lines always kept, the rest drawn with variable density, and an objective in Lagrangian form. But it simplifies in several places, so its numbers cannot be compared with the paper directly. The matrix is only $64\times64$ , with one slice, one coil and 2-fold acceleration. The sparsifying transform is a three-level Haar wavelet without total variation. The solver is an accelerated soft-threshold iteration with 50 steps, not conjugate gradients. The weight is chosen from held-out lines, not from the noise level.

With these simplifications the zero-filled error is 23.9% and the sparse error 13.3%, and the small-dot contrast recovers from 71% to 87.3%. The direction agrees with the paper: sparse reconstruction removes most of the incoherent artifacts, but the contrast is not fully restored. With regular skipping the sparse fit fails in the same way, the error staying at 49.8%, which corresponds to the unrecoverable equispaced case of the paper’s Fig. 2.

Source: M. Lustig · MIT 6.3000 Signal Processing

Try it in the Lab

Repeat the comparison of the paper’s Fig. 2 in the experiment

Choose “Skip every other line: the head folds”, then switch to the sparse wavelet fit. Predict that the fold does not go away. The error reads about 49.8%.
Choose “Random lines, gaps left empty” and note the 23.9% error; switch to “Random lines, sparse fit fills gaps” and the error falls to 13.3%. Both record the same number of lines.
In “Fit weight, centre rows and seed”, set the rows always kept at the centre to 0. Predict that the error grows without the centre rows. The readout should rise from 13.3% to 18.8%, and the zero-filled error from 23.9% to 30.9%.

Exercises with explained answers

Eighty of 192 phase encodes are recorded. What is the reduction factor, and, treating it as 1D random sampling, what interference standard deviation does Eq. [1] give?
Show the answer
$p=192/80=2.4$ and $\sigma=\sqrt{1.4/192}=0.085$ . The paper actually uses variable density, so the interference is not uniform and this number is only a scale.
In the paper’s simulation, what is lost first as the acceleration grows, and what does that mean for clinical use?
Show the answer
The low-contrast features go first; once the interference submerges them they cannot be recovered. Compressed sensing therefore suits high-contrast applications that need speed, such as angiography, and calls for care where diagnosis depends on subtle contrast.
Why do the authors prefer 3D imaging to single-slice 2D imaging?
Show the answer
3D imaging can undersample the $k_y$ – $k_z$ plane at random, spreading the interference in two dimensions for higher incoherence. Single-slice 2D imaging spreads it along one direction only, and the interference is $D^{1/4}$ times larger, 16 times for a $256\times256$ image.

Sources and what we recomputed

MIT 6.3000 Signal Processing, Magnetic Resonance Imaging lecture handout (6 May 2025), data and examples from E. Adalsteinsson
Read in full. Page 9: the proton gyromagnetic ratio rounded to 42.58 MHz/T. Pages 11–13: gradient fields, the phase integral and the definition of k as the time integral of the gradient. Pages 15–18: k-space data, the complex image and its magnitude. Page 19: 256 lines at 2 s each, about 8.5 minutes. Pages 20–21: keeping only even columns of k-space. Pages 24–36: two-coil unfolding, including the coil pair whose sensitivities cannot separate the folded pixels. The handout uses the sign convention of the Fourier kernel that this course follows.
M. Lustig, D. Donoho and J. M. Pauly, Sparse MRI: the application of compressed sensing for rapid MR imaging, Magnetic Resonance in Medicine 58, 1182–1195 (2007), doi:10.1002/mrm.21391
Read in full. Fig. 1 and the three requirements of compressed sensing; Fig. 2 (a 256-sample sparse signal, 8-fold undersampled); Eq. 1 (spread of the interference) and its derivation in Appendix B, which the text calls Appendix II; the remark that single-slice 2DFT interference is D^{1/4} larger; Eq. 2 (transform point spread); Eqs. 3–4 (constrained reconstruction); Appendix A, Eqs. A1–A3 (Lagrangian form and nonlinear conjugate gradients, 80 to 200 iterations); Results and Discussion (2.4-fold brain, 5-fold angiography from 22 s to 4.4 s, loss of low-contrast features). The in-vivo images are described, not reproduced.