Data-driven approaches recently achieved remarkable success in mri reconstruction, but integration into clinical routine remains challenging due to a lack of generalizability and interpretability.
In this paper, we address these challenges in a unified framework based on generative image priors.
We propose a novel deep neural network based regularizer which is trained in a generative an unsupervised setting on reference magnitude images only.
After training, the regularizer encodes higher-level domain statistics which we demonstrate by synthesizing images without data.
Embedding the trained model in a classical variational approach yields high-quality reconstructions irrespective of the sub-sampling pattern.
In addition, the model shows stable behavior when confronted with out-of-distribution data in the form of contrast variation.even if the test data deviate significantly from the training data.
Furthermore, a probabilistic interpretation provides a distribution of reconstructions and hence allows uncertainty quantification.
To reconstruct parallel mri, we propose a fast algorithm to jointly estimate the image and the sensitivity maps.
The results demonstrate competitive performance, on par with state-of-the-art end-to-end deep learning methods, while preserving the flexibility with respect to sub-sampling patterns and allowing for uncertainty quantification.

§ INTRODUCTION

A vast literature is dedicated to reducing examination time in mri while retaining the diagnostic value of the resulting images.
On the hardware side, pi <cit.> exploits spatially varying sensitivity maps of coil arrays.
Such pi hardware has become standard in clinical systems, but noise amplification limits the potential speed up with classical reconstruction techniques <cit.>.
On the algorithmic side, cs theory <cit.> and variational approaches enable greater acceleration under the assumption that the reconstruction has a sparse representation in some basis.
Often, this basis is hand-crafted and only reflects rudimentary aspects of the prior information in the underlying distribution.
As a prominent example, the tv assumes sparsity in the spatial image gradients (which translates to piecewise constant image intensities), and has been successfully applied to parallel mri <cit.>.

Hand-crafted prior information, in general, fails to capture the complexity of the underlying distribution <cit.>.
On the other hand, data-driven approaches have shown promising results mri reconstruction, often leading to superior results compared to classical cs techniques <cit.>.
Here, information is encoded in a deep neural network, which is trained on reference data in an off-line step.
Data-driven approaches have successfully been applied as pre-processing steps in \( k \)-space <cit.> as well as post-processing steps in image-space <cit.>.
vn <cit.> imitate the structure of an iterative reconstruction scheme by unrolling an optimization algorithm.
In a dynamical setting, this approach has been proven useful for quantifying flow in four-dimensional mri <cit.>.
At the pinnacle, purely data-driven methods disregard any physical measurement model and instead learn direct maps from \( k \)-space to image-space (AUTOMAP, <cit.>).

(kspace1) at (0 * + , -1 * + ) ;
at (1 * + , -2 * + ) ;
[inner sep=0, outer sep=0] (kspace3) at (2 * + , -3 * + ) ;
at ($(kspace1.north) + (0.2, 0.1)$) \( z \in \mathbb{C}^m \);
[inner sep=0, outer sep=0] (kspace6)at (5 * + , -6 * + ) ;
ıin 1, 2, 3

at ($(kspace3.south west)!\fr!(kspace6.south west)$) \( \cdot \);
at ($(kspace3.north east)!\fr!(kspace6.north east)$) \( \cdot \);

[count=ı] in 00, 20, 40, 90

[inner sep=0, outer sep=0] (reconstructionı) at (1cm + ı* + ı* , 0.3)

at ($(reconstruction\i.north) + (0, 0.2cm)$)
ı=1\( u^{(0)} \in \R^n_{\geq0} \)
ı<4\( u^{(\which)} \)\( u^{(K)} \eqqcolon \optimal{u} \)

at (1cm + ı* + ı* , -2.5)
ı=1\( \Sigma^{(0)} \)
ı<4\( \Sigma^{(\which)} \)\( \Sigma^{(K)} \eqqcolon \Sigma^* \)

at (0.9cm + ı* + ı* , -1.95) \( \ddots \);
[count=] in 00, 01, 02, 14

[inner sep=0, outer sep=0] (sensitivity mapsı) at (0.7cm + ı* + ı* + * , -* - 1.0cm) ;

[rectangle, draw, minimum width=0.6cm, minimum height=0.6cm] (initial guess box) at ($(reconstruction1.west) + (-1, 0)$) (<ref>);
[latex-] (initial guess box) – ++(-1, 0);
[-latex] (initial guess box) – (reconstruction1);
[-latex] (initial guess box) |- (sensitivity maps14.west);
ıin 1, 2, 3

ȷı+ 1
[inner sep=0,outer sep=0] (midı) at ($(reconstruction\i)!.5!(reconstruction\j)$) \( \square \);
(reconstructionı) – (midı) [-latex] – (reconstructionȷ);
[latex-,dashed] (midı) – ++(0, 0.5) node [above] \( z \);

[on background layer]
[thick, draw=green!20!black, opacity=0.5, fill=green!50!black] ($(reconstruction1.north west)+(-0.3, 0.5)$) rectangle ($(reconstruction4.south east) + (+0.3, -1.5)$);

[shift=(16, 1.6)]
(upper) [draw, rounded corners] at (0, 0) \( \bigl( \, \cdot \, - L^{-1}_u \nabla H(\,\cdot\,, \Sigma^{(k)}) \bigr) \);
[->] (-1.1, 0.8) – (upper.north -| -1.1, 2);
at (-1.3, 1.) \( u^{(k)} \);
(proxU) [draw, rounded corners] at (0, -1) \( u^{(k+1)} = (\, \cdot \,)_+ \);
(nablaS) [draw, rounded corners] at (0, -2.3) \( \bigl( \, \cdot \, - L^{-1}_\Sigma \nabla H(u^{(k+1)}, \, \cdot \,)\bigr) \);
(1.1, 0.8) – ($(upper.north -| 1.1, 0.8) + (0.,0.2)$);
[loosely dashed] ($(upper.north -| 1.1, 0.8) + (0.,0.15)$) – ++(0, -2.1);
[->] ($(upper.north -| 1.1, 0.8) + (0.,-1.95)$) – (nablaS.north -| 1.1,0);
at (1.1, 1) \( \Sigma^{(k)} \);
(lower) [draw, rounded corners] at (0, -3.3) \( \mathcal{S}^{-1} \bigl( \operatorname{diag}(\xi_i + \mu)^{-1} \mathcal S(\mu\,\cdot\,) \bigr) \);
[on background layer]
[thick, draw=black, fill opacity=0.5, fill=green!50!black] ($(upper.north west -| lower.south west)+(-0.3, 0.3)$) rectangle ($(lower.east |- lower.south)+(0.3, -0.3)$);

[on background layer]
[thick, draw=blue!20!black, opacity=0.5, fill=blue!30!black] ($(upper.north west -| lower.south west)+(-0.15, 0.15)$) rectangle ($(lower.east |- proxU.south)+(0.15, -0.15)$);

[on background layer]
[thick, draw=red!20!black, opacity=0.5, fill=red!30!black] ($(nablaS.north west -| lower.west)+(-0.15, 0.15)$) rectangle ($(lower.east |- lower.south)+(0.15, -0.15)$);

[<-, dashed] (upper.east) – (2.5, 0) node [right] \( z \);
[<-, dashed] (nablaS.east) – (2.5, 0 |- nablaS.east) node [right] \( z \);
[->] (upper) – (proxU);
[->] (nablaS) – (lower);
(proxU) -| (-1.3, -1.7);
[loosely dashed] (-1.3, -1.8) – ++(0, -2);
[->] (-1.3, -3.9) – (-1.3, -4.2) node [below] \( u^{(k + 1)} \);
[->] (lower.south -| 1.1, 0) – (1.1, -4.2) node [below] \( \Sigma^{(k + 1)} \);

[dashed, opacity=0.5] (mid3.north west) to [in=-110,out=20] ($(upper.north west -| lower.south west)+(-0.3, 0.3)$);
[dashed, opacity=0.5] (mid3.south west) to [in=110,out=-20] ($(lower.south west |- lower.south)+(-0.3, -0.3)$);

Sketch of the reconstruction algorithm:
To jointly reconstruct the spin density \( u \) and the sensitivity maps \( \Sigma \), we impose data-fidelity, image-regularity, and coil-regularity in the iterations of ipalm <cit.>.
The function \( H \) incorporates our learned regularizer acting on \( u \).
Details are discussed in <ref>.

Data-driven methods typically require large training datasets.
In particular, state-of-the-art networks mimicking iterative schemes like the end-to-end vn <cit.> as well as AUTOMAP <cit.> (and related methods such as <cit.>) require image-data pairs for training.
Such image-data pairs are scarce in medical applications, whereas reconstructed images are much more abundantly available in the form of DICOM data <cit.>.
Data-driven methods that only assume access to reference images (but still allow for some kind of data-fidelityconsistency) include gan <cit.> and score-based diffusion models <cit.>.
Both methods share the problem that it is not obvious how the encoded prior information is best used:
gan suffer from the range-dilemma <cit.>

and authors have proposed to optimize the parameters of the gan at inference time <cit.>, effectively turning it into a deep image prior <cit.>.
Diffusion models have shown remarkable results in mri reconstruction <cit.>, but it is still an open question how to optimally incorporate data-fidelity in the reverse diffusion process (see <cit.> and the discussion in <ref> for an overview of proposed methods).

, and for score-based diffusion models it is often not clear which distribution the reconstructions are drawn from.

In addition, (with the exception of score-based diffusion models) all of the mentioned models effectively act as a point estimator, mapping a \( k \)-space datum to an image.
However, having access to a distribution of reconstructions and, consequently, being able to quantify uncertainty in the reconstruction is of utmost importance in medical applications.
Further, these point estimators typically assume a particular acquisition modality and do not generalize with respect to acquisition masks <cit.>.
In fact, it has been shown that reconstruction quality of some data-driven approaches deteriorates when more data becomes available <cit.>.
This severely hampers adoption of these methods in clinical practice.
In any case, a method that marries the great representation power of data-driven deep neural networks with the interpretability of variational approaches is still sought after.

In this work we pursue a principled approach to parallel mri that combines the merits of modern data-driven approaches with the benefits of classical variational approaches.
In particular, we learn a highly expressive image-space prior in the well-known maximum-likelihood framework.
In contrast to most other data-driven approaches, training our model does not require image-data pairs, but only assumes access to a database of reference reconstructions.
To solve different pi reconstruction problems, we propose a novel algorithm to jointly estimate the image as well as the sensitivity maps in a sense-type <cit.> forward model.

A sketch of our proposed approach is shown in <ref>.

We demonstrate that combining this prior with suitable data-likelihood terms yields competitive performance for a multitude of pi problems.
Specifically, the variational formulation makes our method agnostic to sub-samplingacquisition patterns.
Additionally, the accompanying probabilistic interpretation naturally enables experts to explore the full posterior distribution of any reconstruction problem.
In particular, we analyze the posterior expectation as well as the pixel-wise marginal posterior variance.

§.§ On the Hunt for the Optimal Regularizer

The classical variational approach for solving inverse problems amounts to minimizing an energy functional

\begin{equation}
\min_x E(x, z) \coloneqq D(x, z) + R(x).%
\label{eq:variational}
\end{equation}

Here, \( x \) is the reconstruction given the datum \( z \).
The data-fidelity term \( D \) models the physical acquisition process, including stochasticity.
For example, if \( A \) models the physical acquisition process, \( D(x, z) = \sfrac{1}{2}\norm{Ax - z}_2^2 \) is optimal under the assumption of Gaussian measurement noise.
On the other hand, the regularization term \( R \) encodes prior knowledge about the solution.
Classical regularizers \( R \) include magnitude penalization (e.g. the Tikhonov functional \( x \mapsto \norm{x}_2^2 \)) or encode smoothness assumptions (e.g. the tv \( x \mapsto \sum_{i=1}^n \sqrt{{(\mathrm{D} x)}_{i}^2 + {(\mathrm{D} x)}_{i + n}^2} \)).

Often, <ref> is used without fully acknowledging its probabilistic interpretation:
It computes the map estimate of the posterior \( p(x\mid z) \propto p(z \mid x) p_R(x) \) composed of a likelihood \( -\log p(z \mid x) \propto D(x, z) \) and a prior \( -\log p_R(x) \propto R(x) \).
Hand-crafted regularizers, in general, can not faithfully represent the intricate statistics of the underlying distribution.
Since these regularizers are overly simplistic, methods like early stopping the optimization of <ref> have emerged <cit.>.
Clearly, the need for such heuristics arises only because hand-crafted regularizers are a bad model for natural or medical images.

In this paper we pose the following question:
Can we learn a regularizer such that \( p_R \) is indistinguishable from the underlying reference distribution?
We answer this by encoding the regularizer as a neural network endowed with learnable parameters, which we train on reference data using maximum likelihood (see <ref>).

This method of learning an unnormalized probability density is known as learning an ebm <cit.>, and has recently been used in the context of mri reconstruction by <cit.>.

After training, the network can be used in arbitrary reconstruction tasks by embedding it in <ref>, thereby combining the versatility of variational approaches with the great representation power of deep neural networks.

On one hand, concerning versatility, the strict separation of data likelihood and prior makes this approach agnostic to the
number of receiver coils and sub-sampling pattern.

acquisition protocol.
We can reconstruct images from single- and multi-coil data with arbitrary acquisition masks.

Further, the probabilistic interpretation enables access to a distribution of reconstructions \( p(x \mid z) \), and we can compute different Bayesian estimators.
In addition to the map estimate (\( \delta_x \) is the Dirac measure centered at \( x \))

\begin{equation}
\max_x \int p(x^\prime \mid z)\ \mathrm{d}\delta_x(x^\prime),
\end{equation}

we exploit the distribution to compute the mmse estimate

\begin{equation}
\min_x \int  p(x^\prime \mid z)\norm{x - x^\prime}^2 \ \mathrm{d}x^\prime \eqqcolon \mathbb{E}[x \mid z],%
\label{eq:mmse}
\end{equation}

as well as the pixel-wise marginal variance (at the \( i^\text{th} \) pixel)

\begin{equation}
\int {\big( x^\prime_i - \mathbb{E}{[x \mid z]}_i \big)}^2\ p(x^\prime_i\mid z)\ \mathrm{d}x^\prime_i \eqqcolon (\operatorname{Var}[x\mid z])_i.%
\label{eq:variance}
\end{equation}

On the other hand, the representation power is empirically demonstrated by two experiments:
First, our model is capable of synthesizing realistic images without any data.
Second, we show that we can faithfully reconstruct images from severely ill-posed problems, e.g. when using random acquisition masks, where profound knowledge of the underlying anatomy is required.

§.§ Parallel MRI

In modern parallel mri systems, the view of <ref> is overly simplistic:
These systems typically utilize coil arrays with spatially varying sensitivity maps <cit.> which are, in general, entirely unknown.
In other words, the physical acquisition model is not well specified.
Thus, for any real practical reconstruction problem, a great challenge lies in estimating these sensitivity maps accurately.

Many strategies for estimating the sensitivity maps in an off-line step have been proposed (e.g. ESPIRiT <cit.>).
Typically, these require the acquisition of fully sampled acl in the \( k \)-space center such that, in essence, a low-resolution reconstruction along with the corresponding sensitivity maps can be estimated.
The disadvantage of such methods is two-fold:
First, assuming a non-Cartesian sub-sampling patternacquisition protocol, the acquisition of the acl essentially requires an additional (albeit low-resolution) scan.
This increases examination time and opens the possibility of patient movement and subsequent misalignment artifacts.
Second, high-frequency regions of the \( k \)-space are not taken into account for estimating the sensitivity maps.
Generally, the sensitivity maps can be assumed to be smooth (indeed we also make this assumption).
However, since errors in the sensitivity estimation adversely affect downstream tasks significantly, it is vital to exploit all available data to the best extent.

In this work, similar to <cit.>, we propose to estimate the image and sensitivity maps jointly.
This joint estimation hinges on the observation that the sensitivity maps are much smoother than the imaged anatomy.
To enforce this smoothness and simultaneously resolve ambiguities, we impose a simple quadratic penalization on the spatial gradient of the sensitivity maps.
The joint estimation allows us to utilize all available \( k \)-space data to estimate the image as well as the sensitivity maps.
The resulting optimization problem can be efficiently solved by acceleratedmomentum-based non-convex optimization algorithms.
Reconstructing a parallel mri image takes about 5 on consumer hardware.

§.§ Related Work

In this section we review some related work on learned mri reconstruction, but restrict our attention to methods that share similarities with our proposed approach.
For a more general overview of data-driven mri reconstruction we refer to <cit.>.
The authors of <cit.> give a very broad overview of potential risks of using modern techniques for medical image reconstruction in general.

§.§.§ Energy-based Models

After submission of this paper, we were made aware of concurrent works that also utilize ebm for mri reconstruction.
In <cit.>, the authors propose to learn a regularizer and use proximal gradient descent for inference.
The main difference to our approach is that for parallel mri, their algorithm assumes access to sensitivity maps, which have to be precomputed using, e.g., ESPIRiT <cit.>.
Additionally, their algorithm requires hand-tuning of step-sizes and is relatively slow.
On the contrary, we propose an algorithm for joint reconstruction of image and coil sensitivities which does not require hand-tuning of step-sizes, and is fast due to being accelerated.

The authors extended their work in <cit.>, where they learn an ebm for both image- and k-space.
While this obviates the need for sensitivity estimation, it requires to train two independent networks that need to be balanced at inference time.
Additionally, the k-space ebm requires fully-sampled reference k-space data, which is scarcely available.
In contrast, we only require reference DICOM (magnitude) images to train one network.
In addition, we perform a data-independent analysis of the learned regularizer as well as uncertainty analysis.

§.§.§ Diffusion Models

Diffusion models <cit.> aim to model the gradient of the log-prior while undergoing a diffusion process:
Let \( p_t \) denote the data distribution at diffusion time \( t \).
The aim is to learn a time-conditional score network \( s_\theta \) such that \( s_\theta(\,\cdot\,, t) \approx \nabla \log p_t \) for all \( t > 0 \).
This approach shares many similarities with ebm, where indeed it should hold that \( \nabla R = s_\theta(\,\cdot\,, 0) \).

Diffusion models can be used to generate unconditional samples from the data distribution <cit.>.
This is computationally demanding as it requires solving a sde with high accuracy, which typically needs thousand of gradient evaluations <cit.>.
For inverse problems, it is not clear how to optimally incorporate data-fidelity into the sde.
Proposed approaches include data projection <cit.>, annealed Langevin dynamics <cit.> and diffusion posterior sampling <cit.>.
These approaches require hand-tuning of parameters, sometimes at each step of the reverse diffusion <cit.>.
The recent work of <cit.> showed that none of these methods generate samples from the true posterior distribution and propose to augment the score models with normalizing flows.
While their inference algorithm is parameter-free, the approach is opaque due to the introduction of the normalizing flow and is still computationally demanding.

In contrast, ebm enjoy a natural probabilistic interpretation with access to an analytic expression of the posterior.
In addition, map inference does not require solving an sde, but can be done by efficient optimization algorithms.
Having access to the function value (as opposed to only the gradient) also has practical applications, such as being able to inspect the regularization landscape (see <ref>) and utilizing backtracking in optimization algorithms <cit.>.

In the works of <cit.>, parallel imaging is tackled by off-line sensitivity estimation, which comes with the drawbacks outlined in <ref>.
The authors of <cit.> propose to reconstruct individual coil images using a model trained solely on rss reconstructions.
While the results are impressive, the computational cost consequently depends on the number of coils utilized in the physical scanner, and the authors report reconstruction times of up to 10.
In our joint reconstruction algorithm, the gradient of the network is only evaluated once per iteration, irrespective of the number of coils.
Imposing spatial regularity on the coils is extremely fast, where at each iteration we only have to compute very few fast Fourier transforms (see <ref>).

§.§.§ Joint reconstruction

The common approach to parallel imaging is based on a two-step approach:
Coil sensitivities are estimated (typically from acl regions) in an off-line step, after which the image is reconstructed.
A joint reconstruction was first proposed by <cit.> based on alternating minimization, where they explicitly parametrize the sensitivities with low-order polynomials.
In <cit.>, the authors propose an iteratively regularized Gauss-Newton algorithm that enforces spatial smoothness of the sensitivity maps during the iterations.
Their algorithm is not guaranteed to converge for arbitrary initializations, requires hand-tuning of parameters at each update step and leads to noisy solutions when run for too long.
This algorithm was extended to incorporate classical variational penalties, such as the tv <cit.> or the tgv <cit.>, by the authors of <cit.>, in order to suppress noisy solutions.

Indeed, their approach shares a lot of similarities with our proposed approach.
The differences can be summarized as follows:
Instead of hand-crafted regularizers, we employ modern generative learning techniques to learn an expressive regularizer from data, which leads to state-of-the-art reconstructions.
In addition, instead of the iteratively regularized Gauss-Newton algorithm with additional non-trivial sub-problems for the variational penalties, we employ the ipalm algorithm <cit.> for optimization.
Thus, we can guarantee convergence and only require hand-tuning of two parameters (see <ref>).
A sketch of our reconstruction algorithm is shown in <ref>.

For completeness, we mention that the end-to-end variational network of <cit.> also estimates the sensitivities jointly with the image.
However, their approach differs significantly to ours in that they learn a mapping from \( k \)-space to image-space discriminatively, where they utilize the estimated sensitivity maps to enforce data fidelity.
Thus, the network only works well with a particular sub-sampling pattern and coil configuration.
In contrast, we learn image features generatively and impose hand-crafted spatial regularity onto the sensitivity maps, therefore being agnostic to the sub-sampling pattern as well as coil configuration.

§ METHODS

§.§ Maximum-likelihood training

To learn a regularizer such that its induced distribution is indistinguishable from the underlying reference distribution we proceed as follows:
We equip the regularizer with parameters and denote with \( \{ R_\theta : \theta \in \Theta \} \) the family of \( \theta \)-parametrized functions \( \Function{R_\theta}{\set{X}}{\R} \).
\( \Theta \) is a suitably selected set of parameters, \( \set{X} \) is the space of the underlying distribution.
We discuss our particular choice of the \( \theta \)-parametrized family in <ref>.
\( R_\theta \) induces a Gibbs distribution on \( \set{X} \) with density

\begin{equation}
p_\theta(x) \coloneqq \exp(-R_\theta(x))Z^{-1}_\theta,%
\label{eq:gibbs}
\end{equation}

where \( Z_\theta = \int_\set{X}\exp(-R_\theta(\xi))\,\mathrm{d}\xi \) is the partition function.
We find the optimal parameters by maximizing the likelihood of reference data — drawn from a reference distribution \( p_\text{data} \) — under our model:

\begin{equation}
\min_{\theta \in \Theta} \left\{ \Gamma(\theta) \coloneqq \mathbb{E}_{x\sim p_{\text{data}}}[-\log p_{\theta}(x)] \right\}.%
\label{eq:ml}
\end{equation}

The loss function <ref> admits the gradient <cit.> \( \nabla \Gamma(\theta) = \mathbb{E}_{x\sim p_\text{data}} [\nabla_\theta R_\theta(x)] - \mathbb{E}_{x\sim p_\theta} [\nabla_\theta R_\theta(x)] \).

For any interesting regularizer, computing the partition function \( Z_\theta \) is intractable.
Thus, we resort to mcmc techniques to approximate \( \mathbb{E}_{x \sim p_\theta} [\nabla_\theta R_\theta(x)] \).
In particular, following <cit.> we use the ula <cit.>, which iterates

\begin{equation}
x^{(j)} \sim \mathcal{N}\Bigl(x^{(j - 1)} + \frac{\zeta}{2}\nabla \log p_\theta (x^{(j - 1)}), \zeta \mathrm{Id}\Bigr),\ j \in \llbracket 1, J \rrbracket.
\label{eq:langevin}
\end{equation}

Here, \( \mathcal{N}(\mu, \Sigma) \) denotes the normal distribution on \( \set{X} \) with mean \( \mu \) and covariance \( \Sigma \), and \( \llbracket \,\cdot\,,\,\cdot\,\rrbracket \) denotes an integer interval.
Further, \( \zeta > 0 \) is the discretization time step of the associated continuous-time Langevin diffusion sde, which is known to be \( p_\theta \)-stationary <cit.>.
The time discretization biases ula asymptotically, which could be corrected via metropolization.
However, the poor non-asymptotic performance of metropolized ula <cit.> makes this less useful in practice, and we do not use metropolization.
Note that by <ref> \( \nabla \log p_\theta = -\nabla R_\theta \).

In general, maximum-likelihood based learning of ebm is known to be unstable <cit.> due to the slow convergence of ula (<ref>).
To minimize the burn-in time and aid convergence of our sampler, we utilize persistent initialization <cit.>:
Assuming \( p_\theta \) changes only slightly during each update to \( \theta \), samples from previous learning states are good initial guesses for the current iteration.
We implement this idea using a replay buffer that holds samples from previous iterates.
After iterating (<ref>) starting from \( x^{(0)} \) drawn from the replay buffer, we write \( x^{(J)} \) back into the buffer with \( 1 - \pi_{\text{reinit}} \) chance.
Otherwise, we draw a random sample from \( p_{\text{data}} \) to write it to the buffer.
Upon reinitialization, to help mode coverage, we randomly permute the pixels of (on average) every second sample.
We refer to <cit.> for more information of maximum-likelihood training of ebm.

§.§ Reconstruction algorithm

We utilize an image-space SENSE-type <cit.> forward model for parallel mri.
In detail, we relate the real-valued spin density \( x \in \R^n_{\geq 0} \) to the noisy measurement data \( z = {(z_1,\dotsc,z_C)}^\top \in \mathbb{C}^{m} \).
Here, \( m = CK \) combines the \( K \) acquired \( k \)-space data points (not necessarily on a Cartesian grid) of \( C \in \mathbb{N} \) receiver coils.
The forward operator

\begin{equation}
\Function{A}{\R^n}{\mathbb{C}^m}, \quad x \mapsto \begin{pmatrix}
\mathcal{F}_M(\sigma_1 \odot x) \\
\mathcal{F}_M(\sigma_2 \odot x) \\
\vdotswithin{\mathcal{F}_M(\sigma_1 \odot x)} \\
\mathcal{F}_M(\sigma_{C} \odot x)
\end{pmatrix}
\end{equation}

utilizes the sensitivity maps \( {(\sigma_c)}_{c=1}^{C} \in {(\mathbb{C}^n)}^{C} \) of the \( C \) receiver coils, as well as the sampling operator \( \Function{\mathcal{F}_M}{\mathbb{C}^n}{\mathbb{C}^{K}} \) defined by a \( k \)-space trajectory \( M \).

In our model, we assume a real-valued spin-density and empirically demonstrate good performance on the fastMRI dataset <cit.> in <ref>.
In phase-sensitive imaging however, a complex spin-density needs to be assumed.
Our approach can be generalized to this setting by, e.g. splitting real and imaginary channels as in <cit.>, and training the regularizer on this data.

To ease notation, we define \( \Sigma \coloneqq {(\sigma_c)}_{c=1}^C \), \( \mathcal{C} \coloneqq {(\mathbb{C}^n)}^C \), and \( |\,\cdot\,|_{\mathcal{C}} \colon \mathcal{C} \to \R^n_{\geq 0} \), \( {(\sigma_c)}_{c=1}^C \mapsto \sqrt{\sum_{c=1}^C |\sigma_c|^2} \). (\( |\,\cdot\,| \) is the complex modulus acting element-wise on its argument.)
In the case of sub-sampling on a Cartesian grid, \( \mathcal{F}_M \) encodes the Fourier transform followed by the multiplication with a binary mask.
In such sense-type setups, to quantitatively compare to fully-sampled rss reconstructions, the spin density \( x \) has to be re-weighted by \( |\Sigma|_{\mathcal{C}} \) <cit.>:
Denoting the reconstructed image \( u \in \R^n_{\geq 0} \), it is given by
\(
u = x \odot |\Sigma|_{\mathcal{C}},
\)
which motivates an immediate change of variables

\begin{equation}
x = u \oslash |\Sigma|_{\mathcal{C}}.
\label{eq:change of variables}
\end{equation}

Instead of estimating the spin density (and the sensitivity maps), we propose to directly estimate the image \( u \) (and the sensitivity maps).
Formally, given the noisy measurement data \( z = {(z_1,\dotsc, z_C )}^{\top} \)  we propose to find

\begin{equation}
\begin{aligned}
\argmin_{(u, \Sigma) \in \R^n \times \mathcal{C}} H(u, \Sigma) + \delta_{\mathbb{R}_{\geq 0}^n}(u) + \mu F(\Sigma)
\end{aligned}
\label{eq:mri problem}
\end{equation}

\begin{equation}
\begin{aligned}
\Function{H}{\R^n \times \mathcal{C}&}{\R}, \\
(u, {(\sigma_c)}_{c=1}^C) &\mapsto \frac{1}{2} \norm{%
\begin{pmatrix}
\mathcal{F}_M(\sigma_1 \odot u \oslash |\Sigma|_{\mathcal{C}}) \\
\vdotswithin{%
\mathcal{F}_M(\sigma_1 \dot u \slash |\Sigma|_{\mathcal{C}})%
} \\
\mathcal{F}_M(\sigma_{C} \odot u \oslash |\Sigma|_{\mathcal{C}})
\end{pmatrix} - \begin{pmatrix}
z_1 \\ \vdotswithin{y_1} \\ z_C
\end{pmatrix}%
&\phantom{\mapsto}\ + \lambda R (u)
\end{aligned}%
\label{eq:H}
\end{equation}

combines the data fidelity and image-regularization and

\begin{equation}
\begin{aligned}
\Function{F}{\mathcal{C}&}{\R_{\geq 0}}, \\
{(\sigma_c)}_{c=1}^C &\mapsto \frac{1}{2}\sum_{c=1}^{C} \big( \norm{\mathrm{D} \mathfrak{Re}(\sigma_c)}_2^2 + \norm{\mathrm{D} \mathfrak{Im}(\sigma_c)}_2^2 \big)
\end{aligned}
\end{equation}

encodes the smoothness prior on the sensitivity maps.
In the above, \( \delta_{\R^n_{\geq 0}} \) enforces non-negativity on the spin -density using the indicator function

\begin{equation}
\Function{\delta_{\set{A}}}{\R^n}{\{ 0, \infty \}}, \quad x \mapsto \begin{cases} 0 & \text{ if } x \in \set{A}, \\
\infty & \text{ else.}\end{cases}
\end{equation}

\( \Function{\mathrm{D}}{\R^n}{\R^{2n}} \) is the discrete gradient operator such that \( (\mathrm{D} x) \) contains vertically stacked horizontal and vertical first-order finite differences (see, e.g., <cit.>).
It implements physically motivated Dirichlet boundary conditions (\( \sigma_c = 0 \) outside of the image domain \( \forall c \)).
\( \lambda \in \R_+ \) and \( \mu \in \R_+ \) are scalars trading off the strength of the regularization on the image and the sensitivity maps respectively.

Observe that the division by \( |\Sigma|_{\mathcal{C}} \) in <ref> follows from the identification <ref>.
This formulation has the advantage that the sensitivity maps \( {(\sigma_c)}_{c=1}^C \) are implicitly normalized, and we do not have to impose a normalization constraint explicitly during optimization.
Further, the existence of a minimizer is ensured by coercivity of <ref>, which is a direct result of the enforced Dirichlet boundary conditions.

We solve <ref> using the ipalm algorithm <cit.> with Lipschitz backtracking, summarized in <ref>.

The algorithm utilizes the combined fidelity-regularity functional \( H \) defined in (<ref>).
For our learned regularizer, we would have \( R = R_\theta \) as in (<ref>), whereas for a hand crafted regularizer \( R \) may be TV (<ref>).

All gradients are understood in the \( \mathbb{C}\mathbb{R} \)-sense <cit.>.
Recall that the proximal operator \( \Function{\prox_{\alpha G}}{\set{H}}{\set{H}} \) of a proper extended real-valued function \( \Function{G}{\mathcal{H}}{\interval[open left]{-\infty}{\infty}} \) (\( \mathcal{H} \) is a Hilbert space) is the map \( \bar{x} \mapsto \argmin_x \sfrac{1}{2} \norm{\bar{x} - x}^2 + \alpha G(x) \).
Observe that \( \prox_{\delta_{\R^n_+}} \) retrieves the positive part of its argument, that is \( \prox_{\delta_{\R^n_+}}(x) = {(x)}_+ \).
\( \prox_{\mu F} \) can be solved in closed form by utilizing the discrete sine transform (see <cit.> for a more rigorous discussion) as (denoting \( \hat{\imath} \coloneqq \sqrt{-1} \))

\begin{equation}
{(\sigma_c)}_{c=1}^C \mapsto {\big(Q_\mu(\mathfrak{Re}(\sigma_c)) + \hat{\imath} Q_\mu(\mathfrak{Im}(\sigma_c))\big)}_{i=1}^C.
\end{equation}

Here, \( Q_\mu \colon y \mapsto \mathcal{S}^{-1}\bigl({\operatorname{diag}(\xi_i + \mu)}^{-1}\mathcal{S}(\mu y) \bigr) \) uses the discrete sine transform \( \mathcal{S} \) and the eigenvalues \( \xi_i \) of the discrete Laplace operator, which are of the form \( \xi_i = 2 - 2 \cos \phi_i \) for equally spaced angles \( \phi_i \).
We initialize the algorithm with the zf rss reconstruction and the corresponding sensitivity maps (for all \( c \in \llbracket 1, C \rrbracket \))

\begin{equation}
u^{(0)} = \sqrt{\sum_{c=1}^C|\mathcal{F}_M^*(z_c)|^2}\ \text{ and }\ \sigma_c^{(0)} = \mathcal{F}_M^*(z_c) \oslash u^{(0)}.%
\label{eq:initial guess}
\end{equation}

A sketch of the reconstruction algorithm is shown in <ref>.

\( (u^{(0)}, \Sigma^{(0)}) \in \R^n \times \mathcal{C} \), \( K \in \mathbb{N} \), \( \gamma_1 \in (0, 1) \), \( \gamma_2 \in (0, 1) \), initial \( (L_x, L_\Sigma) \in \R_+ \times \R_+ \)
\( (u^{(K)}, \Sigma^{(K)}) \) solving <ref>
\( (u^{(1)}, \Sigma^{(1)}) = (u^{(0)}, \Sigma^{(0)}) \)
\( k \in \llbracket 1, K - 1 \rrbracket \)
\( \bar{u} = u^{(k)} + \frac{k}{k + 3} (u^{(k)} - u^{(k - 1)})\)
\( (u^{(k + 1)}, L_x) = \operatorname{bt}(H(\,\cdot\,,\Sigma^{(k)}), \delta_{\R^n_+}, \bar{u}, L_u, \gamma_1, \gamma_2) \)
\( \bar{\Sigma} = \Sigma^{(k)} + \frac{k}{k + 3} (\Sigma^{(k)} - \Sigma^{(k - 1)}) \)
\( (\Sigma^{(k+1)}, L_\Sigma) = \operatorname{bt}(H(u^{(k + 1)},\,\cdot\,), \mu F, \bar{\Sigma}, L_\Sigma, \gamma_1, \gamma_2) \)

ipalm <cit.> instantiation to solve (<ref>).
\( \operatorname{bt} \) is <ref>.

\( E \), \( P \), \( x_0 \), \( L_0 \), \( \gamma_1 \), \( \gamma_2 \)
\( (x, L) \)
\( L \leftarrow L_0 \)
\( x = \prox_{L^{-1}P}(x_0 - L^{-1}\nabla E(x_0)) \)
\( d = x - x_0 \)

\( E(x)  \leq E(x_0) + \langle \nabla E(x_0), d \rangle + \frac{L}{2} \norm{d}_2^2 \)

\( L \leftarrow \gamma_1 L \)

\( L \leftarrow L/\gamma_2 \)

Backtracking procedure to find the local Lipschitz constants in <ref>.

§.§ Experimental data

For all experiments we utilize the fastMRI knee dataset <cit.>.
Specifically, the training data are the rss reconstructions of size \( \tilde{n} = \num{320} \times \num{320} \) of the multi-coil corpd training split.
We used the central 11 slices to ensure reasonable training data, resulting in a total of 5324 training slices.
To have consistent intensity ranges during training, we normalized each slice by \( x \mapsto \sfrac{x - \min_i x_i}{\norm{x}_\infty - \min_i x_i} \) individually (this normalization was not performed for any of the reference methods).
For validation and testing, we used the multi-coil corpd validation split, discarding samples with width \( w \notin \{ \num{368}, \num{372} \} \), leaving 91 scans.
The scans were split into 30 validation samples and 61 test samples by lexicographic ordering of the filenames.
To be consistent with training, we again restrict our interest to the central 11 slices, resulting in 330 validation slices and 671 test slices.
For the out-of-distribution experiments, we used the central 11 slices of the corpdfs scans (again excluding width \( \notin \{ 368, 372 \} \)) in the fastMRI knee validation dataset.

§.§ Network architecture and implementation details

The family of \( \theta \)-parametrized functions we consider in this work followsworks follow the simple structure

\begin{equation}
R_\theta = U \circ FC \circ S_L \circ S_{L-1} \circ \ldots \circ S_2 \circ S_1 \circ \crop,
\label{eq:reg}
\end{equation}

where \( \Function{\crop}{\R^n}{\R^{\tilde{n}}} \) is a center-crops the validation images of size \( n = \num{640} \times w \) to the size of the training images.
The layer \( S_l \colon x \mapsto \lrelu(\widetilde{W}_l\lrelu(W_l x + b_l) + \tilde{b}_l) \), \( l \in \llbracket 1, L \rrbracket \) utilizes the leaky ReLU \( \lrelu \colon x \mapsto \max\{ \gamma x, x \} \) with leak coefficient \( \gamma = \num{0.05} \).
The weights and biases \( \{ W_l, b_l, \widetilde{W}_l, \tilde{b}_l \} \) in the \( l^{\text{th}} \) layer encode the convolution kernels of size \( \num{3} \times \num{3} \) where \( \widetilde{W}_{l} \) has stride 2.
The number of features in each of the \( L = 6 \) layers follows a geometric progression with common ratio 1.75, starting at 48 features in the first layer.
Finally, \( FC \colon x \mapsto W_{\text{FC}}x \) is a fully connected layer mapping to a scalar and \( U \) is the absolute value.
The entirety of the learnable parameters can be summarized as \( \theta = {\{ W_l, b_l, \widetilde{W}_l, \tilde{b}_l \}}_{l=1}^L \cup \{ W_{\text{FC}} \} \).
We do not impose any constraints on any of the parameters, thus \( \Theta \cong \R^{n_p} \), where \( n_p = \num{21350640} \) is the total number of parameters.

Although our network is quite sizable, it has significantly less learnable parameters than, e.g., the discriminative end-to-end vn of <cit.> (3e7) or the score-based diffusion models of <cit.> (6.7e7).

We optimize <ref> with AdaBelief <cit.> (\( \beta_1 = 0.9, \beta_2 = 0.999 \)).

In contrast to most previous works <cit.> we did not find it necessary to regularize our model by means of, e.g., Lipschitz regularization, magnitude penalization or similar techniques.

We use a learning rate of \( \num{5e-4} \), exponentially decreasing with rate 0.5 at update steps \( \{ \num{500}, \num{2000}, \num{3000}, \num{5000}, \num{7000} \} \), using a batch size of \( 50 \) for 27000 parameter updates.
To stabilize training, we smooth the data distribution by convolving it with a normal distribution of standard deviation \( \num{1.5e-2} \).
For approximating \( \mathbb{E}_{p_\theta} \), we run ula for \( J_{\text{max}} = 500 \) steps.
To accelerate training in the early stages, we use an exponential schedule, detailed by \( J_h = \lceil J_{\text{max}} (1 - \exp(\sfrac{-h}{\num{1000}})) \rceil \), at the \( h^{\text{th}} \) parameter update.

For persistent initialization, we use a replay buffer holding 8000 images with reinitialization chance \( \pi_{\text{reinit}} = \qty{1}{\percent} \).

Further, to minimize the burn-in time of our sampler, we use a replay buffer <cit.> holding 8000 images.
The reinitialization chance of the replay buffer is 1%, and reinitialization happens with an equal chance of uniform noise or samples from the dataset.

Training took approximately one month on a machine equipped with one NVIDIA Quadro RTX 8000.

§.§ Simulation study and posterior sampling

To construct a real-valued simulation study, we use the rss validation and test data detailed in <ref>, and retrospectively sub-sample the corresponding \( k \)-space data.
To approximately map the reconstructions to the same intensities seen during training, we normalize the data by \( z \mapsto \sfrac{z}{\norm{u^{(0)}}_\infty} \), and normalize the final reconstruction by \( u^{*} \mapsto u^{*}\norm{u^{(0)}}_\infty \).
To run inference on this data, we utilize <ref> with one sensitivity map that is fixed at the identity (\( C = 1 \), \( \sigma_1 = \mathbf{1} \) is fixed as the one-vector).
We found the optimal regularization parameter \( \lambda \) by grid search on the validation dataset.

In addition to computing map estimators in the sense of <ref>, we also examine the posterior distribution.
In particular, we approximate the mmse (<ref>) and variance (<ref>) integrals by mcmc:
We run ula on the Gibbs density of \( H \), where again \( C = 1 \) and \( \sigma_1 = \mathbf{1} \) is fixed as the one-vector, i.e. \( p(u \mid z) \propto \exp\big( -H(u, \mathbf{1}) \big) \).
The algorithm is initialized with uniform noise.
We discard the first 10000 samples to reach a steady state, and then save every 1515 iteration, for a total of 160000 iterations (resulting in 10000 saved samples).
We found that the regularization parameter \( \lambda \) barely influenced the results of the posterior sampling and consequently set it to \( \lambda = 1 \) for all sub-sampling patternsmasks.

§.§ Parallel imaging

For the parallel imaging experiments, we run <ref> with \( K = 100 \).
As in the real-valued simulation study, we normalize the data by \( z \mapsto \sfrac{z}{\norm{u^{(0)}}_\infty} \), and normalize the final reconstruction by \( u^{*} \mapsto u^{*}\norm{u^{(0)}}_\infty \).
To find the optimal regularization parameters, we fix \( \mu = 10 \) and obtain \( \lambda \) by linear least-squares regression of the initial residuum \( \sum_{c=1}^C \norm{\mathcal{F}_M(\sigma_{c}^{(0)}(z^{\text{val},i}) \odot u^{(0)}(z^{\text{val},i})) - {(z^{\text{val},i})}_{c}}_2^2 \) against \( \min_{\lambda} \norm{\optimal{u}(z^{\text{val},i}, \lambda) - u^{\text{val},i}}_2^2 \) (found by grid search) for all image-data pairs \( (u^{\text{val},i}, z^{\text{val},i}) \) in the validation set.
Here, we view the initial image \( u^{(0)} \), the sensitivity maps \( {(\sigma_c)}_{c=1}^C \), and the optimal reconstruction \( u^* \) as maps to make dependencies explicit.
This regression is performed only for 4-fold Cartesian sub-sampling with 8 acl and other sub-sampling patternsschemes use the same fit.
Generalization experiments marked with \( \dagger \) use the linear \( \lambda \)-fit calculated on corpd data.
For experiments marked with \( * \), we re-ran the regression on corpdfs data on a 4-fold Cartesian sub-sampling with 8 acl.

A particular characteristic of our reconstruction approach is that itsit intensities are not quantitatively comparable to the reference.
In detail, although we normalize the reconstruction by the rss of the sensitivity maps, we observed that especially in low-intensity regions (e.g. air) the reconstruction did not match the reference.
To remedy this and allow for fair quantitative evaluation, we utilize the validation data to fit a spline curve (cubic splines, 5 equally spaced knots) against the scatter of reconstructed and reference intensities.
For the generalization experiments, we fit the spline curve again on an independent corpdfs validation dataset.

The spline curves for both CORPD and CORPD-FS are shown in <ref> in the appendix.
The insets show that our reconstructions prefer zero-intensity in background regions, whereas the reference images have non-zero background intensity.

For parallel imaging, we define the mmse as follows:
Let \( (\optimal{u}, {(\optimal{\sigma}_c)}_{c=1}^C) \) be a solution to <ref>.
Then, we fix \( {(\optimal{\sigma}_c)}_{c=1}^C \) and sample the conditional probability \( p(u \mid z, {(\optimal{\sigma}_c)}_{c=1}^C ) \).
This amounts to performing Langevin sampling on the Gibbs density of \( H(\, \cdot \,, {(\optimal{\sigma}_c)}_{c=1}^C) \colon \R_{\geq0}^n \to \R \).

The sensitivities may also be included in the Langevin procedure, but we empirically found no noticeable difference to freezing them.
We believe that this is due to the strong imposed spatial regularity.

We evaluate the quality of the estimated sensitivity maps by computing the null-space residual <cit.>:
Let \( u_c = \mathcal{F}^\ast(z_c) \), \( c = 1,\dotsc,C \), denote the fully-sampled coil images.
The null-space residual

\begin{equation}
\pi_c = \frac{\sigma_c}{|\Sigma|^2_{\mathcal{C}}} \sum_{i=1}^{C} \bar{\sigma}_i u_i - u_c,
\label{eqq:null space residual}
\end{equation}

where \( \bar{\,\cdot\,} \) denotes complex conjugation, should only contain noise since \( u_i = \sigma_i u \) when \( \sigma_i \) is exact.
Thus, any residual signal points to sub-optimal sensitivity estimates.

§.§ Comparison and evaluation

We compare our approach to the following methods:
For a hand-crafted prior, we chose the Charbonnier smoothed total variation

\begin{equation}
\text{\gls{tv}}: \R^n_{\geq 0} \ni x \mapsto \sum_{i=1}^n \sqrt{{(\mathrm{D} x)}_{i}^2 + {(\mathrm{D} x)}_{i + n}^2 + \epsilon^2}
\label{eq:charb}
\end{equation}
tv: \( \R^n_{\geq 0} \ni x \mapsto \sum_{i=1}^n \sqrt{{(\mathrm{D} x)}_{i}^2 + {(\mathrm{D} x)}_{i + n}^2 + \epsilon^2} \),
with \( \epsilon = \num{e-3} \).
In the real-valued simulation study, we compare against the fastMRI baseline method <cit.>

as well as the diffusion-based approach of <cit.>.
However, due to time and computational constraints, we limit our comparison to one arbitrarily picked image per sub-sampling pattern.
The implementation as well as the trained model are taken from their github repository, thus the training database differs to ours as it includes the CORPD-FS data.
We use 2000 steps in the reverse diffusion.

As a state-of-the-art discriminative approach for parallel mri, we compare against the end-to-end vn approach from <cit.>.
The implementation was taken from the fastMRI github repository with default parameters.
The fastMRI baseline method as well as the vnmodels were trained on the subset of the fastMRI dataset detailed <ref>, with masks generated using random 4-fold Cartesian sub-sampling with 8 acl.

We compare the reconstructions quantitatively using the psnr, nmse and ssim <cit.>.
ssim uses a \( \num{7} \times \num{7} \) uniform filter and parameters \( K_1 = \num{0.01} \), \( K_2 = \num{0.03} \).
We define the acceleration factor (denoted Acc. in tables showing quantitative results) as the ratio of the image size and the acquired \( k \)-space points on the Cartesian grid.

§ RESULTS

§.§ Data-independent analysis

The local regularization landscape is a proxy for many interesting properties, relating to (local) convexity, generalization capabilities and adversarial robustness <cit.>.
To effectively visualize the \( \tilde{n} = (\num{320} \times \num{320}) \)-dimensional landscape, we follow <cit.>:
Let \( {(x^{(k)})}_{k=1}^{K} \), \( x^{(k)} \in \R^{\tilde{n}} \) be samples drawn from the Langevin process <ref>, starting from uniform noise.
Denote with \( (v_1, v_2) \) the first two principal components of the matrix \( [x^{(0)} - x^{(K)};x^{(1)} - x^{(K)};\dotsc;x^{(K-1)} - x^{(K)} ] \) and denote \( \bar{x} \coloneqq x^{(K)} + \sfrac{1}{K-1} \sum_{k=1}^{K-1} (x^{(k)} - x^{(K)}) \).
For the experiments, we set \( K = \num{10000} \).

In <ref> we show \( \R^2 \ni (\xi_1, \xi_2) \mapsto \exp \big( -R(\bar{x} + \xi_1 v_1 + \xi_2 v_2) / \widetilde{T} \big) \), where \( \widetilde{T} = 7 \) was chosen to yield visually pleasing results, along with the Langevin trajectory.
The figure shows two interesting properties:
First, the samples from the Langevin process are almost indistinguishable from the reference data, which demonstrates that our prior is extremely strong and faithfully represents the distribution of the training data.

In particular, this implies that the learned regularizer — by construction — should only be applied to reconstruction tasks where the underlying data is also drawn from the same distribution.
The experiments in <ref> show that our regularizer can reasonably reconstruct knee images of different contrast by adapting regularization strength, but prior work <cit.> demonstrates that performance quickly degrades when confronted with, e.g., rotated images.

Second, in two dimensions, the landscape appears smooth on the considered domain and almost log-concave around modes of the induced distribution.
Although not applicable directly, these findings should be taken as evidence that the high-dimensional landscape is also reasonably well-behaved.
This is corroborated empirically by the ease of optimization: For all reconstruction tasks, we only need in the order of 10 iterations of ipalm.

[count=] in 0, 1

at (* 13-5, -* 2) ;
[count=] /in
0/"9.49, -2.94", "22.32, -4.76",
20/"10.6, -.86", "19.6, -.84",
50/"10.6, -.86", "19.6, -.84",
200/"10.6, -.86", "19.6, -.84",
500/"7.85, 0.98", "16.85, 0.93",
999/"6.8, -.55", "19.8, -2.80"

((- 2) * 2 + 1) * 3 - 3
(sample) [draw, rectangle, inner sep=2, gray] at (* 9 + * 2.2 - 3 * 2.2, ) ;
=1 =6
(anno) [draw, gray, rectangle, minimum width=0.3cm, minimum height=0.3cm] at () ;
[gray] (sample) – (anno);

[ultra thick] ([xshift=-1cm,yshift=1.3cm]sample21.west) to [out=270,in=90] ([xshift=1cm,yshift=-1.3cm]sample16.east);

Surface plots of the local regularization landscape \( \R^2 \ni (\xi_1, \xi_2) \mapsto \exp(-R(\bar{x} + \xi_1v_1 + \xi_2v_2)/\widetilde{T}) \).
The Langevin trajectory (<ref>) is shown in gold along with some representative states.

§.§ Simulation Study

We show qualitative results of the retrospectively sub-sampled real-valued data in <ref> (left), where the results are shown in increasing acceleration.
In the second row we show the results from 4-fold Cartesian sub-sampling with 8 acl.
Note that this is the setting on which the U-Net was trained, and indeed it yields satisfactory reconstructions.
The tv reconstruction for this task is not able to fully remove the sub-sampling artifacts, but increasing regularization would lead to significant loss of detail.
Our approach is at least on par qualitatively with the discriminative U-Net, and the quantitative analysis in <ref> shows superiority over all reference methods.

The third row details the reconstructions for spiral sub-sampling.
Here, the acceleration factor is approximately \( \num{5} \), with a more densely sampled \( k \)-space center, which is typically advantageous for traditional reconstruction techniques.
The tv reconstruction removes most of the sub-sampling artifacts, although some are still visible, especially in the background.
The discriminative U-Net approach in this task (and indeed all tasks apart from Cartesian sub-sampling, on which it was trained) struggles to discriminate between details in the anatomy and sub-sampling artifacts.
Thus, it hallucinates details into the reconstruction that are not reflected in the data.
In contrast, our approach is able to faithfully reconstruct the knee with no visible artifacts.
The results are similar for the pseudo-radial sub-sampling pattern with 45 spokes shown in the fourth row.

The 3-fold random sub-sampling shown in the first row is particularly interesting.
Here, the \( k \)-space center is not more densely sampled than any other region, which manifests in the zero-filling reconstruction by large-scale intensity shifts.
None of the reference methods are able to correct this, since they do not have knowledge of the anatomy of the human knee.
On the other hand, our approach, due to its generative nature, can restore the general shape of the knee well, and due to the variational approach also faithfully keeps details that are present in the data.
The generative approach alleviates the requirement for sampling schemes to densely sample the \( k \)-space center, which is in stark contrast to the general theory.

We show a quantitative comparison against the diffusion-based approach of <cit.> in <ref> that includes reconstruction time and number of trainable parameters.
The accompanying qualitative results are shown in <ref> in the appendix.
The results are separated since we only evaluated on one image per sub-sampling pattern due to time and computational constraints, and thus it does not constitute a comprehensive comparison.
Our mmse estimate consistently beats the diffusion-based approach, and our map estimate is only inferior for the random sub-sampling pattern.
To generate the 10000 samples for the mmse estimate takes only about 30 more time compared to the one sample generated by the reverse sde.
On the other hand, computing our map estimate is about 80 times faster.

Quantitative results for the real-valued simulation study using different \( k \)-space trajectories.
Numbers in parenthesis indicate the acceleration factor, bold typeface indicates the best method.

2*zf     2*tv     2*U-Net     2cOurs

map     mmse

3*Random (3)
\( \mathparagraph \)     12.931758880615234     20.806852340698242     19.52480125427246     30.8723201751709     32.44416046142578

\( \ddagger \)           100 * 0.6423648595809937     100 * 0.10220915079116821     100 * 0.12830260396003723     100 * 0.013434510678052902     100 * 0.010225472040474415

\( \mathsection \)       0.4817352294921875     0.7321547269821167     0.5721480846405029     0.8576399087905884     0.9027918577194214

3*Cartesian (4)
\( \mathparagraph \)     24.160869598388672     31.469541549682617     34.161155700683594     35.53403091430664     36.168399810791016

\( \ddagger \)           100 * 0.059642061591148376     100 * 0.008760045282542706     100 * 0.004405752755701542     100 * 0.003225221298635006     100 * 0.002839048160240054

\( \mathsection \)       0.704010546207428     0.853325366973877     0.8865730166435242     0.8949872851371765     0.9105123281478882

3*Spiral (\( \approx \num{5} \))
\( \mathparagraph \)     21.210956573486328     31.248300552368164     27.76167869567871     35.35053253173828     36.209625244140625

\( \ddagger \)           100 * 0.1024087592959404     100 * 0.009196978993713856     100 * 0.019011199474334717     100 * 0.0033660277258604765     100 * 0.0028191949240863323

\( \mathsection \)       0.617068350315094     0.84514981508255     0.7843644618988037     0.8829594850540161     0.9014105200767517

3*Radial (\( \approx \num{6} \))
\( \mathparagraph \)     27.015140533447266     32.86052322387695     31.76373863220215     35.04103088378906     35.46697998046875

\( \ddagger \)           100 * 0.03293876722455025     100 * 0.006426490843296051     100 * 0.007553369738161564     100 * 0.0036013659555464983     100 * 0.0033337108325213194

\( \mathsection \)       0.4817352294921875     0.7103990912437439     0.6549839973449707     0.8577146530151367     0.8871539831161499

\( \mathparagraph \): psnr () \( \uparrow \), \( \ddagger \): nmse (\( \times \num{e2} \)) \( \downarrow \), \( \mathsection \): ssim \( \uparrow \)

12.931758880615234, 20.806852340698242, 19.52480125427246, 30.8723201751709, 32.44416046142578,
24.160869598388672, 31.469541549682617, 34.161155700683594, 35.53403091430664, 36.168399810791016,
21.210956573486328, 31.248300552368164, 27.76167869567871, 35.35053253173828, 36.209625244140625,
27.015140533447266, 32.86052322387695, 31.76373863220215, 35.04103088378906, 35.46697998046875

5, 5, 5, 5

[count=] in random, cartesian, spiral, radial

[count=] /in zero_filled/zf, tv/TV, unet/U-Net <cit.>, ours/Ours (map), mean/Ours (mmse), ground_truth/Reference

at (3.1*, -2.75) ;

(onn) at (* 3.1 + , -*4.5 + );
[spy using outlines=rectangle, magnification=3, width=1.48cm, height=1.48cm, connect spies]
at (3.1*, -*4.5) ;
[red] on (onn) in node [left] at (* 3.1 + 1.5, -*4.5-2.25);

[spy using outlines=rectangle, magnification=3, width=1.48cm, height=1.48cm, connect spies]
at (3.1*, -*4.5) ;
[blue] on (onn) in node [left] at (* 3.1-0.01, -*4.5-2.25);

[inner sep=0.1em, outer sep=0.1em, fill=] at (3.1*+0.95, -*4.5+1.15) [round-mode=places,round-precision=2];

at (4, -*4.5-0.9) ;

[thick, black!50] (1.1, -4) – ++(0, -16.5);
[count=] in random, cartesian, spiral, radial

[count=] /in mean/\(\mathbb{E}[x\mid z]\), var/\( \operatorname{Var}[x\mid z] \)
at (3.1*, -2.75) ;
[spy using outlines=rectangle, magnification=3, width=1.48cm, height=1.48cm, connect spies]
at (3.1*, -*4.5) ;
[red] on (* 3.1 + 0.8, -*4.5 + 0.7) in node [left] at (* 3.1 + 1.5, -*4.5-2.25);

[spy using outlines=rectangle, magnification=3, width=1.48cm, height=1.48cm, connect spies]
at (3.1*, -*4.5) ;
[blue] on (* 3.1 + 0.8, -*4.5 + 0.7) in node [left] at (* 3.1 - 0.01, -*4.5-2.25);

Real-valued simulation study results:

1 row:
Random sub-sampling with acceleration factor 3,
2 row:
4-fold Cartesian sub-sampling with 8 acl,
3 row:
Spiral sub-sampling (5-fold acceleration),
4 row:
Radial sub-sampling using 45 spokes.
The inlays show the sub-sampling patternmask (white), a detail zoom (blue), and the magnitude of the difference to the reference (red, 0  0.2).
Left: Reconstruction results.
The numbers show psnr (over the entire test dataset), with the best method emphasized in red.
Right: Uncertainty quantification through marginal posterior variance (\( \num{0}~\protect\drawcolorbar~\num{0.0025} \)).

Quantitative comparison (psnr) against the diffusion-based method of <cit.>.
Bold typeface indicates the best method.

Random     Cartesian     Spiral     Radial     Time (s)     Parameters
2*Ours     map     31.78     34.97     36.07     36.02     3.13     2*2.1e7
mmse     34.76     35.55     37.00     36.56     333.3
<cit.>         33.31     34.65     35.66     34.98     251.6     6.8e7

§.§ Uncertainty Quantification through Posterior Sampling

The natural probabilistic interpretation of our variational approach gives rise to a distribution of reconstructions for any particular reconstruction problem.
We exploit this distribution to compute pixel-wise variance maps as well as the mmse estimate.
The results in <ref> (right) show large variance around small structures, and more variance when less data is available:
The three-fold random sub-sampling shows the least variance, while the approximately six-fold radial sub-sampling shows the most.
The mmse estimator also yields visually pleasing reconstructions.
In fact, <ref> and <ref> establishes quantitative superiority over the map estimate (see <ref> for a possible explanation).

§.§ Parallel Imaging

The first row of <ref> shows a classical Cartesian sub-sampling patternmask, with an acceleration factor of 4 and using 8 acl.
This coincides with the training setup of the discriminative end-to-end variational network approach of <cit.>.
Consequently, it shows the best performance quantitatively as well as qualitatively.
Our method also yields competitive results, achieving a psnr of [round-mode=places,round-precision=2]35.22513961791992.
This is in line with our expectations, as in general we can not expect a generative approach to beat a discriminative counterpart.
However, the strength of our approach becomes apparent when we slightly change the sub-sampling patternacquisition mask:
The performance of the end-to-end vn deteriorates significantly when the phase-encoding direction is swapped, or when less acl are acquired.
We emphasize that between these tasks the number of acquired data points is fixed.
With these minimal changes in the sub-sampling patternacquisition protocol, the method of <cit.> is no longer able to fully remove back-folding artifacts or introduces severe hallucinations, while our method yields comparable results for all three tasks.
This is also reflected in the quantitative evaluation in <ref>, where our method beats the vn approach by a significant margin.

Shifting towards radial and 2D Gaussian sub-sampling patternsmasks, we observe that the vn approach introduces severe artifacts.
We believe that this is a combination of the sensitivity estimation sub-network failing as well as the image sub-network being confronted with previously unseen sub-sampling artifacts.
Indeed, for this task the more general tv approach is superior to the vn, although reconstructions appear over-smoothed.
In line with expectations, our approach reconstructs the image satisfactorily, yielding the best performance quantitatively as well as qualitatively.

Quantitative results for parallel imaging with different sub-sampling patternsmasks on in- and out-of-distribution data.
The \( \dagger \) column shows results using the corpd \( \lambda \)-fit, while the \( * \) column has corpdfs-adapted parameters (see <ref>).

black    2*Acc.     2*acl

black9*Cartesian     9*4     3*8     \( \mathparagraph \)

\( \ddagger \)

\( \mathsection \)

3*8 (horiz.)
\( \mathparagraph \)

\( \ddagger \)

\( \mathsection \)

\( \mathparagraph \)

\( \ddagger \)

\( \mathsection \)

3*Radial     3*\( \approx \num{11} \)     3*—
\( \mathparagraph \)

\( \ddagger \)

\( \mathsection \)

3*2D Gaussian     3*8     3*—     \( \mathparagraph \)

\( \ddagger \)

\( \mathsection \)

5cIn-distribution (corpd)

2*zf     2*tv     2*vn     2cOurs

map     mmse

27.18990707397461     31.867795944213867     36.92045211791992     35.22513961791992     35.27544403076172

100 * 0.0223829485476017     100 * 0.007890134118497372     100 * 0.002435196889564395     100 * 0.0036337128840386868     100 * 0.003618051065132022

0.7375055551528931     0.8079418540000916     0.9163976311683655     0.88604736328125     0.8878628015518188

31.12578582763672     33.03041076660156     24.720752716064453     36.230690002441406     36.01150894165039

100 * 0.00931438896805048     100 * 0.005917737260460854     100 * 0.04011896252632141     100 * 0.0028249204624444246     100 * 0.0029919317457824945

0.8110671043395996     0.8260417580604553     0.670006275177002     0.898271918296814     0.8988670110702515

24.141496658325195     25.81011390686035     32.16292190551758     35.33468246459961     35.224971771240234

100 * 0.04508989304304123     100 * 0.034832630306482315     100 * 0.0069937678053975105     100 * 0.0035170933697372675     100 * 0.0036184450145810843

0.6923396587371826     0.6951988935470581     0.8927095532417297     0.8880932331085205     0.8893769383430481

28.76137351989746     32.46963882446289     20.558855056762695     34.45677947998047     34.1624641418457

100 * 0.015698831528425217     100 * 0.006666331551969051     100 * 0.10128451883792877     100 * 0.004226132296025753     100 * 0.0045379577204585075

0.7465711236000061     0.8061649203300476     0.6916904449462891     0.8565337657928467     0.8555713295936584

32.0957145690918     34.1363525390625     20.741687774658203     35.34770965576172     35.412174224853516

100 * 0.007424597162753344     100 * 0.004535387735813856     100 * 0.09953409433364868     100 * 0.003436274826526642     100 * 0.003413973143324256

0.8387000560760498     0.8521340489387512     0.6813638210296631     0.8809531927108765     0.8852419853210449

5cOut-of-distribution (corpdfs)

2*zf     2*tv     2*vn     2cOurs

\( \dagger \)     \( * \)

26.091176986694336     31.302940368652344     29.998945236206055     30.602453231811523     31.712961196899414

100 * 0.05352313444018364     100 * 0.01478771585971117     100 * 0.023820199072360992     100 * 0.019470086321234703     100 * 0.013532079756259918

0.6760171055793762     0.7300366759300232     0.7710118293762207     0.731391966342926     0.7327521443367004

26.57789421081543     31.557456970214844     28.57035255432129     30.890687942504883     31.652633666992188

100 * 0.051486801356077194     100 * 0.014020383358001709     100 * 0.028979137539863586     100 * 0.022375022992491722     100 * 0.013702597469091415

0.7056003212928772     0.7315869927406311     0.707637369632721     0.7516953945159912     0.7289502024650574

24.9786319732666     29.91111183166504     29.117982864379883     29.933029174804688     31.255582809448242

100 * 0.06647732108831406     100 * 0.020912175998091698     100 * 0.027873607352375984     100 * 0.029169296845793724     100 * 0.015278858132660389

0.6545588970184326     0.7075069546699524     0.7572644352912903      0.7270150184631348     0.7248818874359131

25.058237075805664     31.15333366394043     26.257211685180664     31.433210372924805     31.355825424194336

100 * 0.07370038330554962     100 * 0.015259397216141224     100 * 0.049717504531145096     100 * 0.014507188461720943     100 * 0.01571197621524334

0.6186882853507996     0.712673008441925     0.6861987709999084     0.7252629399299622     0.7009592056274414

26.754220962524414     31.52347755432129     23.461977005004883     31.874441146850586     32.09153366088867

100 * 0.05457334592938423     100 * 0.014161386527121067     100 * 0.09929478168487549     100 * 0.01432089228183031     100 * 0.0124507499858737

0.7052920460700989     0.7458847165107727     0.6833181977272034     0.7555358409881592     0.7461521029472351

\( \mathparagraph \): psnr () \( \uparrow \), \( \ddagger \): nmse (\( \times \num{e2} \)) \( \downarrow \), \( \mathsection \): ssim \( \uparrow \)

27.18990707397461, 31.867795944213867, 36.92045211791992, 35.22513961791992, 35.27544403076172,
31.12578582763672, 33.03041076660156, 24.720752716064453, 36.230690002441406, 36.01150894165039,
24.141496658325195, 25.81011390686035, 32.16292190551758, 35.33468246459961, 35.224971771240234,
28.76137351989746, 32.46963882446289, 20.558855056762695, 34.45677947998047, 34.1624641418457,
32.0957145690918, 34.1363525390625, 20.741687774658203, 35.34770965576172, 35.412174224853516

3, 4, 4, 4, 5

[count=] in cartesian, cartesian_rot90, cartesian_4, radial, gaussian2d

[count=] /in zero_filled/zf, tv/TV, vn/vn <cit.>, ours/Ours (map), ours_mmse/Ours (mmse),ground_truth/Reference

at (3.1*, -2.8) ;

(onn) at (* 3.1 + , -*4.5 + );
[spy using outlines=rectangle, magnification=3, width=1.47cm, height=1.47cm, connect spies]
at (3.1*, -*4.5) ;
[red] on (onn) in node [left] at (* 3.1 + 1.5 - 0.015, -*4.5-2.25);

[spy using outlines=rectangle, magnification=3, width=1.47cm, height=1.47cm, connect spies]
at (3.1*, -*4.5) ;
[blue] on (onn) in node [left] at (* 3.1 - 0.015, -*4.5-2.25);

[inner sep=0.1em, outer sep=0.1em, fill=] at (3.1*+0.95, -*4.5+1.15) [round-mode=places,round-precision=2];

at (0.6, -*4.5) ;

Parallel imaging on in-distribution data:
1 row:
4-fold Cartesian sub-sampling with 8 acl,
2 row:
as previous with swapped phase encoding direction,
3 row:
4-fold Cartesian sub-sampling with 4 acl,
4 row:
Pseudo-radial sampling with 45 equidistant spokes (\( \approx \num{11} \)-fold sub-sampling).
5 row:
2D Gaussian sampling (8-fold sub-sampling).
The inlays show the sub-sampling mask, a detail zoom (blue), and the magnitude of the difference to the reference (red, 0  0.2).

26.091176986694336, 31.302940368652344, 29.998945236206055, 30.602453231811523, 31.712961196899414 ,
26.57789421081543, 31.557456970214844, 28.57035255432129, 30.890687942504883, 31.652633666992188,
24.9786319732666, 29.91111183166504, 29.117982864379883, 29.933029174804688, 31.255582809448242,
25.058237075805664, 31.15333366394043, 26.257211685180664, 31.433210372924805, 31.355825424194336,
26.754220962524414, 31.52347755432129, 23.461977005004883, 31.874441146850586, 32.09153366088867

5, 5, 5, 4, 5

[count=] in cartesian, cartesian_rot90, cartesian_4, radial, gaussian2d

[count=] /in zero_filled/zf, tv/TV, vn/vn <cit.>, ours/Ours \( \dagger \), ours/Ours \( * \),ground_truth/Reference

at (3.1*, -2.8) ;

(onn) at (* 3.1 + , -*4.5 + );
[spy using outlines=rectangle, magnification=3, width=1.47cm, height=1.47cm, connect spies]
at (3.1*, -*4.5) ;
[red] on (onn) in node [left] at (* 3.1 + 1.5 - 0.015, -*4.5-2.25);

[spy using outlines=rectangle, magnification=3, width=1.47cm, height=1.47cm, connect spies]
at (3.1*, -*4.5) ;
[blue] on (onn) in node [left] at (* 3.1 - 0.015, -*4.5-2.25);

[inner sep=0.1em, outer sep=0.1em, fill=] at (3.1*+0.95, -*4.5+1.15) [round-mode=places,round-precision=2];

at (0.6, -*4.5) ;

Parallel imaging on out-of-distribution data (\( \dagger \): regularization parameters from non-fat-suppressed data, \( * \): adapted parameters):
1 row:
4-fold Cartesian sub-sampling with 8 acl,
2 row:
as previous with swapped phase encoding direction,
3 row:
4-fold Cartesian sub-sampling with 4 acl,
4 row:
Pseudo-radial sampling with 45 equidistant spokes (\( \approx \num{11} \)-fold sub-sampling).
5 row:
2D Gaussian sampling (8-fold sub-sampling).
The inlays show the sub-sampling mask, a detail zoom (blue), and the magnitude of the difference to the reference (red, 0  0.2).

To analyze our reconstruction algorithm also with respect to the estimated sensitivity maps, we show example estimations along with ESPIRiT <cit.> in <ref>.
The figure shows estimations for the 4-fold Cartesian sub-sampling with 8 acl reconstruction problem shown in the first row of <ref>.
Due to the re-parametrization <ref>, our reconstruction algorithm does not require pixel-wise normalized sensitivity maps.
Hence, they look very physically plausible and match the reference well.

We additionally analyze the estimation qualitatively by visualizing the rss null-space residual <cit.> \( |(\pi_c)_{c=1}^C|_\mathcal{C} \).
In <ref>, we compare the sensitivity maps from our joint estimation with ESPIRiT for Cartesian sub-sampling patterns with 8 acl and 4 acl (first and third row in <ref> respectively).
For both tasks, the sensitivity maps from our joint estimation algorithm can reproduce the data very well.
In particular, while the ESPIRiT estimation leads to slightly better results when 8 acl are available, the estimation deteriorates for 4.
In contrast, our estimation remains stable, as it can exploit data that is not in the \( k \)-space center.

[count=] /in
ESPIRiT <cit.>/espirit,
\( \frac{\mathcal{F}^\ast(z_c)}{|(\mathcal{F}^\ast(z_c))_{c=1}^C|_\mathcal{C}} \)/fully-sampled

[rotate=90] at (-13.6,-*3.2) ;
at (0, -*3.2) in 0,...,14

Magnitude of the estimated sensitivities using our joint estimation algorithm (top) versus the ESPIRiT <cit.> estimation (middle) for the reconstruction problem shown in <ref> (1 row).
The bottom row shows the reference coil sensitivities computed with the fully-sampled data.

at (-2.1, 4.8) 8 acl;
[rotate=90] at (-4.3, 2.3) Ours;
at (-2.1, 2.3) ;
at (2.1, 4.8) 4 acl;
at (2.1, 2.3) ;
[rotate=90] at (-4.3, -2.3) ESPIRiT <cit.>;
at (-2.1, -2.3) ;
at (2.1, -2.3) ;

[-latex, line width=1mm] (-2.1, 1) – (-2.9, .3);
[-latex, line width=1mm] (-2.1, -3.6) – (-2.9, -4.3);

[-latex, line width=1mm] (1.8, 2.7) – (1.0, 3.3);
[-latex, line width=1mm] (1.8, -1.9) – (1.0, -1.3);

rss null-space residual \(|(\pi_c)_{c=1}^C|_\mathcal{C} \) using our joint estimation algorithm (top) and ESPIRiT <cit.> (bottom).
We consider the two Cartesian sub-sampling reconstruction problems in the first and third row of <ref>: 8 acl (left) and 4 acl (right).

§.§ Generalization

By construction, the learned model encodes the distribution of the training data.
In the previous sections we have demonstrated thoroughly that the prior is agnostic to shifts in the sub-sampling pattern, which is an expected consequence of the generative approach.
However, a natural question is whether the approach is also robust with respect to shifts in the underlying distribution.
To study this, we apply the regularizer to an underlying distribution of fat-suppressed images.

The quantitative and qualitative results in <ref> and <ref> indicate that our method generalizes better to unseen data thanthat the end-to-end vn approach of <cit.>, although the performance degrades significantly in both cases.
To highlight the advantage of the tuneable regularization parameters in our approach, we show results using the \( \lambda \)-fit (see <ref>) calculated on non-fat-suppressed data, as well as using parameters adapted to the task:
The ability to tune the influence of the regularizer leads to improved performance when confronted with previously unseen data.

§ DISCUSSION

<ref> demonstrates that our learned regularizer encodes the training distribution.
As a consequence, we can reconstruct images from severely ill-posed problems (such as random sub-sampling pattern) satisfactorily.
However, this also means that the performance of our learned regularizer is expected to significantly decline when applied to, e.g., different anatomy (rotation is explored in <cit.>).
Diffusion models <cit.> are another instance of models that are capable of generating realistic-looking images from the data-distribution, but for inverse problems it is not clear if this is needed.
Since our network is purely convolutional, it can easily be localized by removing the deeper layers, effectively recovering the fields-of-experts model <cit.>.
This would remove the majority of the trainable parameters (and consequently the training cost), and we hypothesize that it would lead to better generalization as local features are largely shared between different anatomical structures.

The results in <ref> demonstrate that the regularizer can reconstruct images from multi-coil data satisfactorily, irrespective of the acquisition mask.
This is a very significant advantage over the reference methods, since in practice a particular situation might necessitate a different sub-sampling patternacquisition protocol.
As an example, phase encoding direction is often swapped when a blood vessel is located in a particularly disadvantageous position.
In addition, our regularizer does not need re-training when a new sub-sampling pattern is discovered to be advantageous.

The results on multi-coil data also empirically demonstrate that a real-valued regularizer suffices to reconstruct the underlying spin-density, at least for the acquisition sequences used in our dataset <cit.>.
Similarly to <cit.>, our approach could easily be extended to account for complex images if needed.
In particular, this is the case for contrasts that inherently rely on complex-valued information, such as phase-contrast mri.

Generalization is not limited to the acquisition mask:
The regularizer also performs well in out-of-distribution experiments, which highlight the importance of the ability to control its influence.
<ref> shows that adapting the strength of the regularization to the underlying data strongly improves performance.
We note that the regularization strength was always tuned to a 4-fold Cartesian sub-sampling task with 8 acl (see <ref>).
The results for radial sub-sampling indicate that this is a sub-optimal fit for the out-of-distribution data:
The fit calculated on the corpd data (\( \dagger \) column) performs better than the fit calculated on corpdfs data.
<ref> (right) suggests that the adapted parameters lead to an over-smoothed reconstruction.
Clearly, the performance could be improved by adapting the regularization strength to the corpdfs data and the radial sub-sampling.
In general, adapting the regularization strength to acquisition mask as well as the data is greatly beneficial.
However, this is typically not possible with other data-driven approaches.

The accompanying probabilistic interpretation is significant in two ways:
First, the regularizer can be inspected by data-independent analysis as in <ref>, where experts can easily visualize preferred structures.
This aids in interpreting the encoded information and thus improves the confidence with which clinicians may view the reconstruction.
Second, the variance maps provide the clinician with additional information about the inherent uncertainty in the reconstruction.
<ref> (right) clearly shows high variance around small anatomic structures.
If the clinician's decision making is based on such regions, they might deem it necessary to acquire additional information.

In the single-coil simulation study, quantitative and qualitative analysis shows superiority of the mmse estimate over the map estimate.
We believe that this is related to the training procedure:
During training, the regularizer is only ever confronted with slightly noisy images in the Langevin process (see <ref>), and injecting the same noise in the reconstruction algorithm improves overall performance.
We are unsure why this performance improvement does not translate to the parallel imaging experiments, but hypothesize that the our joint estimation biases the reconstruction such that this effect is no longer observable.

Finally, we emphasize that our algorithm for reconstructing parallel mri is fast:
<ref> converges in around 100 iterations, which takes around 5 on an NVIDIA Titan RTX using approximately 2 of memory.
This is in stark contrast to score-based diffusion models, where the reverse diffusion process typically has to be solved with high accuracy and consequently reconstruction time is in the order of 10 <cit.>.
Speed is advantageous in practice, as images can be viewed while the patient is still in the scanner.
This allows fast changes to the sub-sampling patternacquisition protocol, should they be needed.

§ CONCLUSION

We utilize modern generativeunsupervised learning techniques to train a regularizer such that it encodes the underlying distribution of the training data faithfully.
By embedding the regularizer in a variational reconstruction framework, we can satisfactorily reconstruct single-coil and parallel mri by adapting the data likelihood term.
Quantitative and qualitative analysis indicate competitive reconstruction performance, on-par with (or sometimes superior to) fully supervised methods.
In addition, our approach is agnostic to changes in the acquisition mask, while the supervised reference methods introduce severe hallucinations when confronted with previously unseen data.
The knowledge of the human anatomy encoded in the regularizer even allows us to reconstruct high quality images from random sampling masks.
Finally, on out-of-distribution data, our method is still on-par or better than hand-crafted regularizers such as the tv or other supervised methods.
This indicates that our method also generalizes better with respect to shifts in the underlying distribution.

In addition to a competitive reconstruction performance, our method has a natural accompanying probabilistic interpretation through statistical modeling and Bayesian inference.
This allows experts to explore a distribution of reconstructions, while other data-driven approaches only act as point estimators.
Our experiments indicate that this distribution encodes important diagnostic information, such as high uncertainty around small anatomical structures.
We believe that this additional information can aid the practitioner's decision making.
Further, we can perform a data-independent analysis of the model by visualizing the induced distribution.
This allows to interpret the information encoded in the regularizer, whereas other data-driven methods solely act as black-boxes in this respect.

To reconstruct parallel mri, we propose a fast algorithm to jointly estimate the image and sensitivity maps.
In contrast to off-line sensitivity estimation, our approach does not require auto-calibration lines, can easily be applied to non-Cartesian sampling trajectories, and uses all available data.
The resulting optimization problem can be solved efficiently using non-convex optimization algorithms.
High-fidelity reconstructions are available in approximately 5 on consumer hardware.

In general, we believe that reconstruction approaches based on generative priors have huge potential.
A natural extension of our work would be to combine the image-prior with a learned sensitivity-prior.
In addition, future research could focus on extending the simple architecture used in this paper with more modern building blocks, such as attention layers.

A related direction is to research if local convolutional models can replicate the performance of our generative prior.

§ APPENDIX

§.§ Spline fit

<ref> shows the scatter plot of reconstructed versus reference intensities along with the spline fit discussed in <ref>.

Spline fit computed on a validation set for CORPD (top) and CORPD-FS data (bottom).
The figures show the log-histogram of reconstructed versus reference intensities with the insets show the region around zero, where they deviate the strongest.

§.§ Additional results

A qualitative comparison against the diffusion-based method of <cit.> is presented in <ref>.
The accompanying quantitative evaluation is shown in <ref> in the main body.

[count=] in random, cartesian, spiral, radial

[count=] /in recon/ <cit.>, ours/Ours (map), mean/Ours (mmse), ground_truth/Reference

=1at (2.05*, -1.8) ;
(onn) at (* 2.05 + , -*3 + );
[spy using outlines=rectangle, magnification=3, width=0.98cm, height=0.98cm, connect spies]
at (2.05*, -*3) ;
[red] on (onn) in node [left] at (* 2.05 + 1, -*3-1.5);

[spy using outlines=rectangle, magnification=3, width=0.98cm, height=0.98cm, connect spies]
at (2.05*, -*3) ;
[blue] on (onn) in node [left] at (* 2.05-0.01, -*3-1.5);

at (2.65 + 2.05*3, -*3-.6) ;

Qualitative comparison against <cit.> on in-distribution data:
1 row:
Random sub-sampling with acceleration factor 3,
2 row:
4-fold Cartesian sub-sampling with 8 acl,
3 row:
Spiral sub-sampling (5-fold acceleration),
4 row:
Radial sub-sampling using 45 spokes.
The inlays show the sub-sampling pattern (white), a detail zoom (blue), and the magnitude of the difference to the reference (red, 0  0.2).

§ ACKNOWLEDGMENT

We thank the reviewers for their time and effort spent on reviewing the manuscript.
Their comments and suggestions helped us to greatly improve the quality the manuscript.

Missing 'biblatex' package
The bibliography requires the 'biblatex' package.

journaltitleMagn. Reson. Med.
titleSENSE: Sensitivity Encoding for Fast MRI

journaltitleMagn. Reson. Med.
titleComprehensive quantification of signal-to-noise ratio and g-factor for image-based and k-space-based parallel imaging reconstructions

journaltitleIEEE Trans. Inf. Theory
titleCompressed Sensing

journaltitleMagn. Reson. Med.
titleParallel imaging with nonlinear reconstruction using variational penalties

journaltitleProc. CVPR
titleStatistics of natural images and models
pages541547 Vol. 1

journaltitleMagn. Reson. Med.
titleLearning a variational network for reconstruction of accelerated MRI data

journaltitleMagn. Reson. Med.
titleScan-specific robust artificial-neural-networks for k-space interpolation (RAKI) reconstruction: Database-free deep learning for fast imaging

journaltitleArXiv e-prints
titlefastMRI: An Open Dataset and Benchmarks for Accelerated MRI

Curran Associates, Inc.

booktitleAdv. Neural Inf. Process. Syst.
titleInvert to Learn to Invert

booktitleProc. CVPR
titleDuDoRNet: Learning a Dual-Domain Recurrent Network for Fast MRI Reconstruction with Deep T1 Prior

booktitleWavelets and Sparsity XVIII
titleInverse GANs for Accelerated MRI Reconstruction

abstractScore-based diffusion models provide a powerful way to model images using the gradient of the data distribution. Leveraging the learned score function as a prior, here we introduce a way to sample data from a conditional distribution given the measurements, such that the model can be readily used for solving inverse problems in imaging, especially for accelerated MRI. In short, we train a continuous time-dependent score function with denoising score matching. Then, at the inference stage, we iterate between the numerical SDE solver and data consistency step to achieve reconstruction. Our model requires magnitude images only for training, and yet is able to reconstruct complex-valued data, and even extends to parallel imaging. The proposed method is agnostic to sub-sampling patterns and has excellent generalization capability so that it can be used with any sampling schemes for any body parts that are not used for training data. Also, due to its generative nature, our approach can quantify uncertainty, which is not possible with standard regression settings. On top of all the advantages, our method also has very strong performance, even beating the models trained with full supervision. With extensive experiments, we verify the superiority of our method in terms of quality and practicality.
journaltitleMed. Image Anal.
titleScore-based diffusion models for accelerated MRI
Score-based models,Diffusion models,Inverse problems,MRI

journaltitleIEEE Trans. Pattern Anal. Mach. Intell.
titleTrainable Nonlinear Reaction Diffusion: A Flexible Framework for Fast and Effective Image Restoration

Springer International Publishing

booktitleLect. Notes Comput. Sci.
titleVariational Networks: Connecting Variational Methods and Deep Learning

Springer International Publishing

abstractMagnetic resonance imaging (MRI) is known to be a slow imaging modality and undersampling in k-space has been used to increase the imaging speed. However, image reconstruction from undersampled k-space data is an ill-posed inverse problem. Iterative algorithms based on compressed sensing have been used to address the issue. In this work, we unroll the iterations of the primal-dual hybrid gradient algorithm to a learnable deep network architecture, and gradually relax the constraints to reconstruct MR images from highly undersampled k-space data. The proposed method combines the theoretical convergence guarantee of optimization methods with the powerful learning capability of deep networks. As the constraints are gradually relaxed, the reconstruction model is finally learned from the training data by updating in k-space and image domain alternatively. Experiments on in vivo MR data demonstrate that the proposed method achieves superior MR reconstructions from highly undersampled k-space data over other state-of-the-art image reconstruction methods.
booktitleProc. MICCAI
titleModel Learning: Primal Dual Networks for Fast MR Imaging

Springer Science
Business Media LLC

journaltitleNat. Mach. Intell.
titleDeep variational network for rapid 4D flow MRI reconstruction

Springer Science
Business Media LLC

issn0028-0836, 1476-4687
titleImage Reconstruction by Domain-Transform Manifold Learning
otero/storage/2XT8LY8H/Zhu et al. - 2018 - Image reconstruction by domain-transform manifold .pdf:application/pdf

Society for Industrial & Applied Mathematics (SIAM)

journaltitleSIAM J. Imaging Sci.
titleInertial Proximal Alternating Linearized Minimization (iPALM) for Nonconvex and Nonsmooth Problems

booktitleProc. MICCAI
titleEnd-to-End Variational Networks for Accelerated MRI Reconstruction

abstractThe goal of compressed sensing is to estimate a vector from an underdetermined system of noisy linear measurements, by making use of prior knowledge on the structure of vectors in the relevant domain. For almost all results in this literature, the structure is represented by sparsity in a well-chosen basis. We show how to achieve guarantees similar to standard compressed sensing but without employing sparsity at all. Instead, we suppose that vectors lie near the range of a generative model $G: \mathbb{R}^k \to \mathbb{R}^n$. Our main theorem is that, if $G$ is $L$-Lipschitz, then roughly $\mathcal{O}(k \log L)$ random Gaussian measurements suffice for an $\ell_2/\ell_2$ recovery guarantee. We demonstrate our results using generative models from published variational autoencoder and generative adversarial networks. Our method can use $5$-$10$x fewer measurements than Lasso for the same accuracy.
booktitleProc. ICML
titleCompressed Sensing using Generative Models

booktitleProc. CVPR
titleDeep Image Prior

booktitleInternational Conference on Learning Representations
titleSolving Inverse Problems in Medical Imaging with Score-Based Generative Models

journaltitleMagnetic Resonance in Medicine
titleBayesian MRI reconstruction with joint uncertainty estimation using diffusion models

titleScore-Based Diffusion Models as Principled Priors for Inverse Imaging

Proceedings of the National Academy of Sciences

journaltitleProc. Natl. Acad. Sci.
titleOn instabilities of deep learning in image reconstruction and the potential costs of AI

booktitleProc. CVPR
titleTotal Deep Variation for Linear Inverse Problems

Curran Associates, Inc.

abstractEnergy based models (EBMs) are appealing due to their generality and simplicity in likelihood modeling, but have been traditionally difﬁcult to train. We present techniques to scale MCMC based EBM training on continuous neural networks, and we show its success on the high-dimensional data domains of ImageNet32x32, ImageNet128x128, CIFAR-10, and robotic hand trajectories, achieving better samples than other likelihood models and nearing the performance of contemporary GAN approaches, while covering all modes of the data. We highlight some unique capabilities of implicit generation such as compositionality and corrupt image reconstruction and inpainting. Finally, we show that EBMs are useful models across a wide variety of tasks, achieving state-of-the-art out-of-distribution classiﬁcation, adversarially robust classiﬁcation, state-of-the-art continual online class learning, and coherent long term predicted trajectory rollouts.
booktitleAdv. Neural Inf. Process. Syst.
journaltitlearXiv:1903.08689 [cs, stat]
titleImplicit Generation and Modeling with Energy Based Models
u und Mordatch - 2019 - Implicit Generation and Generalization in Energy-B.pdf:application/pdf
Computer Science - Computer Vision and Pattern Recognition,Computer Science - Machine Learning,Statistics - Machine Learning

abstractThis study investigates the effects of Markov chain Monte Carlo (MCMC) sampling in unsupervised Maximum Likelihood (ML) learning. Our attention is restricted to the family of unnormalized probability densities for which the negative log density (or energy function) is a ConvNet. We ﬁnd that many of the techniques used to stabilize training in previous studies are not necessary. ML learning with a ConvNet potential requires only a few hyper-parameters and no regularization. Using this minimal framework, we identify a variety of ML learning outcomes that depend solely on the implementation of MCMC sampling.
journaltitlearXiv:1903.12370 [cs, stat]
titleOn the Anatomy of MCMC-Based Maximum Likelihood Learning of Energy-Based Models
Computer Science - Computer Vision and Pattern Recognition,Computer Science - Machine Learning,Statistics - Machine Learning

San Diego, CA, USA

abstractWe develop a framework for learning generic, expressive image priors that capture the statistics of natural scenes and can be used for a variety of machine vision tasks. The approach extends traditional Markov Random Field (MRF) models by learning potential functions over extended pixel neighborhoods. Field potentials are modeled using a Products-of-Experts framework that exploits nonlinear functions of many linear ﬁlter responses. In contrast to previous MRF approaches all parameters, including the linear ﬁlters themselves, are learned from training data. We demonstrate the capabilities of this Field of Experts model with two example applications, image denoising and image inpainting, which are implemented using a simple, approximate inference scheme. While the model is trained on a generic image database and is not tuned toward a speciﬁc application, we obtain results that compete with and even outperform specialized techniques.
booktitleIEEE Conference on Computer Vision and Pattern Recognition
eventtitle2005 IEEE Computer Society Conference on Computer Vision and Pattern Recognition (CVPR'05)
shorttitleFields of Experts
titleFields of Experts: A Framework for Learning Image Priors
pages860867 vol. 2
oth und Black - 2005 - Fields of Experts A Framework for Learning Image .pdf:application/pdf

journaltitleNMR in Biomedicine
titleMagnetic resonance imaging reconstruction using a deep energy-based model

abstractDecreasing magnetic resonance (MR) image acquisition times can potentially make MR examinations more accessible. Prior arts including the deep learning models have been devoted to solving the problem of long MRI imaging time. Recently, deep generative models have exhibited great potentials in algorithm robustness and usage flexibility. Nevertheless, none of existing schemes can be learned from or employed to the k-space measurement directly. Furthermore, how do the deep generative models work well in hybrid domain is also worth being investigated. In this work, by taking advantage of the deep energy-based models, we propose a k-space and image domain collaborative generative model to comprehensively estimate the MR data from under-sampled measurement. Equipped with parallel and sequential orders, experimental comparisons with the state-of-the-arts demonstrated that they involve less error in reconstruction accuracy and are more stable under different acceleration factors.
journaltitleMagnetic Resonance Imaging
titleK-space and image domain collaborative energy-based model for parallel MRI reconstruction
Parallel magnetic resonance imaging,Image reconstruction,Energy-based model,Hybrid domain

journaltitleMagn. Reson. Med.
titleESPIRiT-an eigenvalue approach to autocalibrating parallel MRI: Where SENSE meets GRAPPA

journaltitleMagn. Reson. Med.
titleSelf-calibrating parallel imaging with automatic coil sensitivity extraction

abstractAbstract The use of parallel imaging for scan time reduction in MRI faces problems with image degradation when using GRAPPA or SENSE for high acceleration factors. Although an inherent loss of SNR in parallel MRI is inevitable due to the reduced measurement time, the sensitivity to image artifacts that result from severe undersampling can be ameliorated by alternative reconstruction methods. While the introduction of GRAPPA and SENSE extended MRI reconstructions from a simple unitary transformation (Fourier transform) to the inversion of an ill-conditioned linear system, the next logical step is the use of a nonlinear inversion. Here, a respective algorithm based on a Newton-type method with appropriate regularization terms is demonstrated to improve the performance of autocalibrating parallel MRI—mainly due to a better estimation of the coil sensitivity profiles. The approach yields images with considerably reduced artifacts for high acceleration factors and/or a low number of reference lines. Magn Reson Med, 2008. © 2008 Wiley-Liss, Inc.
journaltitleMagn. Reson. Med.
titleImage reconstruction by regularized nonlinear inversion—Joint estimation of coil sensitivities and image content

journaltitleMagn. Reson. Med.
titleJoint image reconstruction and sensitivity estimation in SENSE (JSENSE)

Springer Science
Business Media LLC

journaltitleBMC Med. Imag.
titleA review on deep learning MRI reconstruction without fully sampled k-space

booktitleProc. ICLR
titleScore-Based Generative Modeling through Stochastic Differential Equations

Curran Associates, Inc.

booktitleAdv. in Neural Inf. Proc. Syst.
titleRobust Compressed Sensing MRI with Deep Generative Priors

booktitleProc. ICLR
titleDiffusion Posterior Sampling for General Noisy Inverse Problems

abstractA constrained optimization type of numerical algorithm for removing noise from images is presented. The total variation of the image is minimized subject to constraints involving the statistics of the noise. The constraints are imposed using Lanrange multipliers. The solution is obtained using the gradient-projection method. This amounts to solving a time dependent partial differential equation on a manifold determined by the constraints. As t the solution converges to a steady state which is the denoised image. The numerical algorithm is simple and relatively fast. The results appear to be state-of-the-art for very noisy images. The method is noninvasive, yielding sharp edges in the image. The technique could be interpreted as a first step of moving each level set of the image normal to itself with velocity equal to the curvature of the level set divided by the magnitude of the gradient of the image, and a second step which projects the image back onto the constraint set.
journaltitlePhysica D: Nonlinear Phenomena
shortjournalPhysica D: Nonlinear Phenomena
titleNonlinear total variation based noise removal algorithms
BDR9D7X/Rudin et al. - 1992 - Nonlinear total variation based noise removal algo.pdf:application/pdf

Society for Industrial & Applied Mathematics (SIAM)

abstractThe novel concept of total generalized variation of a function u is introduced, and some of its essential properties are proved. Diﬀerently from the bounded variation seminorm, the new concept involves higher-order derivatives of u. Numerical examples illustrate the high quality of this functional as a regularization term for mathematical imaging problems. In particular this functional selectively regularizes on diﬀerent regularity levels and, as a side eﬀect, does not lead to a staircasing eﬀect.
journaltitleSIAM Journal on Imaging Sciences
shortjournalSIAM J. Imaging Sci.
titleTotal Generalized Variation
redies et al. - 2010 - Total Generalized Variation.pdf:application/pdf

MIT Press - Journals

issn0899-7667, 1530-888X
journaltitleNeural Comput.
shortjournalNeural Comput.
titleTraining Products of Experts by Minimizing Contrastive Divergence
ZDGFCG/Hinton - 2002 - Training Products of Experts by Minimizing Contras.pdf:application/pdf

Bernoulli Society for Mathematical Statistics

titleExponential Convergence of Langevin Distributions and their Discrete Approximations
Diffusions,discrete approximations,geometric ergodicity,Hastings algorithms,irreducible Markov processes,Langevin models,Markov chain Monte Carlo,Metropolis algorithms,posterior distributions

journaltitleJ. R. Stat. Soc. (B)
titleOptimal scaling of discrete approximations to Langevin diffusions

abstractModern imaging methods rely strongly on Bayesian inference techniques to solve challenging imaging problems. Currently, the predominant Bayesian computation approach is convex optimization, which scales very efficiently to high-dimensional image models and delivers accurate point estimation results. However, in order to perform more complex analyses, for example, image uncertainty quantification or model selection, it is necessary to use more computationally intensive Bayesian computation techniques such as Markov chain Monte Carlo methods. This paper presents a new and highly efficient Markov chain Monte Carlo methodology to perform Bayesian computation for high-dimensional models that are log-concave and nonsmooth, a class of models that is central in imaging sciences. The methodology is based on a regularized unadjusted Langevin algorithm that exploits tools from convex analysis, namely, Moreau–Yoshida envelopes and proximal operators, to construct Markov chains with favorable convergence properties. In addition to scaling efficiently to high-dimensions, the method is straightforward to apply to models that are currently solved by using proximal optimization algorithms. We provide a detailed theoretical analysis of the proposed methodology, including asymptotic and nonasymptotic convergence results with easily verifiable conditions, and explicit bounds on the convergence rates. The proposed methodology is demonstrated with four experiments related to image deconvolution and tomographic reconstruction with total-variation and $ℓ_1$ priors, where we conduct a range of challenging Bayesian analyses related to uncertainty quantification, hypothesis testing, and model selection in the absence of ground truth.
journaltitleSIAM J. Imaging Sci.
titleEfficient Bayesian Computation by Proximal Markov Chain Monte Carlo: When Langevin Meets Moreau

Kluwer Academic Publishers

abstractStochastic approximation methods play a central role in maximum likelihood estimation problems involving intractable likelihood functions, such as marginal likelihoods arising in problems with missing or incomplete data, and in parametric empirical Bayesian estimation. Combined with Markov chain Monte Carlo algorithms, these stochastic optimisation methods have been successfully applied to a wide range of problems in science and industry. However, this strategy scales poorly to large problems because of methodological and theoretical difficulties related to using high-dimensional Markov chain Monte Carlo algorithms within a stochastic approximation scheme. This paper proposes to address these difficulties by using unadjusted Langevin algorithms to construct the stochastic approximation. This leads to a highly efficient stochastic optimisation methodology with favourable convergence properties that can be quantified explicitly and easily checked. The proposed methodology is demonstrated with three experiments, including a challenging application to statistical audio analysis and a sparse Bayesian logistic regression with random effects problem.
journaltitleStat. Comput.
titleEfficient Stochastic Optimisation by Unadjusted Langevin Monte Carlo: Application to Maximum Marginal Likelihood and Empirical Bayesian Estimation
Stochastic approximation,Markov chain Monte Carlo methods,Maximum marginal likelihood estimation,Unadjusted Langevin Algorithm,Recursive estimation,Empirical Bayesian inference

booktitleProc. ICML
titleImproved Contrastive Divergence Training of Energy Based Models

Springer Science
Business Media LLC

abstractThis article presents a statistical theory for texture modeling. This theory combines ﬁltering theory and Markov random ﬁeld modeling through the maximum entropy principle, and interprets and clariﬁes many previous concepts and methods for texture analysis and synthesis from a uniﬁed point of view. Our theory characterizes the ensemble of images I with the same texture appearance by a probability distribution f (I) on a random ﬁeld, and the objective of texture modeling is to make inference about f (I), given a set of observed texture examples. In our theory, texture modeling consists of two steps. (1) A set of ﬁlters is selected from a general ﬁlter bank to capture features of the texture, these ﬁlters are applied to observed texture images, and the histograms of the ﬁltered images are extracted. These histograms are estimates of the marginal distributions of f (I). This step is called feature extraction. (2) The maximum entropy principle is employed to derive a distribution p(I), which is restricted to have the same marginal distributions as those in (1). This p(I) is considered as an estimate of f (I). This step is called feature fusion. A stepwise algorithm is proposed to choose ﬁlters from a general ﬁlter bank. The resulting model, called FRAME (Filters, Random ﬁelds And Maximum Entropy), is a Markov random ﬁeld (MRF) model, but with a much enriched vocabulary and hence much stronger descriptive ability than the previous MRF models used for texture modeling. Gibbs sampler is adopted to synthesize texture images by drawing typical samples from p(I), thus the model is veriﬁed by seeing whether the synthesized texture images have similar visual appearances to the texture images being modeled. Experiments on a variety of 1D and 2D textures are described to illustrate our theory and to show the performance of our algorithms. These experiments demonstrate that many textures which are previously considered as from different categories can be modeled and synthesized in a common framework.
journaltitleInternational Journal of Computer Vision
titleFilters, Random Fields and Maximum Entropy (FRAME): Towards a Unified Theory for Texture Modeling
otero/storage/DFUJXDL4/Zhu et al. - 1998 - Filters, Random Fields and Maximum Entropy (FRAME).pdf:application/pdf

Helsinki, Finland

Association for Computing Machinery

abstractA new algorithm for training Restricted Boltzmann Machines is introduced. The algorithm, named Persistent Contrastive Divergence, is different from the standard Contrastive Divergence algorithms in that it aims to draw samples from almost exactly the model distribution. It is compared to some standard Contrastive Divergence and Pseudo-Likelihood algorithms on the tasks of modeling and classifying various types of data. The Persistent Contrastive Divergence algorithm outperforms the other algorithms, and is equally fast and simple.
booktitleProc. ICML
titleTraining Restricted Boltzmann Machines Using Approximations to the Likelihood Gradient

abstractA sensing array for color imaging includes individual luminance- and chrominance-sensitive elements that are so intermixed that each type of element (i.e., according to sensitivity characteristics) occurs in a repeated pattern with luminance elements dominating the array. Preferably, luminance elements occur at every other element position to provide a relatively high frequency sampling pattern which is uniform in two perpendicular directions (e.g., horizontal and vertical). The chrominance patterns are interlaid therewith and fill the remaining element positions to provide relatively lower frequencies of sampling. In a presently preferred implementation, a mosaic of selectively transmissive filters is superposed in registration with a solid state imaging array having a broad range of light sensitivity, the distribution of filter types in the mosaic being in accordance with the above-described patterns.
issn0962-4929, 1474-0508
journaltitleActa Numer.
shortjournalActa Numer.
titleAn introduction to continuous optimization for imaging

titleThe Complex Gradient Operator and the CR-Calculus
Optimization and Control (math.OC),Complex Variables (math.CV),FOS: Mathematics,FOS: Mathematics

Cambridge University Press

titleNumerical Recipes in C (2nd Ed.): The Art of Scientific Computing

journaltitleRadiol. Artif. Intell.
titlefastMRI: A Publicly Available Raw k-Space and DICOM Dataset of Knee Images for Accelerated MR Image Reconstruction Using Machine Learning

journaltitleAdv. Neural Inf. Process. Syst.
titleAdaBelief Optimizer: Adapting Stepsizes by the Belief in Observed Gradients

journaltitleIEEE Trans. Image Process.
titleImage quality assessment: from error visibility to structural similarity

booktitleProc. ICCV
titleRelating Adversarially Robust Generalization to Flat Minima

Curran Associates, Inc.

booktitleAdv. Neural Inf. Process. Syst.
titleVisualizing the Loss Landscape of Neural Nets

Verlag der Technischen Universität Graz

abstractIn the past decades, Computed Tomography (CT) has established itself as one of the most important imaging techniques in medicine. Today, the applicability of CT is only limited by the deposited radiation dose, reduction of which manifests in noisy or incomplete measurements. Thus, the need for robust reconstruction algorithms arises. In this work, we learn a parametric regularizer with a global receptive field by maximizing it’s likelihood on reference CT data. Due to thisunsupervised learning strategy, our trained regularizer truly represents higher-level domain statistics, which we empirically demonstrate by synthesizing CT images. Moreover, this regularizer can easily be applied to different CT reconstruction problems by embedding it in a variational framework, whichincreases flexibility and interpretability compared to feedforward learning-based approaches. In addition, the accompanying probabilistic perspective enables experts to explore the full posterior distribution and may quantify uncertaintyof the reconstruction approach. We apply the regularizer to limited-angle and few-view CT reconstruction problems, where it outperforms traditional reconstruction algorithms by a large margin.
booktitleProc. of the OAGM Workshop
titleComputed Tomography Reconstruction Using Generative Energy-Based Priors