Stochastic Deep Restoration Priors for Imaging Inverse Problems

## Abstract

Abstract Deep neural networks trained as image denoisers are widely used as priors for solving imaging inverse problems. While Gaussian denoising is thought sufficient for learning image priors, we show that priors from deep models pre-trained as more general restoration operators can perform better. We introduce Stochastic deep Restoration Priors (ShaRP) , a novel method that leverages an ensemble of such restoration models to regularize inverse problems. ShaRP improves upon methods using Gaussian denoiser priors by better handling structured artifacts and enabling self-supervised training even without fully sampled data. We prove ShaRP minimizes an objective function involving a regularizer derived from the score functions of minimum mean square error (MMSE) restoration operators, and theoretically analyze its convergence. Empirically, ShaRP achieves state-of-the-art performance on tasks such as magnetic resonance imaging reconstruction and single-image super-resolution, surpassing both denoiser- and diffusion-model-based methods without requiring retraining.

## 1 Introduction

Many problems in computational imaging, biomedical imaging, and computer vision can be viewed as *inverse problems*, where the goal is to recover an unknown image from its noisy and incomplete measurements. Inverse problems are typically ill-posed, thus requiring additional prior information for accurate image reconstruction. While many approaches have been proposed for implementing image priors, the current research focuses on methods based on deep learning (DL) (McCann et al., ; Ongie et al., ; Kamilov et al., ; Wen et al., ).

Deep neural networks trained as image denoisers are widely-used for specifying image priors for solving *general* inverse problems (Romano et al., ; Kadkhodaie & Simoncelli, ; Zhang et al., ). The combination of pre-trained Gaussian denoisers with measurement models has been shown to be effective in many inverse problems, including image super-resolution, deblurring, and medical imaging (Metzler et al., ; Zhang et al., ; Meinhardt et al., ; Dong et al., ; Zhang et al., ; Wei et al., ; Zhang et al., ) (see also the recent reviews
(Ahmad et al., ; Kamilov et al., )). This success has led to active research on novel methods based on denoiser priors, their theoretical analyses, statistical interpretations, as well as connections to related approaches such as score matching and diffusion models (Venkatakrishnan et al., ; Chan et al., ; Romano et al., ; Buzzard et al., ; Reehorst & Schniter, ; Sun et al., , ; Ryu et al., ; Xu et al., ; Liu et al., ; Cohen et al., ; Hurault et al., , ; Laumont et al., ; Gan et al., ; Renaud et al., ).

The mathematical relationship between denoising and the score function (the gradient of the log of the image distribution), known as the Tweedie’s formula (Robbins, ; Efron, ) seemingly implies that Gaussian denoising alone might be sufficient for learning priors, independent of the specific characteristics of an inverse problem. However, there is limited research exploring whether broader classes of priors based on pre-trained restoration models could outperform those based on Gaussian denoisers. In this paper, we present evidence that priors derived from deep models pre-trained as general restoration operators can surpass those trained exclusively for Gaussian denoising. We introduce a novel method called *Stochastic deep Restoration Priors (ShaRP)*, which provides a principled approach to integrate an ensemble of general restoration models as priors to regularize inverse problems. ShaRP is conceptually related to several recent papers exploring priors specified using other types of image restoration operators, such as, for example, image super-resolution models (Zhang et al., ; Liu et al., ; Gilton et al., ; Hu et al., ). However, unlike these methods, ShaRP provides a richer and more flexible representation of image priors by considering restoration models trained on a range of degradation types (e.g., various undersampling masks in MRI or blur kernels in image deblurring). By using more versatile restoration models, ShaRP improves upon traditional methods using Gaussian denoiser priors in two key ways: (a) ShaRP improved performance by using restoration models better suited to mitigating the structured artifacts that arise during inference, which are often linked to the characteristics of the underlying inverse problem. (b) Unlike Gaussian denoisers, the restoration models in ShaRP can often be directly trained in a self-supervised manner without fully sampled measurement data.

We present new theoretical and numerical results highlighting the potential of using an ensemble of restoration models as image priors. Our first theoretical result introduces a novel notion of regularization for inverse problems corresponding to the average of likelihoods associated with the degraded observations of an image. The proposed regularizer has an intuitive interpretation as promoting solutions whose multiple degraded observations resemble realistic degraded images. We show that ShaRP seeks to minimize an objective function containing this regularizer. Our second theoretical result analyzes the convergence of ShaRP iterations when using both exact and inexact minimum mean squared error (MMSE) restoration operators. Numerically, we show the practical relevance of ShaRP by applying it to MRI reconstruction with varying undersampling patterns and rates, using a fixed-rate pre-trained MRI reconstruction network as a prior. We also show that ShaRP can use a pre-trained image deblurring model to perform single image super-resolution (SISR). Our numerical experiments show that ShaRP effectively adapts the pre-trained restoration model as a prior, outperforming existing methods based on image denoisers and diffusion models, and achieving state-of-the-art results. Our experiments additionally highlight the benefit of using restoration models as priors by considering a setting where only undersampled and noisy MRI data is available for pre-training the prior. In such cases, self-supervised training of a restoration model is feasible, whereas training a Gaussian denoiser requires fully sampled data.

## 2 Background

Inverse Problems. Many computational imaging tasks can be formulated as inverse problems, where the goal is to reconstruct an unknown image ${\bm{x}}\in\mathbb{R}^{n}$ from its corrupted measurement

$$ ${\bm{y}}={\bm{A}}{\bm{x}}+{\bm{e}},$ (1) $$

where ${\bm{A}}\in\mathbb{R}^{m\times n}$ is a measurement operator and ${\bm{e}}\in\mathbb{R}^{m}$ is the noise. A common approach to addressing inverse problems is to formulate them as an optimization problem

$$ ${\widehat{\bm{x}}}\in\mathop{\mathsf{arg\,min}}_{{\bm{x}}\in\mathbb{R}^{n}}f({ \bm{x}})\quad\text{with}\quad f({\bm{x}})=g({\bm{x}})+h({\bm{x}})\ ,$ (2) $$

where $g$ is the data-fidelity term that quantifies the fit to the measurement ${\bm{y}}$ and $h$ is a regularizer that incorporates prior information on ${\bm{x}}$. For instance, typical functions used in imaging inverse problems are the least-squares term $g({\bm{x}})=\frac{1}{2}\left\lVert{\bm{A}}{\bm{x}}-{\bm{y}}\right\rVert_{2}^{2}$ and the total variation (TV) regularizer $h({\bm{x}})=\tau\left\lVert{\bm{D}}{\bm{x}}\right\rVert_{1}$, where ${\bm{D}}$ is the image gradient and $\tau>0$ is a regularization parameter.

Deep Learning. DL has emerged as a powerful tool for addressing inverse problems (McCann et al., ; Ongie et al., ; Wen et al., ). Instead of explicitly defining a regularizer, DL methods use
deep neural networks (DNNs) to map the measurements to the desired images (Wang et al., ; Jin et al., ; Kang et al., ; Chen et al., ; Delbracio et al., ; Delbracio & Milanfar, ). Model-based DL (MBDL) is a widely-used sub-family of DL algorithms that integrate physical measurement models with priors specified using CNNs (see reviews by Ongie et al. (); Monga et al. ()). The
literature of MBDL is vast, but some well-known examples include plug-and-play priors (PnP), regularization by denoising (RED), deep unfolding (DU), compressed sensing using generative models (CSGM), and deep equilibrium models (DEQ) (Bora et al., ; Romano et al., ; Zhang & Ghanem, ; Hauptmann et al., ; Gilton et al., ; Liu et al., ; Hu et al., ). These approaches come with different trade-offs in terms of imaging performance, computational and memory complexity, flexibility, need for supervision, and theoretical understanding.

Denoisers as Priors. Score-based models (SBMs) are a powerful subset of DL methods for solving inverse problems that use deep Gaussian denoisers as imaging priors. Plug-and-Play (PnP) methods can be viewed as SBMs that incorporate denoisers within iterative optimization algorithms (see recent reviews (Ahmad et al., ; Kamilov et al., )). These approaches construct a cost function by combining an explicit likelihood with a score function implicitly defined by the denoiser prior. Over the past few years, numerous variants of PnP have been developed (Venkatakrishnan et al., ; Romano et al., ; Metzler et al., ; Zhang et al., ; Meinhardt et al., ; Dong et al., ; Zhang et al., ; Wei et al., ; Hurault et al., ), which has motivated an extensive research into their theoretical properties and empirical effectiveness (Chan et al., ; Buzzard et al., ; Ryu et al., ; Sun et al., ; Tirer & Giryes, ; Teodoro et al., ; Xu et al., ; Sun et al., ; Cohen et al., ; Hu et al., ; Laumont et al., ; Hurault et al., ; Gan et al., ; Cohen et al., ; Fang et al., ; Renaud et al., ; Hu et al., ; Renaud et al., ; Terris et al., ).
Diffusion Models (DMs) represent another category of SBMs; they are trained to learn the score function of the underlying probability distribution governed by stochastic differential equations (SDEs) (Ho et al., ; Song et al., ). Once trained, these models can be used as powerful priors for inverse problems by leveraging their learned score functions. Specifically, pre-trained DMs facilitate posterior sampling by guiding the denoising process to generate data consistent with observed measurements. This approach enables DMs to address inverse problems, often achieving impressive perceptual performance even for highly ill-posed inverse problems (Chung et al., ; Zhu et al., ; Wang et al., ; Feng et al., ; Sun et al., ; Wu et al., ; Song et al., ; Hu et al., ; Alçalar & Akçakaya, ; Zhao et al., ; Rout et al., ; Bai et al., ).

Restoration Networks as Priors. In addition to denoiser-based methods, recent work has also considered using restoration models as implicit priors for solving inverse problems (Zhang et al., ; Liu et al., ; Gilton et al., ; Hu et al., ). It has been observed that pre-trained restoration models can be effective priors for addressing unseen inverse problems, sometimes surpassing traditional denoiser-based approaches (Hu et al., ). However, existing methods present two main limitations.
First, restoration models considered so far have relied on a fixed prior tailored to a specific degradation. Although deep restoration models can be trained in various settings—such as different blur kernels for image deblurring or diverse undersampling masks for MRI reconstruction—current approaches do not leverage this capability, limiting their robustness to diverse artifacts. Second, prior work has not explored the potential of learning restoration priors directly from undersampled measurements, without access to fully sampled data. Unlike Gaussian denoisers, training without fully sampled data is natural for restoration models (Yaman et al., ; Liu et al., ; Tachella et al., ; Chen et al., ; Millard & Chiew, ; Gan et al., ). It is also worth highlighting the related work that has explored using corrupt measurements for training Ambient DMs (Daras et al., ; Aali et al., ). Ambient DMs seek to sample from $p_{\bm{x}}$ using DMs trained directly on undersampled measurements. Thus, during inference Ambient DMs assume access to the image prior $p_{\bm{x}}$, while ShaRP only assumes access to the ensemble of likelihoods of multiple degraded observations.

Our contribution. (1) We propose ShaRP, a new framework for solving inverse problems leveraging a set of priors implicit in a pre-trained deep restoration network. ShaRP generalizes Regularization by Denoising (RED) (Romano et al., ) and Stochastic Denoising Regularization (SNORE) (Renaud et al., ) by using more flexible restoration operators and extends Deep Restoration Priors (DRP) (Hu et al., ) by using multiple restoration priors instead of relying on a single one. (2) We introduce a novel regularization concept for inverse problems that encourages solutions that produce degraded versions closely resembling real degraded images. For example, our regularizer favors an MR image solution only if its various degraded versions are consistent with the characteristics of actual degraded MR images. (3) We show that ShaRP minimizes a composite objective that incorporates our proposed regularizer. We provide a theoretical analysis of its convergence for both exact and approximate MMSE restoration operators. (4) We implement ShaRP with both supervised and self-supervised restoration models as priors and test it on two inverse problems: compressed sensing MRI (CS-MRI) and single-image super-resolution (SISR). Our results highlight the capability of restoration models to achieve state-of-the-art performance. Notably, in the MRI context, we show that restoration networks trained directly on subsampled and noisy MRI data can serve as effective priors, a scenario where training traditional Gaussian denoisers is infeasible.

Figure: Figure 1: A restoration network trained on a set of tasks $\{{\mathbf{H}}_{i}\}$ can be used as a prior within ShaRP to address different target inverse problems without the need for retraining.
Refer to caption: x1.png

## 3 Stochastic Deep Restoration Priors

ShaRP is presented in Algorithm 1. It considers a prior based on a deep restoration model ${\mathsf{R}}({\bm{s}},{\mathbf{H}})$ pre-trained using the family of $b$ degradation operators $\{{\mathbf{H}}_{1},\cdots,{\mathbf{H}}_{b}\}$, such as blur kernels or MRI masks. More specifically, the deep restoration model ${\mathsf{R}}$ is trained to solve the following set of restoration problems

$$ ${\bm{s}}_{i}={\mathbf{H}}_{i}{\bm{x}}+{\bm{n}}_{i}\quad\text{with}\quad{\bm{x} }\sim p_{\bm{x}},\quad{\bm{n}}_{i}\sim{\mathcal{N}}(0,\sigma^{2}{\mathbf{I}}), \quad i\in\{1,\cdots,b\},$ (3) $$

where ${\bm{n}}_{i}$ are the AWGN vectors with variance $\sigma^{2}$ and $p_{\bm{x}}$ denotes the probability distribution of the target images. Importantly, the restoration problems (3) are used exclusively for training ${\mathsf{R}}$ and do not need to match the target inverse problem (1), which involves the measurement operator ${\bm{A}}$.

Figure: Algorithm 1 Stochastic deep Restoration Priors (ShaRP)

Our analysis below shows that ${\widehat{\nabla}}f$ corresponds to a stochastic approximation of an objective function of form $f=g+h$. Similar to traditional stochastic gradient methods, ShaRP can be implemented using various selection strategies for the degradation operators. Although the algorithm above is described using a finite set of $b$ operators, our theoretical analysis adopts a more general approach by considering the degradation operator ${\mathbf{H}}$ to be sampled from a distribution ${\mathbf{H}}\sim p_{{\mathbf{H}}}$.

Each iteration of ShaRP has an intuitive interpretation, where the next solution is obtained by combining the gradient of the data-fidelity term $\nabla g$ and the residual of restored image corresponding to the selected degradation operator. When ${\mathbf{H}}_{i}={\mathbf{I}}$ for all $i\in\{1,\cdots,b\}$, then the restoration prior reduces to the Gaussian denoiser, and ShaRP can be viewed as an instance of the Regularization by Denoising (RED) method (Romano et al., ) and Stochastic Denoising Regularization (SNORE) (Renaud et al., ). On the other hand, when $b=1$, then ShaRP can be viewed as the instance of the Deep Restoration Priors (DRP) method (Hu et al., ). Thus, ShaRP can be viewed as a more versatile generalization of both that can account for various degradation operators.

## 4 Theoretical analysis of ShaRP

We present two theoretical results on ShaRP. The first establishes the regularizer minimized by ShaRP, while the second analyzes its convergence with inexact MMSE operators.

Consider a restoration model that perform MMSE estimation of ${\bm{x}}\in\mathbb{R}^{n}$ for problems (3)

$$ ${\mathsf{R}}^{\ast}({\bm{s}},{\mathbf{H}})=\mathbb{E}\left[{\bm{x}}|{\bm{s}},{ \mathbf{H}}\right]=\int{\bm{x}}\,p({\bm{x}}|{\bm{s}},{\mathbf{H}}){\mathsf{\,d }}{\bm{x}}=\frac{1}{p({\bm{s}}|{\mathbf{H}})}\int{\bm{x}}\,G_{\sigma}({\bm{s}} -{\mathbf{H}}{\bm{x}}){p_{{\bm{x}}}}({\bm{x}}){\mathsf{\,d}}{\bm{x}}.$ (4) $$

where we used the probability density of the observation ${\bm{s}}$ conditioned on the operator ${\mathbf{H}}$

$$ $p({\bm{s}}|{\mathbf{H}})=\int G_{\sigma}({\bm{s}}-{\mathbf{H}}{\bm{x}})p_{\bm{ x}}({\bm{x}}){\mathsf{\,d}}{\bm{x}}.$ (5) $$

The function $G_{\sigma}$ in (5) denotes the Gaussian density function with the standard deviation $\sigma>0$.

We define the ShaRP regularizer as

$$ $h({\bm{x}})=\tau\mathbb{E}_{{\bm{s}}\sim G_{\sigma}({\bm{s}}-{\mathbf{H}}{\bm{ x}}),{\mathbf{H}}\sim p_{\mathbf{H}}}\left[-{\mathsf{log}}\,p({\bm{s}}|{ \mathbf{H}})\right],$ (6) $$

where $\tau>0$ is the regularization parameter and $p_{{\mathbf{H}}}$ is the distribution of all considered degradation operators. The regularizer $h$ is minimized when the degraded versions of ${\bm{x}}$ have a high likelihood under the probability density of degraded observations, $p({\bm{s}}|{\mathbf{H}})$. In other words, $h$ identifies an image as a valid solution if its degraded observations resemble realistic degraded images.

We are now ready to state our first theoretical result.

###### Theorem 1 .

Assume that the prior density $p_{\bm{x}}$ is non-degenerate over $\mathbb{R}^{n}$ and let ${\mathsf{R}}^{\ast}$ be the MMSE restoration operator (4) corresponding to the restoration problems (3).
Then, we have that

$$ $\nabla h({\bm{x}})=\frac{\tau}{\sigma^{2}}\left(\mathbb{E}_{{\bm{s}}\sim G_{ \sigma}({\bm{s}}-{\mathbf{H}}{\bm{x}}),{\mathbf{H}}\sim p_{\mathbf{H}}}\left[{ \mathbf{H}}^{\mathsf{T}}{\mathbf{H}}({\bm{x}}-{\mathsf{R}}^{\ast}({\bm{s}},{ \mathbf{H}}))\right]\right),$ (7) $$

where $h$ is the ShaRP regularizer in (6).

The proof is provided in the appendix. Note that the expression within the square parenthesis in (7) matches the ShaRP update in Line 4 of Algorithm 1, which directly implies that ShaRP using the exact MMSE restoration operator ${\mathsf{R}}^{\ast}$ is a stochastic gradient method for minimizing $f=g+h$, where $g$ is the data-fidelity term and $h$ is the ShaRP regularizer in (6).

We now present the convergence analysis of ShaRP under a restoration operator ${\mathsf{R}}$ that *approximates* the true MMSE restoration operator ${\mathsf{R}}^{\ast}$. For a given degraded observation ${\bm{s}}={\mathbf{H}}{\bm{x}}+{\bm{n}}$ with ${\mathbf{H}}\sim p_{{\mathbf{H}}}$ and ${\bm{n}}\sim{\mathcal{N}}(0,\sigma^{2}{\mathbf{I}})$, we define the stochastic gradient used by ShaRP

$$ ${\widehat{\nabla}}f({\bm{x}})=\nabla g({\bm{x}})+{\widehat{\nabla}}h({\bm{x}}) \quad\text{with}\quad{\widehat{\nabla}}h({\bm{x}})\,\coloneqq\,\frac{\tau}{ \sigma^{2}}{\mathbf{H}}^{\mathsf{T}}{\mathbf{H}}({\bm{x}}-{\mathsf{R}}({\bm{s} },{\mathbf{H}})).$ (8) $$

Since ${\mathsf{R}}$ is an inexact MMSE restoration operator, we also define the bias vector

$$ ${\bm{b}}({\bm{x}})=\frac{\tau}{\sigma^{2}}\mathbb{E}_{{\bm{s}}\sim G_{\sigma}( {\bm{s}}-{\mathbf{H}}{\bm{x}}),{\mathbf{H}}\sim p_{\mathbf{H}}}\left[{\mathbf{ H}}^{\mathsf{T}}{\mathbf{H}}({\mathsf{R}}^{\ast}({\bm{s}},{\mathbf{H}})-{ \mathsf{R}}({\bm{s}},{\mathbf{H}}))\right],$ (9) $$

which quantifies the average difference between the exact and inexact MMSE restoration operators.
Our analysis requires three assumptions that jointly serve as sufficient conditions for our theorem.

###### Assumption 1 .

The function $f$ has a finite minimum $f^{\ast}>-\infty$ and the gradient $\nabla f$ is Lipschitz continuous with constant $L>0$.

This is a standard assumption used in the analysis of gradient-based algorithms (see (Nesterov, ), for example). It is satisfied by a large number of functions, including the traditional least-squares data-fidelity function.

###### Assumption 2 .

The stochastic gradient has a bounded variance for all ${\bm{x}}\in\mathbb{R}^{n}$, which means that there exists a constant $\nu>0$ such that

$$ $\mathbb{E}\left[\left\|{\widehat{\nabla}}f({\bm{x}})-\mathbb{E}\left[{\widehat {\nabla}}f({\bm{x}})\right]\right\|_{2}^{2}\right]\leq\nu^{2},$ $$

where expectations are with respect to ${\mathbf{H}}\sim p_{{\mathbf{H}}}$ and ${\bm{s}}\sim G_{\sigma}({\bm{s}}-{\mathbf{H}}{\bm{x}})$.

This is another standard assumption extensively used in the analysis of online or stochastic optimization algorithms (Bertsekas, ; Ghadimi & Lan, ).

###### Assumption 3 .

The bias is bounded, which means that there exists $\varepsilon>0$ such that for all ${\bm{x}}\in\mathbb{R}^{n}$

$$ $\|{\bm{b}}({\bm{x}})\|_{2}\leq\varepsilon.$ $$

Note that our only assumption on the bias is that it is bounded, which is a relatively mild assumption in the context of biased stochastic gradient methods (Demidovich et al., ).

###### Theorem 2 .

Run ShaRP for $t\geq 1$ iterations using the step-size $0<\gamma\leq 1/L$ under Assumptions 1-3. Then, the sequence ${\bm{x}}^{k}$ generated by ShaRP satisfies

$$ $\mathbb{E}\left[\frac{1}{t}\sum_{k=1}^{t}\|\nabla f({\bm{x}}^{k-1})\|_{2}^{2} \right]\leq\frac{2}{\gamma t}(f({\bm{x}}^{0})-f^{\ast})+\gamma L\nu^{2}+ \varepsilon^{2}.$ $$

The proof is provided in the appendix. This theorem shows that *in expectation*, ShaRP minimizes the norm of the gradient $\nabla f$ up to an error term that has two components, $\gamma L\nu^{2}$ and $\epsilon^{2}$. Since the first component depends on $\gamma$, it can be made as small as desired by controlling the step-size $\gamma$. The second component only depends on the magnitude of the bias $\varepsilon$, which, in turn, directly depends on the accuracy of the restoration operator relative to the true MMSE restoration operator ${\mathsf{R}}^{\ast}$.

###### Theorem 1 .

###### Assumption 1 .

###### Assumption 2 .

###### Assumption 3 .

###### Theorem 2 .

## 5 Numerical Results

We numerically validate ShaRP on two inverse problems of the form ${\bm{y}}={\bm{A}}{\bm{x}}+{\bm{e}}$: (*Compressive Sensing MRI (CS-MRI)* and (b) *Single Image Super Resolution (SISR)*. In both cases, ${\bm{e}}$ represents additive white Gaussian noise (AWGN). For the data-fidelity term in eq. (2), we use the $\ell_{2}$-norm loss for both problems. Quantitative performance is evaluated using Peak Signal-to-Noise Ratio (PSNR) and Structural Similarity Index (SSIM). Additionally, for the SISR task, we include the Learned Perceptual Image Patch Similarity (LPIPS) metric to evaluate perceptual quality. Additional numerical results are provided in the supplementary material.

Figure: Figure 2: Convergence of ShaRP for $4\times$ accelerated MRI reconstruction on the fastMRI dataset. (a)-(b) depict the convergence behavior of ShaRP using restoration operators trained in a supervised manner, while (c)-(d) correspond to those trained in a self-supervised manner. The plots illustrate the average distance $\|{\bm{x}}^{k}-{\bm{x}}^{k-1}\|^{2}_{2}$ and PSNR relative to the ground truth, as a function of the iteration number, with shaded regions representing the standard deviation. Note the stable convergence of ShaRP with both types of priors.
Refer to caption: x2.png

### 5.1 CS-MRI setting

The measurement of CS-PMRI can be modeled as ${\bm{y}}={\bm{P}}{\bm{F}}{\bm{S}}{\bm{x}}+{\bm{e}}$, where ${\bm{P}}$ is the k-space subsampling pattern, ${\bm{F}}$ is the Fourier transform operator, ${\bm{S}}=({\bm{S}}_{1},\cdots,{\bm{S}}_{n_{c}})$ are the multi-coil sensitivity maps, and ${\bm{e}}$ is the noise vector.

Dataset. We simulated multi-coil subsampled measurements using T2-weighted human brain MRI data from the open-access fastMRI dataset, which comprises 4,912 fully sampled multi-coil slices for training and 470 slices for testing. Each slice has been cropped into a complex-valued image with dimensions $320\times 320$. The coil sensitivity maps for each slice are precomputed using the ESPIRiT algorithm (Uecker et al., ). We simulated a Cartesian sampling pattern that subsamples along the $k_{y}$ dimension while fully sampling along the $k_{x}$ dimension.

Ensemble of Restoration Priors for CS-MRI. Recent methods, such as InDI (Delbracio & Milanfar, ) and I^2SB (Liu et al., ), introduce controllable processes for training an ensemble of restoration priors, where each prior functions as an MMSE restoration operator tailored to a specific setting. Building on this approach, we trained an $8\times$ uniform subsampling CS-MRI model with 8 distinct masks as our restoration prior. Similar to InDI, we decompose the original MRI degradation operator ${\bm{M}}$ into several convex combinations of the original task ${\bm{M}}$ and the identity mapping ${\mathbf{I}}$, represented by the new degradation operator ${\mathbf{H}}_{\alpha}=(1-\alpha){\mathbf{I}}+\alpha{\bm{M}}$, with $\alpha$ controlling the degradation level. By selecting a range of $\alpha$ values, we create an ensemble of restoration tasks. Training the restoration network ${\mathsf{R}}$ to handle all these tasks allows it to function as an ensemble of MMSE restoration operators, ${\mathsf{R}}({\bm{s}},{\mathbf{H}}_{\alpha})=\mathbb{E}\left[{\bm{x}}|{\bm{s}}
,{\mathbf{H}}_{\alpha}\right]$. We used the MSE loss to train the restoration model.

Training restoration priors without groundtruth.
When fully-sampled ground truth images are not available for training restoration priors, MRI restoration priors can be trained in self-supervised manner (Yaman et al., ; Millard & Chiew, ; Gan et al., ; Hu et al., ). In self-supervised training, rather than using the ground-truth image as the label, a separate subsampled measurement serves as the label. In such cases, we can train our priors using a weighted $\ell_{2}$
loss function, following the self-supervised approach in (Gan et al., ). We thus train the $8\times$ uniform subsampling CS-MRI model to handle eight distinct restoration operators, each corresponding to a different sampling mask.

Additional details on supervised and self-supervised restoration network training and our CS-MRI sampling masks are in Section B.1 of the appendix.

With the pre-trained $8\times$ models as ensembles of restoration priors, we evaluate ShaRP’s performance across a variety of configurations, including two sub-sampling rates (4$\times$ and 6$\times$), two mask types (uniform and random), and three noise levels ($\sigma=0.005$, $0.01$, and $0.015$).

**Table 1: Quantitative comparison of ShaRP with several baselines for CS-MRI using uniform masks at undersampling rates of 4 and 6 on fastMRI dataset. The best and second best results are highlighted. Notably, ShaRP outperforms SOTA methods based on denoisers and diffusion models.**
|  | 4$\times$ Uniform | 6$\times$ Uniform |  |  |  |  |  |  |  |  |  |  |
| --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- |
| Noise level | $\sigma$ = 0.005 | $\sigma$ = 0.010 | $\sigma$ = 0.015 | $\sigma$ = 0.005 | $\sigma$ = 0.010 | $\sigma$ = 0.015 |  |  |  |  |  |  |
| Metrics | PSNR | SSIM | PSNR | SSIM | PSNR | SSIM | PSNR | SSIM | PSNR | SSIM | PSNR | SSIM |
| Zero-filled | 26.93 | 0.848 | 26.92 | 0.847 | 26.90 | 0.848 | 22.62 | 0.728 | 22.60 | 0.726 | 22.59 | 0.721 |
| TV | 31.17 | 0.923 | 31.08 | 0.921 | 30.91 | 0.915 | 25.00 | 0.806 | 24.94 | 0.803 | 24.33 | 0.755 |
| PnP-FISTA | 35.88 | 0.938 | 31.14 | 0.894 | 30.32 | 0.846 | 26.30 | 0.822 | 25.97 | 0.786 | 25.46 | 0.747 |
| PnP-ADMM | 35.76 | 0.941 | 32.36 | 0.878 | 30.66 | 0.838 | 26.13 | 0.808 | 25.90 | 0.776 | 25.51 | 0.742 |
| DRP | 35.52 | 0.936 | 32.32 | 0.914 | 30.57 | 0.901 | 29.51 | 0.872 | 28.52 | 0.882 | 28.35 | 0.876 |
| DPS | 32.62 | 0.888 | 31.39 | 0.870 | 30.29 | 0.856 | 30.53 | 0.862 | 29.41 | 0.843 | 28.63 | 0.830 |
| DDS | 35.21 | 0.937 | 35.03 | 0.935 | 34.51 | 0.925 | 31.02 | 0.889 | 30.84 | 0.888 | 30.79 | 0.888 |
| ShaRP | 37.59 | 0.963 | 35.81 | 0.951 | 34.92 | 0.942 | 33.42 | 0.940 | 32.86 | 0.932 | 32.09 | 0.922 |

Baselines. ShaRP was compared against several baseline methods, including denoiser-based approaches (PnP-FISTA (Kamilov et al., ), PnP-ADMM (Chan et al., )) and diffusion model-based methods (DPS (Chung et al., ), DDS (Chung et al., )). To highlight the advantages of using a stochastic set of restoration operators, we also compared ShaRP with the DRP method (Hu et al., ), which applies only a single restoration operator. Additional details related to the baseline methods can be found in Section B.1 of the appendix.

Results with supervised MMSE Restoration operator. Figure 2 illustrates the convergence behavior of ShaRP on the test set with an acceleration factor of $R=6$ and additional noise $\sigma=0.01$. Table 1 provides a quantitative comparison of reconstruction performance across different acceleration factors and noise levels using a uniform sub-sampling mask. In all configurations, ShaRP consistently outperforms the baseline methods. The use of a set of restoration operators clearly enhances ShaRP’s performance, highlighting the effectiveness of employing multiple operators to maximize the regularization information provided by the restoration model. Figure 3 presents visual reconstructions for two test scenarios, where ShaRP accurately recovers fine brain details, particularly in the zoomed-in regions, while baseline methods tend to oversmooth or introduce artifacts. These results highlight ShaRP’s superior ability to manage structured artifacts and preserve fine details, outperforming both denoiser-based and diffusion model-based methods.

Figure: Figure 3: Visual comparison of ShaRP with baseline methods on CS-MRI. The top row shows results for a $4\times$ random mask with noise $\sigma=0.005$, and the bottom row for a $6\times$ random mask with noise $\sigma=0.015$. PSNR and SSIM values are in the top-left corner of each image. Error maps and zoomed-in areas highlight differences. Notably, ShaRP with stochastic priors outperforms state-of-the-art methods using denoiser and diffusion model priors.
Refer to caption: x3.png

Results with self-supervised MMSE Restoration operator.
We further evaluate ShaRP’s performance using an restoration model, learned in a self-supervised manner, as introduced in (Gan et al., ). In this setting, we compare ShaRP against two classical methods for CS-MRI reconstruction without groundtruth: TV (Block et al., ) and GRAPPA (Griswold et al., ) and a recent state-of-the-art self-supervised deep unrolling method: SPICER (Hu et al., ). As shown in Table 2, even in scenarios where only incomplete measurements ($8\times$ subsampled measurement) are available, ShaRP can effectively apply a self-supervised trained restoration prior to various reconstruction tasks. ShaRP using self-supervised restoration prior even outperforms DPIR and DPS that use Gaussian denoisers trained using fully-sampled ground truth images (see Table 4 in the appendix). Note that given only undersampled measurements, training Gaussian denoisers is not feasible.

**Table 2: PSNR (dB) and SSIM values for ShaRP with a self-supervised pre-trained restoration operator, compared to several baselines for CS-MRI with random undersampling masks at rates of 4 and 6 on the fastMRI dataset. The best and second best results are highlighted. For reference, the highlighted row presents a PnP method using a Gaussian denoiser, which requires fully sampled data for training. Note the excellent performance of ShaRP even using priors trained without fully-sampled ground-truth data.**
|  | 4$\times$ Random | 6$\times$ Random |  |  |  |  |  |  |  |  |  |  |
| --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- |
| Noise level | $\sigma$ = 0.005 | $\sigma$ = 0.010 | $\sigma$ = 0.015 | $\sigma$ = 0.005 | $\sigma$ = 0.010 | $\sigma$ = 0.015 |  |  |  |  |  |  |
| Metrics | PSNR | SSIM | PSNR | SSIM | PSNR | SSIM | PSNR | SSIM | PSNR | SSIM | PSNR | SSIM |
| PnP-ADMM | 28.83 | 0.842 | 28.39 | 0.816 | 27.70 | 0.786 | 25.59 | 0.776 | 25.19 | 0.740 | 24.93 | 0.728 |
| ADMM-TV | 28.14 | 0.866 | 28.06 | 0.863 | 27.96 | 0.859 | 24.55 | 0.782 | 24.33 | 0.750 | 24.28 | 0.736 |
| GRAPPA | 28.09 | 0.792 | 25.39 | 0.699 | 23.94 | 0.649 | 25.67 | 0.737 | 23.72 | 0.646 | 22.51 | 0.595 |
| SPICER | 31.87 | 0.901 | 31.67 | 0.889 | 31.50 | 0.887 | 30.18 | 0.871 | 30.05 | 0.863 | 30.01 | 0.860 |
| $\text{ShaRP}^{\textbf{self}}$ | 33.87 | 0.909 | 33.64 | 0.900 | 33.21 | 0.892 | 30.87 | 0.899 | 30.36 | 0.890 | 30.21 | 0.875 |

Figure: Figure 4: Visual comparison of ShaRP with several well-known methods on $2\times$ SISR. The top row shows results for SISR with gaussian blur kernel with $\sigma=1.25$, while the bottom row shows results for SISR with gaussian blur kernel with $\sigma=1.5$. The quantities in the top-left corner of each image provide PSNR and SSIM values for each method. The squares at the bottom of each image visualize the zoomed area in the image.
Refer to caption: x4.png

### 5.2 Single Image Super Resolution (SISR)

The measurement operator in SISR can be written as ${\bm{A}}={\bm{S}}{\bm{K}}$, where ${\bm{K}}$ represents convolution with the blur kernel, and ${\bm{S}}$ performs standard $d$-fold down-sampling. In our experiments, we use two Gaussian blur kernels ${\bm{k}}$ , each with distinct standard deviations (1.25 and 1.5), and with down-sampling factor of 2. Both noisy and noise-free cases are considered to evaluate the noisy robustness of ShaRP. We randomly selected 100 images from the ImageNet test set, as provided in DiffPIR(^1^11[https://github.com/yuanzhi-zhu/DiffPIR/tree/main/testsets](https://github.com/yuanzhi-zhu/DiffPIR/tree/main/testsets)).

Ensemble of Restoration Priors for Image Deblurring. Following the approach used to train our CS-MRI restoration prior, we decompose the original deblurring task represented by the Gaussian blur operator ${\bm{K}}$ into convex combinations of the original task and the identity mapping ${\mathsf{I}}$. This results in a new degradation operator defined as ${\mathbf{H}}_{\alpha}=(1-\alpha){\mathsf{I}}+\alpha{\bm{K}}$, where $\alpha$ controls the degradation level. By varying $\alpha$, we generate multiple degradation operators, allowing us to train the restoration network ${\mathsf{R}}$ to handle all these operators. This training enables ${\mathsf{R}}$ to function as an ensemble of MMSE restoration operators, expressed as ${\mathsf{R}}({\bm{s}},{\mathbf{H}}_{\alpha})=\mathbb{E}[{\bm{x}}\mid{\bm{s}},{
\mathbf{H}}_{\alpha}]$, where ${\bm{s}}$ is the degraded image and ${\bm{x}}$ is the original image. The original deblurring operator ${\bm{K}}$ corresponds to convolution with a Gaussian blur kernel of size $31\times 31$ and standard deviation 3. More details on the deblurring restoration network training are in Section B.2 of the Appendix.

Baselines. We compared ShaRP with several baseline methods, including DPIR (Zhang et al., ), DPS (Chung et al., ), DDNM (Wang et al., ), and DiffPIR (Zhu et al., ). DPIR represents the state-of-the-art (SOTA) PnP method that uses pre-trained denoisers as priors to address SISR. In contrast, DPS, DDNM, and DiffPIR use different sampling strategies to leverage pre-trained unconditional diffusion models for solving SISR. More baseline details can be found in Section B.2 of the Appendix.

Results on SISR with deblurring prior. Figure 4 shows the visual reconstruction results for two settings with different blur kernels. As demonstrated, ShaRP successfully recovers most features and maintains high data consistency with the available measurements. Table 3 provides a quantitative comparison of ShaRP against other baseline methods, evaluated across various blur kernels and noise levels. ShaRP achieves the highest PSNR and SSIM values but ranks second in perceptual performance (LPIPS). This is consistent with the SOTA perceptual performance of DMs on SISR. However, note how the use of a deblurring prior within ShaRP enables it to recover fine details, ensuring overall competitiveness of the perceptual quality of the ShaRP solutions.

**Table 3: Quantitative comparison of ShaRP with several baselines for SISR based on two different blur kernels on ImageNet dataset. The best and second best results are highlighted. Notably, ShaRP outperforms SOTA methods based on denoisers and diffusion models.**
|  |  |  |  |  |  |  |  |  |  |  |  |  |
| --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- |
| Noise level | Noiseless | $\sigma$ = 0.01 | Noiseless | $\sigma$ = 0.01 |  |  |  |  |  |  |  |  |
| Metrics | PSNR | SSIM | LPIPS | PSNR | SSIM | LPIPS | PSNR | SSIM | LPIPS | PSNR | SSIM | LPIPS |
| DPIR | 28.10 | 0.809 | 0.305 | 28.05 | 0.807 | 0.308 | 27.90 | 0.803 | 0.314 | 27.87 | 0.800 | 0.314 |
| DDNM | 27.53 | 0.786 | 0.240 | 27.49 | 0.784 | 0.246 | 27.02 | 0.764 | 0.264 | 27.01 | 0.763 | 0.267 |
| DPS | 24.68 | 0.661 | 0.395 | 24.60 | 0.657 | 0.399 | 24.50 | 0.657 | 0.403 | 24.44 | 0.655 | 0.406 |
| DiffPIR | 28.92 | 0.852 | 0.152 | 28.63 | 0.839 | 0.169 | 28.59 | 0.834 | 0.172 | 28.02 | 0.819 | 0.185 |
| DRP | 29.28 | 0.868 | 0.207 | 28.87 | 0.848 | 0.248 | 28.24 | 0.836 | 0.235 | 28.01 | 0.822 | 0.278 |
| ShaRP | 30.09 | 0.891 | 0.179 | 29.03 | 0.852 | 0.223 | 29.28 | 0.872 | 0.209 | 28.06 | 0.821 | 0.268 |

## 6 Conclusion

The work presented in this paper proposes a new ShaRP method for solving imaging inverse problems by using pre-trained restoration network as a prior, presents its theoretical analysis, and applies the method to two well-known inverse problems. Unlike previous approaches that relied on Gaussian denoisers or a single restoration prior, our method uses a set of restoration priors, each corresponding to different degradation settings. The numerical validation shows that ShaRP benefits from stochastically using multiple degradation priors, leading to better results. A key conclusion is the potential effectiveness of exploring priors beyond those defined by traditional Gaussian denoisers and a specific restoration operator.

## Ethics Statement

To the best of our knowledge this work does not give rise to any significant ethical concerns.

## Acknowledgments

The authors extend their sincere gratitude to Regev Cohen and Ashok Popat for their valuable feedback. This research was supported by NSF CCF-2043134. We also thank Google for their generous contribution.

## Appendix A Theoretical Analysis of ShaRP

### A.1 Proof of Theorem 1

###### Theorem .

Assume that the prior density $p_{\bm{x}}$ is non-degenerate over $\mathbb{R}^{n}$ and let ${\mathsf{R}}^{\ast}$ be the MMSE restoration operator (4) corresponding to the restoration problems (3).
Then, we have that

$$ $\nabla h({\bm{x}})=\frac{\tau}{\sigma^{2}}\left(\mathbb{E}_{{\bm{s}}\sim G_{ \sigma}({\bm{s}}-{\mathbf{H}}{\bm{x}}),{\mathbf{H}}\sim p_{\mathbf{H}}}\left[{ \mathbf{H}}^{\mathsf{T}}{\mathbf{H}}({\bm{x}}-{\mathsf{R}}^{\ast}({\bm{s}},{ \mathbf{H}}))\right]\right),$ $$

where $h$ is the ShaRP regularizer in (6).

###### Proof.

The ShaRP regularizer $h({\bm{x}})$ is defined as

$$ $\displaystyle h({\bm{x}})$ $\displaystyle=\tau\mathbb{E}_{{\bm{s}}\sim G_{\sigma}({\bm{s}}-{\mathbf{H}}{ \bm{x}}),{\mathbf{H}}\sim p_{\mathbf{H}}}\left[-{\mathsf{log}}p({\bm{s}}|{ \mathbf{H}})\right]$ $\displaystyle=-\tau\int p({\mathbf{H}})\left[\int G_{\sigma}({\bm{s}}-{\mathbf {H}}{\bm{x}}){\mathsf{log}}p({\bm{s}}|{\mathbf{H}}){\mathsf{\,d}}{\bm{s}} \right]{\mathsf{\,d}}{\mathbf{H}},$ (10) $$

where $G_{\sigma}$ is the Gaussian probability density with variance $\sigma^{2}$ and $p({\bm{s}}|{\mathbf{H}})$ is the likelihood function for the degraded observation given the operator ${\mathbf{H}}$. The expectation over $p({\mathbf{H}})$ accounts for the randomness of the restoration operator ${\mathbf{H}}$.

We start by relating the MMSE restoration operator to the score of the degraded observation

$$ $\nabla p({\bm{s}}|{\mathbf{H}})=\frac{1}{\sigma^{2}}\int\left({\mathbf{H}}{\bm {x}}-{\bm{s}}\right)G_{\sigma}({\bm{s}}-{\mathbf{H}}{\bm{x}})p_{\bm{x}}({\bm{x }}){\mathsf{\,d}}{\bm{x}},$ $$

where $p_{\bm{x}}$ is the prior. By using the definition of the MMSE estimator, we obtain the relationship

$$ $\nabla{\mathsf{log}}p({\bm{s}}|{\mathbf{H}})=\frac{1}{\sigma^{2}}\left({ \mathbf{H}}{\mathsf{R}}^{\ast}({\bm{s}},{\mathbf{H}})-{\bm{s}}\right).$ (11) $$

Consider the function inside the parenthesis in the expression for the ShaRP regularizer (10)

$$ $\rho({\bm{z}})\,\coloneqq\,(G_{\sigma}\ast{\mathsf{log}}p_{{\bm{s}}|{\mathbf{H }}})({\bm{z}})=\int G_{\sigma}({\bm{z}}-{\bm{s}})\,{\mathsf{log}}p({\bm{s}}|{ \mathbf{H}}){\mathsf{\,d}}{\bm{s}},$ $$

where ${\bm{z}}$ has the same dimensions as ${\bm{s}}$ and $\ast$ denotes convolution. The gradient of $\rho$ is given by

$$ $\displaystyle\nabla\rho({\bm{z}})$ $\displaystyle=(\nabla G_{\sigma}*{\mathsf{log}}p_{{\bm{s}}|{\mathbf{H}}})({\bm {z}})=(G_{\sigma}*\nabla{\mathsf{log}}p_{{\bm{s}}|{\mathbf{H}}})({\bm{z}})$ $\displaystyle=\frac{1}{\sigma^{2}}\int G_{\sigma}({\bm{z}}-{\bm{s}})\left[{ \mathbf{H}}{\mathsf{R}}^{\ast}({\bm{s}},{\mathbf{H}})-{\bm{s}}\right]{\mathsf{ \,d}}{\bm{s}}$ $\displaystyle=\frac{1}{\sigma^{2}}\left({\mathbf{H}}\int{\mathsf{R}}^{\ast}({ \bm{s}},{\mathbf{H}})G_{\sigma}({\bm{z}}-{\bm{s}}){\mathsf{\,d}}{\bm{s}}-{\bm{ z}}\right)$ $$

where we used (11). By using ${\bm{z}}={\mathbf{H}}{\bm{x}}$, we write the gradient with respect to ${\bm{x}}$

$$ $\nabla_{\bm{x}}\rho({\mathbf{H}}{\bm{x}})=\frac{1}{\sigma^{2}}{\mathbf{H}}^{ \mathsf{T}}{\mathbf{H}}\left(\int{\mathsf{R}}^{\ast}({\bm{s}},{\mathbf{H}})G_{ \sigma}({\bm{s}}-{\mathbf{H}}{\bm{x}}){\mathsf{\,d}}{\bm{s}}-{\bm{x}}\right)$ $$

By using this expression in (10), we obtain the desired result

$$ $\displaystyle\nabla h({\bm{x}})$ $\displaystyle=-\frac{\tau}{\sigma^{2}}\left[\int p({\mathbf{H}})\int G_{\sigma }({\bm{s}}-{\mathbf{H}}{\bm{x}})\left({\mathbf{H}}^{\mathsf{T}}{\mathbf{H}}({ \mathsf{R}}^{\ast}({\bm{s}},{\mathbf{H}})-{\bm{x}})\right){\mathsf{\,d}}{\bm{s }}{\mathsf{\,d}}{\mathbf{H}}\right]$ $\displaystyle=\frac{\tau}{\sigma^{2}}\mathbb{E}_{{\bm{s}}\sim G_{\sigma}({\bm{ s}}-{\mathbf{H}}{\bm{x}}),{\mathbf{H}}\sim p_{\mathbf{H}}}\left[{\mathbf{H}}^{ \mathsf{T}}{\mathbf{H}}({\bm{x}}-{\mathsf{R}}^{\ast}({\bm{s}},{\mathbf{H}})) \right].$ $$

∎

###### Theorem .

###### Proof.

### A.2 Proof of Theorem 2

###### Theorem .

Run ShaRP for $t\geq 1$ iterations using the step-size $0<\gamma\leq 1/L$ under Assumptions 1-3. Then, the sequence ${\bm{x}}^{k}$ generated by ShaRP satisfies

$$ $\mathbb{E}\left[\frac{1}{t}\sum_{k=1}^{t}\|\nabla f({\bm{x}}^{k-1})\|_{2}^{2} \right]\leq\frac{2}{t}(f({\bm{x}}^{0})-f^{\ast})+\gamma L\nu^{2}+\varepsilon^{ 2}.$ $$

###### Proof.

First note that from the definition of the bias in eq. (9), we have that

$$ $\mathbb{E}\left[{\widehat{\nabla}}f({\bm{x}}^{k-1})\,|\,{\bm{x}}^{k-1}\right]= \nabla f({\bm{x}}^{k-1})+{\bm{b}}({\bm{x}}^{k-1}),$ (12) $$

where the expectation is with respect to ${\bm{s}}\sim G_{\sigma}({\bm{s}}-{\mathbf{H}}{\bm{x}}^{k-1})$ and ${\mathbf{H}}\sim p_{{\mathbf{H}}}$. In order to simplify the notation, we will drop these subscripts from the expectations in the analysis below.

Consider the iteration $k\geq 1$ of ShaRP with inexact MMSE operator

$$ $\displaystyle f({\bm{x}}^{k})$ $\displaystyle\leq f({\bm{x}}^{k-1})+\nabla f({\bm{x}}^{k-1})^{\mathsf{T}}({\bm {x}}^{k}-{\bm{x}}^{k-1})+\frac{L}{2}\|{\bm{x}}^{k}-{\bm{x}}^{k-1}\|_{2}^{2}$ $\displaystyle=f({\bm{x}}^{k-1})-\gamma\nabla f({\bm{x}}^{k-1})^{\mathsf{T}}{ \widehat{\nabla}}f({\bm{x}}^{k-1})+\frac{L\gamma^{2}}{2}\|{\widehat{\nabla}}f( {\bm{x}}^{k-1})\|^{2},$ $$

where in the first line we used the Lipschitz continuity of $\nabla f$. By taking the expectation with respect to ${\bm{s}}\sim G_{\sigma}({\bm{s}}-{\mathbf{H}}{\bm{x}}^{k-1})$ and ${\mathbf{H}}\sim p_{{\mathbf{H}}}$ on both sides of this expression, we get

$$ $\displaystyle\mathbb{E}[f({\bm{x}}^{k})|{\bm{x}}^{k-1}]$ $\displaystyle\leq f({\bm{x}}^{k-1})-\gamma\nabla f({\bm{x}}^{k-1})^{\mathsf{T} }(\nabla f({\bm{x}}^{k-1})+{\bm{b}}({\bm{x}}^{k-1}))+\frac{L\gamma^{2}}{2} \mathbb{E}\left[\|{\widehat{\nabla}}f({\bm{x}}^{k-1})\|_{2}^{2}|{\bm{x}}^{k-1}\right]$ $\displaystyle\leq f({\bm{x}}^{k-1})-\frac{\gamma}{2}\|\nabla f({\bm{x}}^{k-1}) \|_{2}^{2}+\frac{\gamma}{2}\|{\bm{b}}({\bm{x}}^{k-1})\|_{2}^{2}$ $\displaystyle\quad\quad+\frac{L\gamma^{2}}{2}\left(\mathbb{E}\left[\|{\widehat {\nabla}}f({\bm{x}}^{k-1})\|_{2}^{2}|{\bm{x}}^{k-1}\right]-\left(\mathbb{E}[{ \widehat{\nabla}}f({\bm{x}}^{k-1})|{\bm{x}}^{k-1}]\right)^{2}\right)$ $\displaystyle\leq f({\bm{x}}^{k-1})-\frac{\gamma}{2}\|\nabla f({\bm{x}}^{k-1}) \|_{2}^{2}+\frac{\gamma\varepsilon^{2}}{2}+\frac{L\gamma^{2}\nu^{2}}{2}.$ $$

In the second row, we completed the square, applied eq. (12), and used the assumption that $\gamma\leq 1/L$. In the third row, we used the variance and bias bounds in Assumptions 2 and 3. By rearranging the expression, we get the following bound

$$ $\|\nabla f({\bm{x}}^{k-1})\|_{2}^{2}\leq\frac{2}{\gamma}\left(f({\bm{x}}^{k-1} )-\mathbb{E}[f({\bm{x}}^{k})|{\bm{x}}^{k-1}]\right)+L\gamma\nu^{2}+\varepsilon ^{2}$ $$

By taking the total expectation, averaging over $t$ iterations, and using the lower bound $f^{\ast}$, we get the desired result

$$ $\mathbb{E}\left[\frac{1}{t}\sum_{k=1}^{t}\|\nabla f({\bm{x}}^{k-1})\|_{2}^{2} \right]\leq\frac{2}{\gamma t}(f({\bm{x}}^{0})-f^{\ast})+L\gamma\nu^{2}+ \varepsilon^{2}.$ $$

∎

###### Theorem .

###### Proof.

## Appendix B Experiment Details

### B.1 Implementation details of CS-MRI tasks

Subsampling pattern for CS-MRI.
In this paper, we explored two types of subsampling patterns for MRI reconstruction tasks. All undersampling masks were generated by first including a set number of *auto-calibration signal (ACS)* lines, ensuring a fully-sampled central k-space region.

Figure 5 illustrates the k-space trajectories for both random and uniform (equidistant) subsampling at acceleration factors of 4, 6, and 8. Notably, different patterns were used for training and testing. During training, our restoration prior was only trained on a uniform mask with a subsampling rate of 6. However, for inference, we employed both uniform and random masks at subsampling rates of 4 and 6, creating a mismatch between the pre-trained restoration prior and the test configurations.

Figure: Figure 5: Illustration of the undersampling masks used for CS-MRI in this work. (a) The eight different $8\times$ uniform masks used for training the restoration prior. (b) The inference setting for ShaRP, demonstrating how the prior trained on the masks in (a) can be applied to solve other problems without retraining.
Refer to caption: x5.png

Figure: Algorithm 2 Supervised Training of CS-MRI Restoration Network

#### B.1.1 Implementation of Supervised Prior for CS-MRI

Models training for supervised case.
We use the same U-Net architecture as employed in the official implementation of DDS(^2^22[https://github.com/HJ-harry/DDS](https://github.com/HJ-harry/DDS)) for ${\mathsf{R}}(\cdot;\alpha)$. For the supervised learning case, we select 1,000 different $\alpha$ values to train the model, following the $\alpha$ schedule outlined by I^2SB (Liu et al., ). The model is trained with Adam optimizer with a learning rate of $5\times 10^{-5}$. As shown in Algorithm 2, we train our supervised learning model using eight different masks for $8\times$ uniform sampling CS-MRI reconstruction. In the pseudocode, $\{{\bm{M}}_{1},{\bm{M}}_{2},\cdots,{\bm{M}}_{8}\}$ represent the eight different MRI degradation operators, each defined by a unique sampling pattern, as shown in Figure 5 (a). This results in a total of 8,000 possible combinations of $\alpha$ values and sampling masks, effectively creating an ensemble of restoration priors during training.

Inference with a Subset of the Ensemble (Supervised Case). During inference, to simplify computation and focus on the most effective priors, we use only a subset of the supervised trained ensemble. Specifically, we fix the $\alpha$ value to a particular choice (e.g., $\alpha=0.5$) and use the 8 different sampling masks $\{{\bm{M}}_{1},{\bm{M}}_{2},\cdots,{\bm{M}}_{8}\}$, resulting in 8 restoration priors.

Step size and regularization parameter. To ensure fairness, for each problem setting, each method—both proposed and baseline—is fine-tuned for optimal PSNR using 10 slices from a validation set separate from the test set. The same step size $\gamma$ and regularization parameter $\tau$ are then applied consistently across the entire test set.

Baseline details. We compare ShaRP with several variants of denoiser- and diffusion model-based methods. For denoiser-based approaches, we include PnP-FISTA (Kamilov et al., ), PnP-ADMM (Chan et al., ). PnP-FISTA and PnP-ADMM correspond to the FISTA and ADMM variants of PnP, both utilizing AWGN denoisers built on DRUNet (Zhang et al., ). For diffusion model-based methods, we compare with DPS (Chung et al., ) and DDS (Chung et al., ), which use pre-trained diffusion models as priors and apply different posterior sampling strategies to address general inverse problems. We use the same pre-trained diffusion model configuration as outlined in the DDS paper. For all baseline methods, we fine-tuned their parameters to maximize the PSNR value. Notably, both the DRUNet denoiser and the diffusion model were trained using the same dataset employed for training our restoration prior. For a fair comparison, the diffusion model pre-trained for DDS and DPS use the same network architecture as our restoration network . All models are trained from scratch on the fastMRI training set, following the architecture settings provided in DDS(^3^33[https://github.com/HJ-harry/DDS](https://github.com/HJ-harry/DDS)). We also compared with method that also use the deep restoration prior to solve general inverse problem: DRP (Hu et al., ). For DRP, we utilize the same pre-trained restoration network as in ShaRP. However, instead of employing a set of degradation priors, DRP uses a single fixed prior. For a fair comparison, we selected the optimal fixed prior—defined by a fixed $\alpha$ and subsampling mask—based on PSNR performance on the validation set, and applied it accordingly.

#### B.1.2 Implementation of Self-Supervised Prior for CS-MRI

Figure: Algorithm 3 Self-Supervised Training of CS-MRI Restoration Network

Models training for (Self-Supervised Case).
For self-supervised training, the ground truth reference ${\bm{x}}$ is not available as a label. Instead, as shown in Algorithm 3, we work with pairs of subsampled measurements, $y_{i}$ and $y_{j}$, along with their corresponding sampling operators, ${\bm{M}}_{i}$ and ${\bm{M}}_{j}$. These paired measurements exhibit significant overlap within the shared *auto-calibration signal (ACS)* region, which increases the weighting of these overlapping k-space regions. Following the approach proposed by SSDEQ (Gan et al., ), we introduce a diagonal weighting matrix ${\bm{W}}$ to account for the oversampled regions in the loss function. By incorporating this weighted loss, we are able to train our MMSE restoration operator using incomplete measurements alone. Furthermore, unlike the supervised case where we use the combination of $\alpha$ values to form an ensemble, in the self-supervised setting, we construct the ensemble using only eight different sampling masks across the entire dataset.

Inference Using All Restoration Priors (Self-Supervised Case). During inference in the self-supervised setting, we utilize all 8 restoration priors corresponding to the different sampling masks. By incorporating the entire ensemble, we fully leverage its capacity to remove the artifacts and enhance reconstruction performance.

Step size and regularization parameter. To ensure fairness, for each problem setting, each method—both proposed and baseline—is fine-tuned for optimal PSNR using 10 slices from a validation set separate from the test set. The same step size $\gamma$ and regularization parameter $\tau$ are then applied consistently across the entire test set.

Baseline details. In the self-supervised setting, we compared ShaRP with two widely used traditional methods: TV (Block et al., ) and GRAPPA (Griswold et al., ), both of which address the restoration problem without requiring fully-sampled references. Additionally, we included SPICER (Hu et al., ), a recent state-of-the-art self-supervised deep unrolling method designed for MRI reconstruction using only pairs of undersampled measurements. To ensure consistency, we trained the SPICER model on the same amount of paired data used for training our restoration prior in the $8\times$ uniform CS-MRI setting and applied it to other CS-MRI configurations.

### B.2 Implementation details of SISR tasks

Figure: Algorithm 4 Gaussian Deblurring Restoration network training

Restoration Model training.
We use the same U-Net architecture as the Gaussian deblurring model provided by I^2SB(^4^44[https://github.com/NVlabs/I2SB](https://github.com/NVlabs/I2SB)). Utilizing the pre-trained checkpoints from their repository, we fine-tune our model accordingly. Specifically, we align with their codebase and configure the model type to OT-ODE to satisfy our MMSE restoration operator assumption.

To create an ensemble of restoration priors, we consider a family of degradation operators that are convex combinations of the identity mapping ${\mathbf{I}}$ and the Gaussian blur operator ${\bm{K}}$. The blurring operator ${\bm{K}}$ corresponds to convolution with a Gaussian blur kernel of size $31\times 31$ and standard deviation 3. Specifically, we define the degradation operator as ${\mathbf{H}}_{\alpha}=(1-\alpha){\mathsf{I}}+\alpha{\bm{K}}$, where $\alpha\in[0,1]$ controls the degradation level. By varying $\alpha$, we generate multiple degradation operators, allowing us to train the restoration network ${\mathsf{R}}$ to handle all these operators, expressed as ${\mathsf{R}}({\bm{s}},{\mathbf{H}}_{\alpha})=\mathbb{E}\left[{\bm{x}}|{\bm{s}}
,{\mathbf{H}}_{\alpha}\right]$, where ${\bm{s}}$ is the degraded image and ${\bm{x}}$ is the original image.

We select 1,000 different $\alpha$ values from the interval $[0,1]$, following the $\alpha$ schedule outlined by I^2SB (Liu et al., ). This results in 1,000 different degradation operators ${\mathbf{H}}_{\alpha}$, effectively creating an ensemble of restoration priors during training. The model is trained using the Adam optimizer with a learning rate of $5\times 10^{-5}$.

Inference with a Subset of the Ensemble. During inference, to simplify computation and focus on the most effective priors, we use only a subset of the supervised trained ensemble. Specifically, we select 6 $\alpha$ values, resulting in 6 restoration priors.

Step size and regularization parameter. To ensure fairness, for each problem setting, each method—both proposed and baseline—is fine-tuned for optimal PSNR using 5 images from a validation set separate from the test set. The same step size $\gamma$ and regularization parameter $\tau$ are then applied consistently across the entire test set.

Baseline details. We compare ShaRP against several denoiser- and diffusion model-based methods. For denoiser-based approaches, we evaluate DPIR (Zhang et al., ), which relies on half-quadratic splitting (HQS) iterations with DRUNet denoisers. For diffusion model-based methods, we compare with DPS (Chung et al., ), DDNM (Wang et al., ), and DiffPIR (Zhu et al., ). These methods all use the same pre-trained diffusion models as priors, but each employs a distinct posterior sampling strategy to solve general inverse problems. We specifically use the pre-trained diffusion model from DiffPIR. We also compared with method that also use the deep restoration prior to solve general inverse problem: DRP (Hu et al., ). For DRP, we utilize the same pre-trained deblurring network as in ShaRP. However, instead of employing a set of degradation priors, DRP uses a single fixed prior. For a fair comparison, we selected the optimal fixed prior—defined by a fixed $\alpha$ based on PSNR performance on the validation set, and applied it accordingly. For all baselines, we fine-tuned their parameters to maximize PSNR performance. Notably, the diffusion model backbone for all diffusion-based baselines was trained on the same dataset used to train our restoration prior.

## Appendix C Additional results for CS-MRI

### C.1 Performance of ShaRP for random subsampling setting

Due to space constraints, we present only the quantitative performance for the uniform subsampling setting in the main paper. In this section, we further evaluate ShaRP’s performance on random subsampling setting, with two sub-sampling rates (4$\times$ and 6$\times$), and three noise levels ($\sigma=0.005$, $0.01$, and $0.015$).

Table 4 provides a quantitative comparison of reconstruction performance across different acceleration factors and noise levels using a uniform sub-sampling mask. In all configurations, ShaRP consistently outperforms the baseline methods. The use of a set of restoration operators clearly enhances ShaRP’s performance, highlighting the effectiveness of employing multiple operators to maximize the regularization information provided by the restoration model. Figure 6 presents visual reconstructions for two test scenarios, where ShaRP accurately recovers fine brain details, particularly in the zoomed-in regions, while baseline methods tend to oversmooth or introduce artifacts. These results highlight ShaRP’s superior ability to manage structured artifacts and preserve fine details, outperforming both denoiser-based and diffusion model-based methods.

Figure: Figure 6: Visual comparison of ShaRP with baseline methods on CS-MRI for $6\times$ random sampling mask with noise $\sigma=0.015$. PSNR and SSIM values are in the top-left corner of each image. Error maps and zoomed-in areas highlight differences. Notably, ShaRP with stochastic priors outperforms state-of-the-art methods using denoiser and diffusion model priors.
Refer to caption: x6.png

**Table 4: Quantitative comparison of ShaRP with several baselines for CS-MRI using random masks at undersampling rates of 4 and 6 on fastMRI dataset. The best and second best results are highlighted. Notably, ShaRP outperforms SOTA methods based on denoisers and diffusion models.**
|  | 4$\times$ Random | 6$\times$ Random |  |  |  |  |  |  |  |  |  |  |
| --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- |
| Noise level | $\sigma$ = 0.005 | $\sigma$ = 0.010 | $\sigma$ = 0.015 | $\sigma$ = 0.005 | $\sigma$ = 0.010 | $\sigma$ = 0.015 |  |  |  |  |  |  |
| Metrics | PSNR | SSIM | PSNR | SSIM | PSNR | SSIM | PSNR | SSIM | PSNR | SSIM | PSNR | SSIM |
| Zero-filled | 25.83 | 0.815 | 25.81 | 0.812 | 25.76 | 0.807 | 22.68 | 0.724 | 22.67 | 0.722 | 22.67 | 0.719 |
| TV | 28.14 | 0.866 | 28.06 | 0.863 | 27.96 | 0.859 | 24.55 | 0.782 | 24.33 | 0.750 | 24.28 | 0.736 |
| PnP-FISTA | 29.31 | 0.863 | 28.40 | 0.817 | 27.49 | 0.799 | 26.01 | 0.797 | 25.63 | 0.756 | 24.94 | 0.717 |
| PnP-ADMM | 28.83 | 0.842 | 28.39 | 0.816 | 27.70 | 0.786 | 25.59 | 0.776 | 25.19 | 0.740 | 24.93 | 0.728 |
| DRP | 29.97 | 0.880 | 29.37 | 0.839 | 28.31 | 0.794 | 26.98 | 0.866 | 26.78 | 0.853 | 26.49 | 0.821 |
| DPS | 31.72 | 0.874 | 30.45 | 0.857 | 29.50 | 0.843 | 30.32 | 0.856 | 29.36 | 0.824 | 27.99 | 0.810 |
| DDS | 32.41 | 0.910 | 32.37 | 0.906 | 32.25 | 0.901 | 30.59 | 0.876 | 30.35 | 0.874 | 30.31 | 0.879 |
| ShaRP | 34.66 | 0.949 | 33.57 | 0.920 | 33.18 | 0.931 | 31.53 | 0.924 | 31.46 | 0.918 | 31.45 | 0.914 |

### C.2 Performance of additional baseline methods on matched and mismatched settings

In this section, we highlight an important observation: pre-trained restoration networks typically exhibit poor generalization to mismatched settings. We chose two commonly used methods (SwinIR (Liang et al., ) and E2E-VarNet (Sriram et al., )) for the specific setting of CS-MRI. We trained them on the same $8\times$ uniform subsampling setting as our restoration prior and directly applied them to solve both matched and mismatched problems, as ShaRP did. As shown in the Table 5, the baseline method’s performance drops significantly under mismatched conditions, whereas ShaRP maintains stable performance and convergence guarantees. This demonstrates ShaRP’s ability to adapt pre-trained restoration models as priors and use it to solve problems under mismatched settings. As shown in the Figure 7, due to the mismatched settings, the two baseline methods suffer from over-smoothing, lack important details, and exhibit artifacts, whereas ShaRP still provides high-quality reconstruction performance. This indicates that ShaRP can balance data fidelity and the artifact removal capabilities of the prior model, leading to an artifact-free reconstruction that preserves important details.

**Table 5: Quantitative comparison of ShaRP with task-specific baselines trained on the $8\times$ uniform mask. Baselines perform well in matched settings (highlighted in the table) but show a significant drop under mismatched conditions. In contrast, ShaRP remains robust, handling both matched and mismatched scenarios effectively.**
| Settings | $4\times$ Uniform | $4\times$ Random | $6\times$ Uniform | $6\times$ Random | $8\times$ Uniform |  |  |  |  |  |
| --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- |
| Metrics | PSNR | SSIM | PSNR | SSIM | PSNR | SSIM | PSNR | SSIM | PSNR | SSIM |
| SwinIR | 24.78 | 0.849 | 25.09 | 0.841 | 29.55 | 0.907 | 27.98 | 0.819 | 29.37 | 0.898 |
| E2E-VarNet | 35.40 | 0.957 | 33.48 | 0.945 | 32.79 | 0.936 | 31.02 | 0.913 | 32.59 | 0.919 |
| ShaRP | 37.59 | 0.963 | 34.66 | 0.949 | 33.42 | 0.940 | 31.53 | 0.924 | 32.37 | 0.907 |

Figure: Figure 7: Visual comparison of ShaRP with task-specific baseline methods on CS-MRI for $6\times$ random sampling mask with noise $\sigma=0.015$. PSNR and SSIM values are in the top-left corner of each image. Error maps and zoomed-in areas highlight differences. Notably, ShaRP with stochastic priors outperforms state-of-the-art methods using denoiser and diffusion model priors.
Refer to caption: x7.png

## Appendix D Additional visual results for SISR

In this section, we present additional visual results to numerical comparisons for the SISR task.

### D.1 Additional visual results against baselines

As illustrated in Figure 8 and Figure 9, ShaRP outperforms all baseline approaches under both blur kernel settings, achieving higher PSNR and SSIM values. Moreover, we maintain superior data consistency with the measurements while achieving enhanced perceptual quality. The use of an ensemble of deblurring priors enables our method to recover fine details at varying corruption levels, contributing to the improved performance.

Figure: Figure 8: Visual comparison of ShaRP with several well-known methods on $2\times$ SISR with gaussian blur kernel with $\sigma=1.5$. The quantities in the top-left corner of each image provide PSNR and SSIM values for each method. The squares at the bottom of each image visualize the zoomed area in the image.
Refer to caption: x8.png

Figure: Figure 9: Visual comparison of ShaRP with several well-known methods on $2\times$ SISR with gaussian blur kernel with $\sigma=1.5$. The quantities in the top-left corner of each image provide PSNR and SSIM values for each method. The squares at the bottom of each image visualize the zoomed area in the image.
Refer to caption: x9.png

### D.2 Additional visual results against DRP

To further emphasize the necessity and advantages of using an ensemble of deblurring priors, as opposed to a fixed prior like in DRP (Hu et al., ), we provide additional visual comparison results. As shown in Figure 10, ShaRP consistently recovers finer details, resulting in improved PSNR and SSIM scores, along with enhanced perceptual performance.

Figure: Figure 10: Visual comparison of ShaRP with DRP on $2\times$ SISR with gaussian blur kernel with $\sigma=1.5$. The quantities in the bottom-left corner of each image provide PSNR and SSIM values for each method. The squares at the bottom of each image visualize the zoomed area in the image.
Refer to caption: x10.png

## Peer Reviews

confidence: 3
contribution: 1
initial_rating: 3
presentation: 3
questions: Please take a look at the comments under weaknesses
soundness: 3
strengths: 1. The experiments show that training a deep restoration network on multiple tasks can improve its regularization capability.

2. Theoretical analysis is also provided.
summary: This paper extends the work of [Hu et al., 2024c] for the regularization of inverse problems using a deep restoration prior trained on multiple restoration tasks. The authors theoretically analyze the induced objective function and the convergence of the proposed algorithm. Empirical results demonstrate that the framework outperforms Gaussian denoiser priors in image reconstruction tasks, including MRI and super-resolution.
weaknesses: 1. The paper seems limited in novelty and closely resembles [Hu et al., 2024c]. The main contribution is the training of a restoration prior across the same task with varying levels of ill-posedness. The authors should clearly explain how the proposed methodology differs from [Hu et al., 2024c].

2. The authors claim that the proposed prior trained on a general restoration task outperforms Gaussian denoisers; however, there is a lack of sufficient experiments to support this claim. Gaussian denoisers can be used for solving general inverse problems, but the experiments presented here are confined to the same inverse problem and forward operator on which the prior was trained, although with different levels of ill-posedness. For example, in Section 5.1, the restoration prior is trained on MRI with an x8 acceleration and applied to the same problem with x4 and x6 acceleration factors. Although Algorithm 1 suggests that the proposed prior can be used to solve any other inverse problems, the experiments are limited to the same problem type. To validate the claimed superiority over Gaussian denoisers, the authors are requested to include experiments where the trained prior is applied to a different inverse problem from the one it was originally trained on. For example, training the prior on super-resolution and using it for solving MRI reconstruction. You can also train your prior on multiple distinct tasks as suggested by Algorithm I, for example training the prior on super-resolution and image in-painting tasks, and use it for solving a different task like MRI. This comparison is crucial as Gaussian denoisers can solve any other different inverse problems once trained.

Hu, Yuyang, et al. "A Restoration Network as an Implicit Prior." ICLR 2024.

---

confidence: 3
contribution: 2
initial_rating: 5
presentation: 3
questions: 1. Page 1: “...Tweedie’s formula (Robbins, 1956; Efron, 2011) seemingly implies that Gaussian denoising alone might be sufficient for learning priors,...”: That is indeed true and this work does not disprove this fact. 

2. Page 2: “...ShaRP provides a richer and more flexible representation of image priors…”: At this stage, it is rather vague as to what “richer and more flexible” means. 
3. Page 2: “Unlike Gaussian denoisers, the restoration models in ShaRP can often be directly
trained in a self-supervised manner”: Even Gaussian denoisers can be trained in an unsupervised manner using only noisy images (e.g., using a SURE loss). 

4. It might be good to put some of the math background in the intro, effectively shortening the material in the first three pages. 

5. Page 3: “We introduce a novel regularization concept for inverse problems that encourages solutions that produce degraded versions closely resembling real degraded images.”: Firstly, this would be the case if the degradation is the same for which the restoration network is trained, not the forward operator “A” for which you want to solve the inverse problem. Secondly, can one not promote this property using a simple Gaussian denoiser together with a data-consistency loss? 

6. Theorem 1: “Assume that the prior density $p_x$ is non-degenerate over $\mathbb{R}^n$”: This is almost always not true. 

7. Assumption 3: Define what “$b(x)$” is. 

Requested changes:

1. Provide a comprehensive review of the convergence analysis of biased SGD in terms of the assumptions and results to put your theoretical contributions in perspective. 

2. Clarify whether one can approximate the conditional expectation $E[x|s, H]$ using a restoration operator trained on a self-supervised loss. If not, the interpretation of ShaRP as biased SGD does not hold in this case.
soundness: 2
strengths: 1. Using general restoration operators instead of simple Gaussian denoisers is interesting and novel. Theorem 1 also puts this idea on a solid theoretical foundation. 

2. The experimental results are compelling. Comparison of ShaRP with different baseline methods, especially PnP algorithms that utilize Gaussian denoisers, provides evidence for the claim made by the authors in the abstract.
summary: The paper proposes a new approach for learning a regularizer for linear imaging inverse problems. While Gaussian denoisers have been successfully used as effective image priors, the authors extend the idea by using more general restoration operators as image priors. In the numerical experiments, the authors show the reconstruction performance with their new “ShaRP” regularizers for two important inverse problems (compressive MRI and single-image super-resolution). The numerical experiments are performed with ShaRP regularizers trained in supervised and self-supervised manners.
weaknesses: 1. I am not particularly convinced about the motivation for using a general restoration prior, that arises out of training a deep reconstruction operator for linear inverse problems (with several degradation operators), for solving another linear inverse problem. While learning Gaussian denoisers is a problem that is fundamentally easier than solving a general inverse problem (with an ill-posed operator), this is not the case here. 

2. The novelty in the convergence analysis (namely, the proof of Theorem 2) is unclear to me. To my understanding, the assumptions and the result are along the same lines as in the paper entitled “A Guide Through the Zoo of Biased SGD” by Demidovich et al. (which the authors have cited). 

3. Although the authors show the performance of shaRP trained in an unsupervised manner, I don’t think the training loss used in the unsupervised case (e.g., the loss in Algorithm 3) corresponds to a restoration operator that approximates the conditional expectation $E[x|s, H]$. Consequently, the theoretical analysis does not explain the experimental results with a restoration operator trained in an unsupervised manner.

---

confidence: 5
contribution: 3
initial_rating: 6
presentation: 3
questions: 1. Can you please provide more details about the impact of the choice of $b$ on the numerical results? An ablation study might help here. Likewise, you could numerically evaluate the impact of $\alpha$.
2. What happens if you modify the number of restoration priors in B.1.1?
3. Please rewrite the motivation in Section 4. In particular, I am missing a good motivation and reasoning for the actual definition in equation (6), which might help the reader to better understand ShaRP. In addition, the object $G_\sigma(s-Hx)$ could be better motivated from a mathematical point of view.
4. The results in Table 4 do not show a clear tendency for DiffPIR, DRP, or ShaRP. Is there a reason for this available? Are there specific conditions under which each method performs best?
soundness: 4
strengths: The paper introduces a novel concept for stochastic priors, which is complemented by two important theoretical assertions (see summary). The assumptions can be regarded as mild. Both results are highly relevant for numerical experiments. 
The paper remarkably advances SOTA for CS-MRI and is on par with competing methods in the case of SISR.
summary: This paper introduces a novel concept termed ShaRP for stochastic priors for the regularization of inverse problems. At its core, a set of $b$ restoration problems is considered, for which the MMSE is minimized. In this setting, ShaRP is the expectation of the probability of the degraded version of the images given the probability density of the observations. Under certain assumptions, a closed-form solution of the gradient of the regularizer is derived and the convergence of $\nabla f$ is shown. The results are numerically tested for CS-MRI and SISR.
weaknesses: The presentation of the paper could be improved in some places. In particular, I am missing a more concise introduction to the motivation behind the actual concept of ShaRP in Section 4.
The parameter $b$ representing the number of restoration problems is essential for the numerical performance. However, I am lacking a substantial discussion of the role of this parameter. In particular, no details about the role of $\alpha$ (see Section 5.1) and its impact on the results were provided.

---

confidence: 5
contribution: 1
initial_rating: 1
presentation: 3
questions: I don't have any questions for the authors.
soundness: 3
strengths: The paper is clearly written, for the most part.
summary: This paper extends the "deep restoration priors" approach from Hu et al 2024c from one degradation operator to multiple ones.
weaknesses: This paper is a trivial extension of Hu et al 2024c.  In this paper, the degradation operator is randomly chosen from a set, whereas in Hu et al 2024c there was a single operator.  But this change led to no new challenges or questions.  For example, the "ShaRP" regularizer in (6) is a trivial extension of (9) in Hu et al 2024c that now takes an expectation over H.  Figure 2 in this paper is a direct copy of Figure 1 from Hu et al 2024c.  The experiments are conducted on different linear inverse problems, but there is no intellectual contribution there.

---

confidence: 5
contribution: 2
initial_rating: 6
presentation: 3
questions: What is the practical inference time of the proposed method in comparison to state-of-the-art (SOTA) methods? Additionally, the visual comparisons presented do not clearly demonstrate significant improvements. It would be beneficial to include more compelling visual examples to better illustrate the advantages of the proposed approach.
soundness: 3
strengths: - The paper is well-written and easy to follow.
- The theoretical analysis of the convergence is thorough and well-explained.
summary: This paper introduces Stochastic deep Restoration Priors (ShaRP), a method that leverages an ensemble of deep restoration models to regularize imaging inverse problems. ShaRP improves upon Gaussian denoiser-based methods by handling structured artifacts more effectively, enabling self-supervised training without fully sampled data.
weaknesses: 1. **Comparison to Diffusion-Based Methods:**
   The paper lacks a comprehensive comparison to diffusion-based methods, such as DiffIR [1] and DDRM [2]. Including these comparisons would strengthen the results and provide clarity on how ShaRP performs relative to other leading methods in the field.

2. **Supervised vs. Self-Supervised ShaRP:**
   In line 154, the authors mention a key contribution: "We implement ShaRP with both supervised and self-supervised restoration models as priors and test it on two inverse problems: CS-MRI and SISR." However, the experimental section does not provide a direct comparison between the self-supervised and supervised versions of ShaRP for the same task. This comparison is necessary to assess the benefits of the self-supervised approach.

3. **Self-Supervised Nature:**
   For a restoration network to be trained on a set of tasks, such as a set of blur kernels \(H_i\), access to ground truth data is still required. This approach, which involves sampling multiple times from the ground truth, raises questions about whether the method can truly be considered self-supervised. Clarification is needed regarding the self-supervised claim.

4. **Use of Multiple Degradation Operators:**
   The rationale for using a set of degradation operators \(H_1, H_2, \ldots, H_k\) in cases where the target problem involves only a single fixed operator (e.g., \(H_1\)) is unclear. It would be helpful if the authors could explain why introducing multiple degradation operators is necessary or beneficial when solving a fixed-task problem.

**References:**

- [1] DiffIR: Efficient Diffusion Model for Image Restoration
- [2] Denoising Diffusion Restoration Models

Meta-review: This paper addresses deep learning reconstruction for computational imaging. It builds on the work on deep restoration priors by Hu et al. from 2024 by introducing multiple degradation operator in place of one, and using a different optimization algorithm which is an adaptation of SNORE (Renaud et al 2024b). This new combination results in improved numerical results, in particular on subsampled MRI.

Reviewers made positive remarks about the clarity and the numerical results, but there has been a broad agreement  that the novelty in this paper relative to Hu et al. 24 and other prior work is limited and not sufficient for acceptance at ICLR. While the rebuttal phase has been intense, no fundamentally new information emerged from the discussions, while the reviewers made precise arguments to support their position.

The manuscript may benefit from a "full-disclosure" rewrite that proactively addresses the criticism of N69Z, XTiX, and E3Gi. A focus on domain-specific numerics (MRI, medical imaging in general) and practical impact may further strengthen the paper. The reviewers provided precise, detailed remarks that should help guide these revisions effectively.

While the leading-order criticisms were about novelty,  part of the dialogue between the authors and the reviewers revolved around the meaning of "prior" and its relation to denoisers via Tweedie's formula. A prior distribution doesn't depend on a specific observation process. Priors learned by Gaussian denoisers (via Tweedie) and "priors" that result from restoration networks which, as the authors explain, are tuned to a particular problem, are different objects. The authors are justified in stating that "In practice, the effectiveness of a denoiser is tied to how well the noise model it's trained on matches the actual noise in the image". But this is a statement about structured denoising more than about priors. A bona fide prior and a known forward model allow to compute any estimator, at least in principle. They can be used to generate a dataset and do any kind of ML. The question is how to build effective algorithms that exploit this (with finite samples) and at this point one may obtain improved results by schemes such as the one proposed by the authors. This distinction might benefit from further clarification in the revised manuscript.

For these reasons I am recommending a rejection. Nonetheless, I hope that the authors will find the reviewers' feedback and the discussions valuable as they work on refining their paper.

Decision: Reject