PaperlayerView exact revision on arXiv

Source-linked AI summary

Bayesian Nonparametric Dictionary Learning for Compressed Sensing MRI

Yue Huang, John Paisley, Qin Lin, Xinghao Ding, Xueyang Fu, Xiao-ping Zhang

arXiv:1302.2712v3cs.CVphysics.med-phstat.AP

TL;DR

The paper addresses CS-MRI reconstruction from highly undersampled k-space, where fixed dictionary settings may not suit the image. It uses BPFA to learn image-specific sparse patch dictionaries, optionally combined with TV. The authors report favorable comparisons with other algorithms across MRI datasets and sampling patterns.

  • Problem

    Dictionary-learning methods often require preset dictionary size and patch sparsity, although mismatched settings can degrade performance in CS-MRI reconstruction.

    Source: §I · p. 2

  • Method

    BPFA learns a sparse patch dictionary during reconstruction, inferring its effective size and patch-specific sparsity, and can be combined with TV regularization.

    Source: 1) Generative model · p. 3, §III · p. 4

  • Results

    The proposed framework compares favorably with other algorithms on real and complex-valued MRI across sampling trajectories and rates.

    Source: §V · p. 12

  • Takeaways & Limitations

    In noisy settings, the denoising property of dictionary learning allows the measured k-space values to be fixed while BPFA and TV fill in the reconstruction.

    Source: B. Discussion on λ · p. 6

Abstract

from arXiv · show

We develop a Bayesian nonparametric model for reconstructing magnetic resonance images (MRI) from highly undersampled k-space data. We perform dictionary learning as part of the image reconstruction process. To this end, we use the beta process as a nonparametric dictionary learning prior for representing an image patch as a sparse combination of dictionary elements. The size of the dictionary and the patch-specific sparsity pattern are inferred from the data, in addition to other dictionary learning variables. Dictionary learning is performed directly on the compressed image, and so is tailored to the MRI being considered. In addition, we investigate a total variation penalty term in combination with the dictionary learning model, and show how the denoising property of dictionary learning removes dependence on regularization parameters in the noisy setting. We derive a stochastic optimization algorithm based on Markov Chain Monte Carlo (MCMC) for the Bayesian model, and use the alternating direction method of multipliers (ADMM) for efficiently performing total variation minimization. We present empirical results on several MRI, which show that the proposed regularization framework can improve reconstruction accuracy over other methods.

I. INTRODUCTION

Compressed sensing can accelerate MRI by reconstructing images from undersampled k-space, but image-level fixed bases may not fit the varied features in MR images. The paper proposes in situ Bayesian nonparametric dictionary learning for CS-MRI, inferring a compact dictionary and patch sparsity from the data.

I. INTRODUCTION

  • Undersampled k-space enables faster MRI acquisition, and compressed sensing can reconstruct high-quality images when signals are sparse in a transform domain.The inverse problem is ill-posed, but CS permits reconstruction from fewer measurements than Nyquist sampling requires.Source: §I · p. 1
  • Image-level CS-MRI methods often use fixed bases that capture only selected features, although MR images contain varied structures such as edges and textures.Sparse MRI, for example, constrains both wavelet coefficients and total variation.Source: §I · p. 1
  • Dictionary learning adapts a sparse patch-level basis to the image, and CS theory links signal sparsity in the chosen transform domain to the required measurement count.Adaptive dictionary learning has also been effective in applications including denoising, super-resolution, interpolation, and inpainting.Source: §I · pp. 1–2
  • BPFA replaces preset dictionary size and patch sparsity with a Bayesian nonparametric model that infers these from the data and learns a compact dictionary.The paper uses MCMC for stochastic optimization of dictionary variables and ADMM for efficient TV minimization.Source: §I · p. 2

A. Two approaches to CS-MRI inversion

CS-MRI regularization can impose sparsity globally in a transform domain or locally through sparse dictionary representations of patches. The distinction is whether the representation is predefined for the whole image or learned from patch data.

A. Two approaches to CS-MRI inversion

1) Image-level sparse regularization:

  • Image-level regularization enforces sparsity in a transform domain over the entire image, as in Sparse MRI's wavelet-plus-TV penalty.The wavelet term encourages sparse coefficients while TV encourages neighboring pixels to be similar but permits sharp edge changes; μ controls the trade-off.Source: 1) Image-level sparse regularization · p. 2

2) Patch-level sparse regularization:

  • Patch-level regularization instead requires image patches to have sparse representations in a dictionary that can be learned from data.Unlike a predefined transform, the dictionary domain can be unknown; K-SVD is one example that learns the dictionary offline and uses OMP.Source: 2) Patch-level sparse regularization · p. 2

B. Dictionary learning with beta process factor analysis

BPFA addresses dictionary learning's need to preset dictionary size and patch sparsity by inferring them from data. Its generative model represents each image patch as a sparse weighted combination of dictionary elements.

B. Dictionary learning with beta process factor analysis

  • Typical dictionary-learning methods preset dictionary size and patch sparsity or an error threshold, and mismatched settings can significantly degrade performance.BPFA is introduced as a Bayesian nonparametric alternative that can infer these values.Source: B. Dictionary learning with beta… · pp. 2–3

1) Generative model:

  • BPFA assigns each dictionary element a probability and models each patch with a sparse weighted combination of elements, using binary indicators and Gaussian weights.Most element probabilities are expected to be small, so unused elements can be identified from the data; K is set above the expected dictionary size.Source: 1) Generative model · p. 3
  • The model learns which dictionary elements are used across patches, with patch-specific sparsity determined by the element-wise product of indicators and weights.In the infinite-dictionary limit, the number of elements used per patch follows a Poisson distribution with parameter γ.Source: 1) Generative model · p. 3

2) Relationship to K-SVD:

BPFA differs from K-SVD by pruning unused dictionary elements and updating sparsity through a Bernoulli-process posterior. In the reconstruction framework, BPFA supplies local regularization and TV supplies a global smoothness penalty, with MCMC and ADMM used for optimization.

2) Relationship to K-SVD:

  • Unlike K-SVD's OMP-based sparsity selection and fixed supplied dictionary size, BPFA prunes unused elements and updates sparsity using a Bernoulli-process posterior.BPFA also samples weights and dictionary elements from Gaussian posteriors and learns auxiliary variance variables.Source: 2) Relationship to K-SVD · p. 3

3) Example denoising problem:

  • BPFA denoises by learning patch representations and averaging their dictionary reconstructions to form the image.The example uses 6 × 6 patches and truncates the dictionary at K = 108; fewer than 100 elements are used by the data.Source: 3) Example denoising problem · pp. 3–4
  • BPFA adaptively infers the mismatched noise parameter in the denoising example, whereas K-SVD's performance suffers when its error parameter is set incorrectly.The example compares K-SVD settings that match or mismatch the ground-truth noise parameter.Source: 3) Example denoising problem · p. 4

III. CS-MRI WITH BPFA AND TV PENALTY

  • The reconstruction framework uses BPFA as local regularization and total variation as global regularization, with TV encouraging homogeneity while allowing sharp edges.The trade-off parameters include λg, γε, and λ; γε is sampled at each iteration, so the trade-off changes adaptively.Source: §III · p. 4
  • ADMM efficiently optimizes the isotropic TV penalty by performing dual ascent on an augmented Lagrangian for the TV coefficients.The coefficients are formed from image differences between neighboring pixels.Source: §III · p. 4

A. Algorithm

The nonconvex objective is optimized by cycling among TV, BPFA, and image-reconstruction subproblems. Gibbs sampling handles BPFA variables that cannot be optimized analytically.

A. Algorithm

  • The objective separates into TV, BPFA, and image-reconstruction subproblems, which are cycled through using ADMM.The algorithm repeats the cycle until convergence.Source: A. Algorithm · pp. 4–5
  • P1 and P3 are globally optimal conditional on the latest values of the other parameters, while P2 has no analytic solution for all BPFA variables.For fixed D, α_1:N, and x, the solution is also globally optimal.Source: A. Algorithm · p. 4
  • The algorithm uses Gibbs sampling to stochastically optimize the BPFA subproblem, with one sample per iteration.A single sample was found sufficient in practice and faster than taking more samples.Source: A. Algorithm · p. 5

1) Algorithm for P1 (total variation):

The algorithm separates total-variation coefficients from image gradients and uses BPFA Gibbs updates for dictionary learning and sparse coding.

1) Algorithm for P1 (total variation):

  • For P1, β_i represents the vertical and horizontal TV coefficients for pixel i, and the algorithm updates the corresponding Lagrange multiplier.The coefficients are defined by splitting β_i from ψ_i x through an augmented Lagrangian.Source: 1) Algorithm for P1 (total variation) · p. 5

2) Algorithm for P2 (BPFA):

  • For P2, Gibbs sampling updates the complex dictionary, patch-specific binary indicators and weights, noise and weight precisions, and element-use probabilities.The weights and indicators form α_i by element-wise multiplication, and one sample was sufficient in practice for a faster algorithm.Source: 2) Algorithm for P2 (BPFA) · p. 5
  • The dictionary update approximates the image patches X with Dα plus complex Gaussian noise, while the element-use probabilities depend on counts of patches using each element.The posterior update adds correlated complex-plane Gaussian noise to a regularized least-squares dictionary estimate.Source: 2) Algorithm for P2 (BPFA) · p. 5
  • Sparse coding samples z before s given z, a direction found more efficient for finding an objective mode; the inclusion probability for z_ik reflects how well dictionary element d_k correlates with the residual.When z_ik = 0, s_ik is sampled from the prior and does not enter the model.Source: 2) Algorithm for P2 (BPFA) · p. 5

3) Algorithm for P3 (MRI reconstruction):

The reconstruction subproblem has a closed-form solution that is computed efficiently in the Fourier domain. The resulting Fourier coefficients are updated element-wise before an inverse transform produces the image.

3) Algorithm for P3 (MRI reconstruction):

  • The image update is a least-squares problem with a closed-form solution, but direct inversion of its N × N matrix is computationally prohibitive.The formulation includes the BPFA-denoised reconstruction x_BPFA alongside the TV and data-fidelity terms.Source: 3) Algorithm for P3 (MRI reconstruction) · p. 6
  • Working in the Fourier domain diagonalizes the left-hand matrix, enabling element-wise updates of the Fourier coefficients followed by an inverse Fourier transform.The diagonal structure follows from the finite-difference operator and the measured k-space locations.Source: 3) Algorithm for P3 (MRI reconstruction) · p. 6
  • The inverse Fourier transform of the updated coefficients gives the reconstructed MR image.The update uses measured k-space values; unmeasured locations do not contribute to the corresponding measurement term.Source: 3) Algorithm for P3 (MRI reconstruction) · p. 6

B. Discussion on λ

In noisy CS-MRI, the fidelity parameter λ can remain large because BPFA supplies a denoised image, avoiding the usual need to tune λ for denoising.

B. Discussion on λ

  • In noisy reconstruction, λ controls denoising for many CS-MRI algorithms, but dictionary learning allows it to remain large without parameter selection.A large λ effectively fixes measured k-space values while BPFA and TV fill missing values.Source: B. Discussion on λ · p. 6
  • With large λ, the reconstruction x may be noisy but has subsampling artifacts removed, while x_BPFA is the denoised version of x.The paper attributes this denoising to BPFA.Source: B. Discussion on λ · p. 6
  • The algorithm’s quality therefore depends largely on BPFA’s image-denoising quality.The paper points to synthetic and clinical examples for this behavior.Source: B. Discussion on λ · p. 6

IV. EXPERIMENTAL RESULTS

The experiments compare BPFA with and without TV against established CS-MRI methods across synthetic and real-valued MRI, sampling masks, and noise settings. They report competitive reconstruction performance and examine dictionary learning’s denoising behavior.

IV. EXPERIMENTAL RESULTS

  • The evaluation covers real-valued and complex-valued MRI and a synthetic phantom, using PSNR and SSIM and comparisons with Sparse MRI, DLMRI, PBDW, and BPFA without TV.Sampling includes 2D random, Cartesian with random phase encodes, and pseudo-radial masks.Source: §IV · p. 6

A. Set-up

  • Experiments use 6 × 6 patches, initialize with k-space zero-filling, start with K = 108 dictionary elements, and run 1000 iterations.The final dictionary can be smaller under the sparse BPFA prior; λ = 10^100 effectively enforces equality with measured k-space values.Source: A. Set-up · p. 7

B. Experiments on a GE phantom

  • With 30% Cartesian sampling and σ = 0.1 noise, BPFA’s xBPFA output denoises the reconstruction; roughly 80 dictionary elements are used, averaging 2.28 active elements per patch.The best TV reconstruction required sweeping λg with λ set to 1, whereas BPFA’s noisy-image reconstruction used λ = 10^100.Source: B. Experiments on a GE phantom · p. 7, Figure 3 · p. 7

C. Experiments on real-valued (synthetic) MRI

  • The PSNR results show competitive performance across sampling masks and rates, a slight improvement with TV, and more noise-like residuals for the proposed method.The evaluated masks rank 2D random highest, followed by pseudo-radial and Cartesian sampling; the real-valued MRI experiments are described as idealized rather than clinically representative.Source: C. Experiments on real-valued… · pp. 7–8, Table II · p. 8, Table III · p. 8

D. Experiments on complex-valued MRI

The complex-valued MRI experiments report competitive reconstruction performance and assess noisy data with different noise levels. The denoised BPFA reconstruction remains stable as the fidelity parameter varies, unlike the reconstruction that enforces fit to noisy measurements.

D. Experiments on complex-valued MRI

E. Experiments in the noisy setting

  • For noisy brain MRI with σ = 0.01, 0.02, and 0.03, the experiments compare BPFA with and without TV against the noise-free ground truth.The paper reports PSNR by noise level for 35% Cartesian sampling.Source: E. Experiments in the noisy setting · p. 10, Table V · p. 10
  • As λ increases, the reconstruction x eventually degrades from enforcing a closer fit to noisy measurements, while xBPFA levels off without degrading.At small λ, both reconstructions degrade because dictionary learning smooths x too much; the authors report robustness to λ and no stability issues with extremely large values.Source: Figure 8 · p. 9, E. Experiments in the noisy setting · p. 10

F. Dictionary learning and further discussion

The conclusion reports favorable comparisons across MRI settings and a noisy-measurement reconstruction that avoids dependence on λ. The discussion shows that BPFA adapts dictionary use to sampling rate and patch content, while noting implementation-dependent runtime comparisons.

F. Dictionary learning and further discussion

  • At 10%, 20%, and 30% sampling, roughly 70, 80, and 95 of 108 dictionary elements were significantly used, respectively, with more elements used per patch at higher sampling rates.The learned dictionaries were similar in shape, while the distribution of per-patch usage showed separate background and foreground modes; the foreground mode increased with sampling rate.Source: F. Dictionary learning and further… · p. 11, Figure 10 · p. 11
  • Results were not very sensitive to patch size, and BPFA generally reached high quality in half the iterations; its speedup over DLMRI was attributed to omitting OMP.The runtime comparison used an Intel Xeon CPU E5-1620 at 3.60GHz with 16.0G RAM; the authors note that they may have run more iterations than necessary.Source: F. Dictionary learning and further… · p. 11
  • BPFA automatically adapts dictionary size and per-patch sparsity, while MCMC also adjusts noise variance and score-weight variances.The authors connect this denoising-related modeling to CS-MRI with BPFA.Source: F. Dictionary learning and further… · p. 11
  • MCMC sampling is presented as a potential advantage in highly non-convex Bayesian models, where such methods are reported to rarely perform worse than deterministic methods.The discussion notes that sampling techniques are readily derived for BPFA because it is Bayesian.Source: F. Dictionary learning and further… · pp. 11–12
  • The reported speed advantage over K-SVD in Matlab may be smaller in other coding languages.The authors attribute BPFA's Matlab speed over K-SVD to requiring fewer loops.Source: F. Dictionary learning and further… · p. 11

V. CONCLUSION

  • The framework compares favorably with other algorithms across sampling trajectories and rates, and handles noisy MRI without dependence on λ by using a near-equality measurement constraint and denoised dictionary output.The conclusion reports experiments on real and complex-valued MRI.Source: §V · p. 12

VI. APPENDIX

The appendix explains how ADMM separates a smooth squared-error term from a nonsmooth penalty by splitting variables and relaxing their equality constraint.

VI. APPENDIX

  • ADMM separates the smooth squared-error term from a nonsmooth convex penalty, such as an ℓ1 penalty, by introducing a second vector v.The equality v = x is then relaxed using an augmented Lagrangian term.Source: §VI · p. 12
  • The augmented formulation seeks a minimax saddle point by minimizing over x and v while performing dual ascent for η.Source: §VI · p. 12
  • Cycling through x, v, and u updates simplifies optimization: the x objective is quadratic with a simple analytic solution, while the v update is a proximity operator.For a TV penalty, the v solution is analytical.Source: §VI · p. 12
Loading 1302.2712v3…